Nektar++
StdHexExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdHexExp.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: Heaxhedral methods
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <StdRegions/StdHexExp.h>
36 
37 #ifdef max
38 #undef max
39 #endif
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 namespace StdRegions
46 {
47 
48 StdHexExp::StdHexExp()
49 {
50 }
51 
52 StdHexExp::StdHexExp(const LibUtilities::BasisKey &Ba,
53  const LibUtilities::BasisKey &Bb,
54  const LibUtilities::BasisKey &Bc)
55  : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), 3,
56  Ba, Bb, Bc),
57  StdExpansion3D(Ba.GetNumModes() * Bb.GetNumModes() * Bc.GetNumModes(), Ba,
58  Bb, Bc)
59 {
60 }
61 
62 StdHexExp::StdHexExp(const StdHexExp &T) : StdExpansion(T), StdExpansion3D(T)
63 {
64 }
65 
66 StdHexExp::~StdHexExp()
67 {
68 }
69 
70 bool StdHexExp::v_IsBoundaryInteriorExpansion()
71 {
72  return (m_base[0]->GetBasisType() == LibUtilities::eModified_A &&
73  m_base[1]->GetBasisType() == LibUtilities::eModified_A &&
74  m_base[2]->GetBasisType() == LibUtilities::eModified_A) ||
75  (m_base[0]->GetBasisType() == LibUtilities::eGLL_Lagrange &&
76  m_base[1]->GetBasisType() == LibUtilities::eGLL_Lagrange &&
77  m_base[1]->GetBasisType() == LibUtilities::eGLL_Lagrange);
78 }
79 
80 ///////////////////////////////
81 /// Differentiation Methods
82 ///////////////////////////////
83 /**
84  * For Hexahedral region can use the PhysTensorDeriv function defined
85  * under StdExpansion. Following tenserproduct:
86  */
87 void StdHexExp::v_PhysDeriv(const Array<OneD, const NekDouble> &inarray,
88  Array<OneD, NekDouble> &out_d0,
89  Array<OneD, NekDouble> &out_d1,
90  Array<OneD, NekDouble> &out_d2)
91 {
92  StdExpansion3D::PhysTensorDeriv(inarray, out_d0, out_d1, out_d2);
93 }
94 
95 /**
96  * @param dir Direction in which to compute derivative.
97  * Valid values are 0, 1, 2.
98  * @param inarray Input array.
99  * @param outarray Output array.
100  */
101 void StdHexExp::v_PhysDeriv(const int dir,
102  const Array<OneD, const NekDouble> &inarray,
103  Array<OneD, NekDouble> &outarray)
104 {
105  switch (dir)
106  {
107  case 0:
108  {
109  PhysDeriv(inarray, outarray, NullNekDouble1DArray,
110  NullNekDouble1DArray);
111  }
112  break;
113  case 1:
114  {
115  PhysDeriv(inarray, NullNekDouble1DArray, outarray,
116  NullNekDouble1DArray);
117  }
118  break;
119  case 2:
120  {
121  PhysDeriv(inarray, NullNekDouble1DArray, NullNekDouble1DArray,
122  outarray);
123  }
124  break;
125  default:
126  {
127  ASSERTL1(false, "input dir is out of range");
128  }
129  break;
130  }
131 }
132 
133 void StdHexExp::v_StdPhysDeriv(const Array<OneD, const NekDouble> &inarray,
134  Array<OneD, NekDouble> &out_d0,
135  Array<OneD, NekDouble> &out_d1,
136  Array<OneD, NekDouble> &out_d2)
137 {
138  StdHexExp::v_PhysDeriv(inarray, out_d0, out_d1, out_d2);
139 }
140 
141 void StdHexExp::v_StdPhysDeriv(const int dir,
142  const Array<OneD, const NekDouble> &inarray,
143  Array<OneD, NekDouble> &outarray)
144 {
145  StdHexExp::v_PhysDeriv(dir, inarray, outarray);
146 }
147 
148 /**
149  * Backward transformation is three dimensional tensorial expansion
150  * \f$ u (\xi_{1i}, \xi_{2j}, \xi_{3k})
151  * = \sum_{p=0}^{Q_x} \psi_p^a (\xi_{1i})
152  * \lbrace { \sum_{q=0}^{Q_y} \psi_{q}^a (\xi_{2j})
153  * \lbrace { \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{r}^a (\xi_{3k})
154  * \rbrace}
155  * \rbrace}. \f$
156  * And sumfactorizing step of the form is as:\\
157  * \f$ f_{r} (\xi_{3k})
158  * = \sum_{r=0}^{Q_z} \hat u_{pqr} \psi_{r}^a (\xi_{3k}),\\
159  * g_{p} (\xi_{2j}, \xi_{3k})
160  * = \sum_{r=0}^{Q_y} \psi_{p}^a (\xi_{2j}) f_{r} (\xi_{3k}),\\
161  * u(\xi_{1i}, \xi_{2j}, \xi_{3k})
162  * = \sum_{p=0}^{Q_x} \psi_{p}^a (\xi_{1i}) g_{p} (\xi_{2j}, \xi_{3k}).
163  * \f$
164  *
165  * @param inarray ?
166  * @param outarray ?
167  */
168 void StdHexExp::v_BwdTrans(const Array<OneD, const NekDouble> &inarray,
169  Array<OneD, NekDouble> &outarray)
170 {
171  ASSERTL1((m_base[1]->GetBasisType() != LibUtilities::eOrtho_B) ||
172  (m_base[1]->GetBasisType() != LibUtilities::eModified_B),
173  "Basis[1] is not a general tensor type");
174 
175  ASSERTL1((m_base[2]->GetBasisType() != LibUtilities::eOrtho_C) ||
176  (m_base[2]->GetBasisType() != LibUtilities::eModified_C),
177  "Basis[2] is not a general tensor type");
178 
179  if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
180  m_base[2]->Collocation())
181  {
182  Vmath::Vcopy(m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints() *
183  m_base[2]->GetNumPoints(),
184  inarray, 1, outarray, 1);
185  }
186  else
187  {
188  StdHexExp::BwdTrans_SumFac(inarray, outarray);
189  }
190 }
191 
192 /**
193  *
194  */
195 void StdHexExp::v_BwdTrans_SumFac(const Array<OneD, const NekDouble> &inarray,
196  Array<OneD, NekDouble> &outarray)
197 {
198  Array<OneD, NekDouble> wsp(
199  m_base[0]->GetNumPoints() * m_base[2]->GetNumModes() *
200  (m_base[1]->GetNumModes() + m_base[1]->GetNumPoints())); // FIX THIS
201 
202  BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
203  m_base[2]->GetBdata(), inarray, outarray, wsp, true,
204  true, true);
205 }
206 
207 /**
208  * @param base0 x-dirn basis matrix
209  * @param base1 y-dirn basis matrix
210  * @param base2 z-dirn basis matrix
211  * @param inarray Input vector of modes.
212  * @param outarray Output vector of physical space data.
213  * @param wsp Workspace of size Q_x*P_z*(P_y+Q_y)
214  * @param doCheckCollDir0 Check for collocation of basis.
215  * @param doCheckCollDir1 Check for collocation of basis.
216  * @param doCheckCollDir2 Check for collocation of basis.
217  * @todo Account for some directions being collocated. See
218  * StdQuadExp as an example.
219  */
220 void StdHexExp::v_BwdTrans_SumFacKernel(
221  const Array<OneD, const NekDouble> &base0,
222  const Array<OneD, const NekDouble> &base1,
223  const Array<OneD, const NekDouble> &base2,
224  const Array<OneD, const NekDouble> &inarray,
225  Array<OneD, NekDouble> &outarray, Array<OneD, NekDouble> &wsp,
226  bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
227 {
228  int nquad0 = m_base[0]->GetNumPoints();
229  int nquad1 = m_base[1]->GetNumPoints();
230  int nquad2 = m_base[2]->GetNumPoints();
231  int nmodes0 = m_base[0]->GetNumModes();
232  int nmodes1 = m_base[1]->GetNumModes();
233  int nmodes2 = m_base[2]->GetNumModes();
234 
235  // Check if using collocation, if requested.
236  bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
237  bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
238  bool colldir2 = doCheckCollDir2 ? (m_base[2]->Collocation()) : false;
239 
240  // If collocation in all directions, Physical values at quadrature
241  // points is just a copy of the modes.
242  if (colldir0 && colldir1 && colldir2)
243  {
244  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, outarray.get(), 1);
245  }
246  else
247  {
248  // Check sufficiently large workspace.
249  ASSERTL1(wsp.size() >= nquad0 * nmodes2 * (nmodes1 + nquad1),
250  "Workspace size is not sufficient");
251 
252  // Assign second half of workspace for 2nd DGEMM operation.
253  Array<OneD, NekDouble> wsp2 = wsp + nquad0 * nmodes1 * nmodes2;
254 
255  // BwdTrans in each direction using DGEMM
256  Blas::Dgemm('T', 'T', nmodes1 * nmodes2, nquad0, nmodes0, 1.0,
257  &inarray[0], nmodes0, base0.get(), nquad0, 0.0, &wsp[0],
258  nmodes1 * nmodes2);
259  Blas::Dgemm('T', 'T', nquad0 * nmodes2, nquad1, nmodes1, 1.0, &wsp[0],
260  nmodes1, base1.get(), nquad1, 0.0, &wsp2[0],
261  nquad0 * nmodes2);
262  Blas::Dgemm('T', 'T', nquad0 * nquad1, nquad2, nmodes2, 1.0, &wsp2[0],
263  nmodes2, base2.get(), nquad2, 0.0, &outarray[0],
264  nquad0 * nquad1);
265  }
266 }
267 
268 /**
269  * Solves the system
270  * \f$ \mathbf{B^{\top}WB\hat{u}}=\mathbf{B^{\top}Wu^{\delta}} \f$
271  *
272  * @param inarray array of physical quadrature points to be
273  * transformed, \f$ \mathbf{u^{\delta}} \f$.
274  * @param outarray array of expansion coefficients,
275  * \f$ \mathbf{\hat{u}} \f$.
276  */
277 void StdHexExp::v_FwdTrans(const Array<OneD, const NekDouble> &inarray,
278  Array<OneD, NekDouble> &outarray)
279 {
280  // If using collocation expansion, coefficients match physical
281  // data points so just do a direct copy.
282  if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()) &&
283  (m_base[2]->Collocation()))
284  {
285  Vmath::Vcopy(GetNcoeffs(), &inarray[0], 1, &outarray[0], 1);
286  }
287  else
288  {
289  // Compute B^TWu
290  IProductWRTBase(inarray, outarray);
291 
292  // get Mass matrix inverse
293  StdMatrixKey masskey(eInvMass, DetShapeType(), *this);
294  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
295 
296  // copy inarray in case inarray == outarray
297  DNekVec in(m_ncoeffs, outarray);
298  DNekVec out(m_ncoeffs, outarray, eWrapper);
299 
300  // Solve for coefficients.
301  out = (*matsys) * in;
302  }
303 }
304 
305 /**
306  * \f$
307  * \begin{array}{rcl}
308  * I_{pqr} = (\phi_{pqr}, u)_{\delta} & = &
309  * \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2}
310  * \psi_{p}^{a}(\xi_{1i}) \psi_{q}^{a}(\xi_{2j}) \psi_{r}^{a}(\xi_{3k})
311  * w_i w_j w_k u(\xi_{1,i} \xi_{2,j} \xi_{3,k})
312  *
313  * J_{i,j,k}\\ & = & \sum_{i=0}^{nq_0} \psi_p^a(\xi_{1,i})
314  * \sum_{j=0}^{nq_1} \psi_{q}^a(\xi_{2,j})
315  * \sum_{k=0}^{nq_2} \psi_{r}^a
316  * u(\xi_{1i},\xi_{2j},\xi_{3k}) J_{i,j,k}
317  * \end{array} \f$ \n
318  * where
319  * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3)
320  * = \psi_p^a( \xi_1) \psi_{q}^a(\xi_2) \psi_{r}^a(\xi_3) \f$ \n
321  * which can be implemented as \n
322  * \f$f_{r} (\xi_{3k})
323  * = \sum_{k=0}^{nq_3} \psi_{r}^a u(\xi_{1i},\xi_{2j}, \xi_{3k})
324  * J_{i,j,k} = {\bf B_3 U} \f$ \n
325  * \f$ g_{q} (\xi_{3k})
326  * = \sum_{j=0}^{nq_1} \psi_{q}^a(\xi_{2j}) f_{r}(\xi_{3k})
327  * = {\bf B_2 F} \f$ \n
328  * \f$ (\phi_{pqr}, u)_{\delta}
329  * = \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{q} (\xi_{3k})
330  * = {\bf B_1 G} \f$
331  *
332  * @param inarray ?
333  * @param outarray ?
334  */
335 void StdHexExp::v_IProductWRTBase(const Array<OneD, const NekDouble> &inarray,
336  Array<OneD, NekDouble> &outarray)
337 {
338  if (m_base[0]->Collocation() && m_base[1]->Collocation() &&
339  m_base[2]->Collocation())
340  {
341  MultiplyByQuadratureMetric(inarray, outarray);
342  }
343  else
344  {
345  StdHexExp::v_IProductWRTBase_SumFac(inarray, outarray);
346  }
347 }
348 
349 /**
350  * Implementation of the local matrix inner product operation.
351  */
352 void StdHexExp::v_IProductWRTBase_MatOp(
353  const Array<OneD, const NekDouble> &inarray,
354  Array<OneD, NekDouble> &outarray)
355 {
356  int nq = GetTotPoints();
357  StdMatrixKey iprodmatkey(eIProductWRTBase, DetShapeType(), *this);
358  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
359 
360  Blas::Dgemv('N', m_ncoeffs, nq, 1.0, iprodmat->GetPtr().get(), m_ncoeffs,
361  inarray.get(), 1, 0.0, outarray.get(), 1);
362 }
363 
364 /**
365  * Implementation of the sum-factorization inner product operation.
366  */
367 void StdHexExp::v_IProductWRTBase_SumFac(
368  const Array<OneD, const NekDouble> &inarray,
369  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
370 {
371  int nquad0 = m_base[0]->GetNumPoints();
372  int nquad1 = m_base[1]->GetNumPoints();
373  int nquad2 = m_base[2]->GetNumPoints();
374  int order0 = m_base[0]->GetNumModes();
375  int order1 = m_base[1]->GetNumModes();
376 
377  Array<OneD, NekDouble> wsp(nquad0 * nquad1 * (nquad2 + order0) +
378  order0 * order1 * nquad2);
379 
380  if (multiplybyweights)
381  {
382  Array<OneD, NekDouble> tmp(inarray.size());
383  MultiplyByQuadratureMetric(inarray, tmp);
384 
385  StdHexExp::IProductWRTBase_SumFacKernel(
386  m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
387  tmp, outarray, wsp, true, true, true);
388  }
389  else
390  {
391  StdHexExp::IProductWRTBase_SumFacKernel(
392  m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
393  inarray, outarray, wsp, true, true, true);
394  }
395 }
396 
397 /**
398  * Implementation of the sum-factorisation inner product operation.
399  * @todo Implement cases where only some directions are collocated.
400  */
401 void StdHexExp::v_IProductWRTBase_SumFacKernel(
402  const Array<OneD, const NekDouble> &base0,
403  const Array<OneD, const NekDouble> &base1,
404  const Array<OneD, const NekDouble> &base2,
405  const Array<OneD, const NekDouble> &inarray,
406  Array<OneD, NekDouble> &outarray, Array<OneD, NekDouble> &wsp,
407  bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
408 {
409  int nquad0 = m_base[0]->GetNumPoints();
410  int nquad1 = m_base[1]->GetNumPoints();
411  int nquad2 = m_base[2]->GetNumPoints();
412  int nmodes0 = m_base[0]->GetNumModes();
413  int nmodes1 = m_base[1]->GetNumModes();
414  int nmodes2 = m_base[2]->GetNumModes();
415 
416  bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
417  bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
418  bool colldir2 = doCheckCollDir2 ? (m_base[2]->Collocation()) : false;
419 
420  if (colldir0 && colldir1 && colldir2)
421  {
422  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, outarray.get(), 1);
423  }
424  else
425  {
426  ASSERTL1(wsp.size() >= nmodes0 * nquad2 * (nquad1 + nmodes1),
427  "Insufficient workspace size");
428 
429  Array<OneD, NekDouble> tmp0 = wsp;
430  Array<OneD, NekDouble> tmp1 = wsp + nmodes0 * nquad1 * nquad2;
431 
432  if (colldir0)
433  {
434  // reshuffle data for next operation.
435  for (int n = 0; n < nmodes0; ++n)
436  {
437  Vmath::Vcopy(nquad1 * nquad2, inarray.get() + n, nquad0,
438  tmp0.get() + nquad1 * nquad2 * n, 1);
439  }
440  }
441  else
442  {
443  Blas::Dgemm('T', 'N', nquad1 * nquad2, nmodes0, nquad0, 1.0,
444  inarray.get(), nquad0, base0.get(), nquad0, 0.0,
445  tmp0.get(), nquad1 * nquad2);
446  }
447 
448  if (colldir1)
449  {
450  // reshuffle data for next operation.
451  for (int n = 0; n < nmodes1; ++n)
452  {
453  Vmath::Vcopy(nquad2 * nmodes0, tmp0.get() + n, nquad1,
454  tmp1.get() + nquad2 * nmodes0 * n, 1);
455  }
456  }
457  else
458  {
459  Blas::Dgemm('T', 'N', nquad2 * nmodes0, nmodes1, nquad1, 1.0,
460  tmp0.get(), nquad1, base1.get(), nquad1, 0.0,
461  tmp1.get(), nquad2 * nmodes0);
462  }
463 
464  if (colldir2)
465  {
466  // reshuffle data for next operation.
467  for (int n = 0; n < nmodes2; ++n)
468  {
469  Vmath::Vcopy(nmodes0 * nmodes1, tmp1.get() + n, nquad2,
470  outarray.get() + nmodes0 * nmodes1 * n, 1);
471  }
472  }
473  else
474  {
475  Blas::Dgemm('T', 'N', nmodes0 * nmodes1, nmodes2, nquad2, 1.0,
476  tmp1.get(), nquad2, base2.get(), nquad2, 0.0,
477  outarray.get(), nmodes0 * nmodes1);
478  }
479  }
480 }
481 
482 void StdHexExp::v_IProductWRTDerivBase(
483  const int dir, const Array<OneD, const NekDouble> &inarray,
484  Array<OneD, NekDouble> &outarray)
485 {
486  StdHexExp::IProductWRTDerivBase_SumFac(dir, inarray, outarray);
487 }
488 
489 void StdHexExp::v_IProductWRTDerivBase_MatOp(
490  const int dir, const Array<OneD, const NekDouble> &inarray,
491  Array<OneD, NekDouble> &outarray)
492 {
493  ASSERTL0((dir == 0) || (dir == 1) || (dir == 2),
494  "input dir is out of range");
495 
496  int nq = GetTotPoints();
497  MatrixType mtype = eIProductWRTDerivBase0;
498 
499  switch (dir)
500  {
501  case 0:
502  mtype = eIProductWRTDerivBase0;
503  break;
504  case 1:
505  mtype = eIProductWRTDerivBase1;
506  break;
507  case 2:
508  mtype = eIProductWRTDerivBase2;
509  break;
510  }
511 
512  StdMatrixKey iprodmatkey(mtype, DetShapeType(), *this);
513  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
514 
515  Blas::Dgemv('N', m_ncoeffs, nq, 1.0, iprodmat->GetPtr().get(), m_ncoeffs,
516  inarray.get(), 1, 0.0, outarray.get(), 1);
517 }
518 
519 void StdHexExp::v_IProductWRTDerivBase_SumFac(
520  const int dir, const Array<OneD, const NekDouble> &inarray,
521  Array<OneD, NekDouble> &outarray)
522 {
523  ASSERTL0((dir == 0) || (dir == 1) || (dir == 2),
524  "input dir is out of range");
525 
526  int nquad1 = m_base[1]->GetNumPoints();
527  int nquad2 = m_base[2]->GetNumPoints();
528  int order0 = m_base[0]->GetNumModes();
529  int order1 = m_base[1]->GetNumModes();
530 
531  // If outarray > inarray then no need for temporary storage.
532  Array<OneD, NekDouble> tmp = outarray;
533  if (outarray.size() < inarray.size())
534  {
535  tmp = Array<OneD, NekDouble>(inarray.size());
536  }
537 
538  // Need workspace for sumfackernel though
539  Array<OneD, NekDouble> wsp(order0 * nquad2 * (nquad1 + order1));
540 
541  // multiply by integration constants
542  MultiplyByQuadratureMetric(inarray, tmp);
543 
544  // perform sum-factorisation
545  switch (dir)
546  {
547  case 0:
548  IProductWRTBase_SumFacKernel(
549  m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
550  m_base[2]->GetBdata(), tmp, outarray, wsp, false, true, true);
551  break;
552  case 1:
553  IProductWRTBase_SumFacKernel(
554  m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
555  m_base[2]->GetBdata(), tmp, outarray, wsp, true, false, true);
556  break;
557  case 2:
558  IProductWRTBase_SumFacKernel(
559  m_base[0]->GetBdata(), m_base[1]->GetBdata(),
560  m_base[2]->GetDbdata(), tmp, outarray, wsp, true, true, false);
561  break;
562  }
563 }
564 
565 void StdHexExp::v_LocCoordToLocCollapsed(const Array<OneD, const NekDouble> &xi,
566  Array<OneD, NekDouble> &eta)
567 {
568  eta[0] = xi[0];
569  eta[1] = xi[1];
570  eta[2] = xi[2];
571 }
572 
573 void StdHexExp::v_LocCollapsedToLocCoord(
574  const Array<OneD, const NekDouble> &eta, Array<OneD, NekDouble> &xi)
575 {
576  xi[0] = eta[0];
577  xi[1] = eta[1];
578  xi[2] = eta[2];
579 }
580 
581 /**
582  * @note for hexahedral expansions _base[0] (i.e. p) modes run fastest.
583  */
584 void StdHexExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
585 {
586  int nquad0 = m_base[0]->GetNumPoints();
587  int nquad1 = m_base[1]->GetNumPoints();
588  int nquad2 = m_base[2]->GetNumPoints();
589 
590  Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
591  Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
592  Array<OneD, const NekDouble> base2 = m_base[2]->GetBdata();
593 
594  int btmp0 = m_base[0]->GetNumModes();
595  int btmp1 = m_base[1]->GetNumModes();
596  int mode2 = mode / (btmp0 * btmp1);
597  int mode1 = (mode - mode2 * btmp0 * btmp1) / btmp0;
598  int mode0 = (mode - mode2 * btmp0 * btmp1) % btmp0;
599 
600  ASSERTL2(mode == mode2 * btmp0 * btmp1 + mode1 * btmp0 + mode0,
601  "Mode lookup failed.");
602  ASSERTL2(mode < m_ncoeffs,
603  "Calling argument mode is larger than total expansion "
604  "order");
605 
606  for (int i = 0; i < nquad1 * nquad2; ++i)
607  {
608  Vmath::Vcopy(nquad0, (NekDouble *)(base0.get() + mode0 * nquad0), 1,
609  &outarray[0] + i * nquad0, 1);
610  }
611 
612  for (int j = 0; j < nquad2; ++j)
613  {
614  for (int i = 0; i < nquad0; ++i)
615  {
616  Vmath::Vmul(nquad1, (NekDouble *)(base1.get() + mode1 * nquad1), 1,
617  &outarray[0] + i + j * nquad0 * nquad1, nquad0,
618  &outarray[0] + i + j * nquad0 * nquad1, nquad0);
619  }
620  }
621 
622  for (int i = 0; i < nquad2; i++)
623  {
624  Blas::Dscal(nquad0 * nquad1, base2[mode2 * nquad2 + i],
625  &outarray[0] + i * nquad0 * nquad1, 1);
626  }
627 }
628 
629 NekDouble StdHexExp::v_PhysEvaluateBasis(
630  const Array<OneD, const NekDouble> &coords, int mode)
631 {
632  ASSERTL2(coords[0] > -1 - NekConstants::kNekZeroTol, "coord[0] < -1");
633  ASSERTL2(coords[0] < 1 + NekConstants::kNekZeroTol, "coord[0] > 1");
634  ASSERTL2(coords[1] > -1 - NekConstants::kNekZeroTol, "coord[1] < -1");
635  ASSERTL2(coords[1] < 1 + NekConstants::kNekZeroTol, "coord[1] > 1");
636  ASSERTL2(coords[2] > -1 - NekConstants::kNekZeroTol, "coord[2] < -1");
637  ASSERTL2(coords[2] < 1 + NekConstants::kNekZeroTol, "coord[2] > 1");
638 
639  const int nm0 = m_base[0]->GetNumModes();
640  const int nm1 = m_base[1]->GetNumModes();
641  const int mode2 = mode / (nm0 * nm1);
642  const int mode1 = (mode - mode2 * nm0 * nm1) / nm0;
643  const int mode0 = (mode - mode2 * nm0 * nm1) % nm0;
644 
645  return StdExpansion::BaryEvaluateBasis<0>(coords[0], mode0) *
646  StdExpansion::BaryEvaluateBasis<1>(coords[1], mode1) *
647  StdExpansion::BaryEvaluateBasis<2>(coords[2], mode2);
648 }
649 
650 int StdHexExp::v_GetNverts() const
651 {
652  return 8;
653 }
654 
655 int StdHexExp::v_GetNedges() const
656 {
657  return 12;
658 }
659 
660 int StdHexExp::v_GetNtraces() const
661 {
662  return 6;
663 }
664 
665 LibUtilities::ShapeType StdHexExp::v_DetShapeType() const
666 {
667  return LibUtilities::eHexahedron;
668 }
669 
670 int StdHexExp::v_NumBndryCoeffs() const
671 {
672  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
673  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
674  "BasisType is not a boundary interior form");
675  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
676  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
677  "BasisType is not a boundary interior form");
678  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
679  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
680  "BasisType is not a boundary interior form");
681 
682  int nmodes0 = m_base[0]->GetNumModes();
683  int nmodes1 = m_base[1]->GetNumModes();
684  int nmodes2 = m_base[2]->GetNumModes();
685 
686  return (2 * (nmodes0 * nmodes1 + nmodes0 * nmodes2 + nmodes1 * nmodes2) -
687  4 * (nmodes0 + nmodes1 + nmodes2) + 8);
688 }
689 
690 int StdHexExp::v_NumDGBndryCoeffs() const
691 {
692  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
693  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
694  "BasisType is not a boundary interior form");
695  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
696  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
697  "BasisType is not a boundary interior form");
698  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
699  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
700  "BasisType is not a boundary interior form");
701 
702  int nmodes0 = m_base[0]->GetNumModes();
703  int nmodes1 = m_base[1]->GetNumModes();
704  int nmodes2 = m_base[2]->GetNumModes();
705 
706  return 2 * (nmodes0 * nmodes1 + nmodes0 * nmodes2 + nmodes1 * nmodes2);
707 }
708 
709 int StdHexExp::v_GetTraceNcoeffs(const int i) const
710 {
711  ASSERTL2((i >= 0) && (i <= 5), "face id is out of range");
712  if ((i == 0) || (i == 5))
713  {
714  return GetBasisNumModes(0) * GetBasisNumModes(1);
715  }
716  else if ((i == 1) || (i == 3))
717  {
718  return GetBasisNumModes(0) * GetBasisNumModes(2);
719  }
720  else
721  {
722  return GetBasisNumModes(1) * GetBasisNumModes(2);
723  }
724 }
725 
726 int StdHexExp::v_GetTraceIntNcoeffs(const int i) const
727 {
728  ASSERTL2((i >= 0) && (i <= 5), "face id is out of range");
729  if ((i == 0) || (i == 5))
730  {
731  return (GetBasisNumModes(0) - 2) * (GetBasisNumModes(1) - 2);
732  }
733  else if ((i == 1) || (i == 3))
734  {
735  return (GetBasisNumModes(0) - 2) * (GetBasisNumModes(2) - 2);
736  }
737  else
738  {
739  return (GetBasisNumModes(1) - 2) * (GetBasisNumModes(2) - 2);
740  }
741 }
742 
743 int StdHexExp::v_GetTotalTraceIntNcoeffs() const
744 {
745  return 2 * ((GetBasisNumModes(0) - 2) * (GetBasisNumModes(1) - 2) +
746  (GetBasisNumModes(0) - 2) * (GetBasisNumModes(2) - 2) +
747  (GetBasisNumModes(1) - 2) * (GetBasisNumModes(2) - 2));
748 }
749 
750 int StdHexExp::v_GetTraceNumPoints(const int i) const
751 {
752  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
753 
754  if (i == 0 || i == 5)
755  {
756  return m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
757  }
758  else if (i == 1 || i == 3)
759  {
760  return m_base[0]->GetNumPoints() * m_base[2]->GetNumPoints();
761  }
762  else
763  {
764  return m_base[1]->GetNumPoints() * m_base[2]->GetNumPoints();
765  }
766 }
767 
768 LibUtilities::PointsKey StdHexExp::v_GetTracePointsKey(const int i,
769  const int j) const
770 {
771  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
772  ASSERTL2(j == 0 || j == 1, "face direction is out of range");
773 
774  if (i == 0 || i == 5)
775  {
776  return m_base[j]->GetPointsKey();
777  }
778  else if (i == 1 || i == 3)
779  {
780  return m_base[2 * j]->GetPointsKey();
781  }
782  else
783  {
784  return m_base[j + 1]->GetPointsKey();
785  }
786 }
787 
788 int StdHexExp::v_CalcNumberOfCoefficients(
789  const std::vector<unsigned int> &nummodes, int &modes_offset)
790 {
791  int nmodes = nummodes[modes_offset] * nummodes[modes_offset + 1] *
792  nummodes[modes_offset + 2];
793  modes_offset += 3;
794 
795  return nmodes;
796 }
797 
798 const LibUtilities::BasisKey StdHexExp::v_GetTraceBasisKey(const int i,
799  const int k) const
800 {
801  ASSERTL2(i >= 0 && i <= 5, "face id is out of range");
802  ASSERTL2(k >= 0 && k <= 1, "basis key id is out of range");
803 
804  int dir = k;
805  switch (i)
806  {
807  case 0:
808  case 5:
809  dir = k;
810  break;
811  case 1:
812  case 3:
813  dir = 2 * k;
814  break;
815  case 2:
816  case 4:
817  dir = k + 1;
818  break;
819  }
820 
821  return EvaluateQuadFaceBasisKey(k, m_base[dir]->GetBasisType(),
822  m_base[dir]->GetNumPoints(),
823  m_base[dir]->GetNumModes());
824 }
825 
826 void StdHexExp::v_GetCoords(Array<OneD, NekDouble> &xi_x,
827  Array<OneD, NekDouble> &xi_y,
828  Array<OneD, NekDouble> &xi_z)
829 {
830  Array<OneD, const NekDouble> eta_x = m_base[0]->GetZ();
831  Array<OneD, const NekDouble> eta_y = m_base[1]->GetZ();
832  Array<OneD, const NekDouble> eta_z = m_base[2]->GetZ();
833  int Qx = GetNumPoints(0);
834  int Qy = GetNumPoints(1);
835  int Qz = GetNumPoints(2);
836 
837  // Convert collapsed coordinates into cartesian coordinates:
838  // eta --> xi
839  for (int k = 0; k < Qz; ++k)
840  {
841  for (int j = 0; j < Qy; ++j)
842  {
843  for (int i = 0; i < Qx; ++i)
844  {
845  int s = i + Qx * (j + Qy * k);
846  xi_x[s] = eta_x[i];
847  xi_y[s] = eta_y[j];
848  xi_z[s] = eta_z[k];
849  }
850  }
851  }
852 }
853 
854 void StdHexExp::v_GetTraceNumModes(const int fid, int &numModes0,
855  int &numModes1, Orientation faceOrient)
856 {
857  int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
858  m_base[2]->GetNumModes()};
859  switch (fid)
860  {
861  case 0:
862  case 5:
863  {
864  numModes0 = nummodes[0];
865  numModes1 = nummodes[1];
866  }
867  break;
868  case 1:
869  case 3:
870  {
871  numModes0 = nummodes[0];
872  numModes1 = nummodes[2];
873  }
874  break;
875  case 2:
876  case 4:
877  {
878  numModes0 = nummodes[1];
879  numModes1 = nummodes[2];
880  }
881  break;
882  default:
883  {
884  ASSERTL0(false, "fid out of range");
885  }
886  break;
887  }
888 
889  if (faceOrient >= eDir1FwdDir2_Dir2FwdDir1)
890  {
891  std::swap(numModes0, numModes1);
892  }
893 }
894 
895 /**
896  * Expansions in each of the three dimensions must be of type
897  * LibUtilities#eModified_A or LibUtilities#eGLL_Lagrange.
898  *
899  * @param localVertexId ID of vertex (0..7)
900  * @returns Position of vertex in local numbering scheme.
901  */
902 int StdHexExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
903 {
904  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
905  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
906  "BasisType is not a boundary interior form");
907  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
908  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
909  "BasisType is not a boundary interior form");
910  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
911  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
912  "BasisType is not a boundary interior form");
913 
914  ASSERTL1((localVertexId >= 0) && (localVertexId < 8),
915  "local vertex id must be between 0 and 7");
916 
917  int p = 0;
918  int q = 0;
919  int r = 0;
920 
921  // Retrieve the number of modes in each dimension.
922  int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
923  m_base[2]->GetNumModes()};
924 
925  if (useCoeffPacking == true) // follow packing of coefficients i.e q,r,p
926  {
927  if (localVertexId > 3)
928  {
929  if (GetBasisType(2) == LibUtilities::eGLL_Lagrange)
930  {
931  r = nummodes[2] - 1;
932  }
933  else
934  {
935  r = 1;
936  }
937  }
938 
939  switch (localVertexId % 4)
940  {
941  case 0:
942  break;
943  case 1:
944  {
945  if (GetBasisType(0) == LibUtilities::eGLL_Lagrange)
946  {
947  p = nummodes[0] - 1;
948  }
949  else
950  {
951  p = 1;
952  }
953  }
954  break;
955  case 2:
956  {
957  if (GetBasisType(1) == LibUtilities::eGLL_Lagrange)
958  {
959  q = nummodes[1] - 1;
960  }
961  else
962  {
963  q = 1;
964  }
965  }
966  break;
967  case 3:
968  {
969  if (GetBasisType(1) == LibUtilities::eGLL_Lagrange)
970  {
971  p = nummodes[0] - 1;
972  q = nummodes[1] - 1;
973  }
974  else
975  {
976  p = 1;
977  q = 1;
978  }
979  }
980  break;
981  }
982  }
983  else
984  {
985  // Right face (vertices 1,2,5,6)
986  if ((localVertexId % 4) % 3 > 0)
987  {
988  if (GetBasisType(0) == LibUtilities::eGLL_Lagrange)
989  {
990  p = nummodes[0] - 1;
991  }
992  else
993  {
994  p = 1;
995  }
996  }
997  // Back face (vertices 2,3,6,7)
998  if (localVertexId % 4 > 1)
999  {
1000  if (GetBasisType(1) == LibUtilities::eGLL_Lagrange)
1001  {
1002  q = nummodes[1] - 1;
1003  }
1004  else
1005  {
1006  q = 1;
1007  }
1008  }
1009 
1010  // Top face (vertices 4,5,6,7)
1011  if (localVertexId > 3)
1012  {
1013  if (GetBasisType(2) == LibUtilities::eGLL_Lagrange)
1014  {
1015  r = nummodes[2] - 1;
1016  }
1017  else
1018  {
1019  r = 1;
1020  }
1021  }
1022  }
1023  // Compute the local number.
1024  return r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1025 }
1026 
1027 /**
1028  * @param outarray Storage area for computed map.
1029  */
1030 void StdHexExp::v_GetInteriorMap(Array<OneD, unsigned int> &outarray)
1031 {
1032  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1033  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1034  "BasisType is not a boundary interior form");
1035  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1036  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1037  "BasisType is not a boundary interior form");
1038  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1039  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1040  "BasisType is not a boundary interior form");
1041 
1042  int i;
1043  int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
1044  m_base[2]->GetNumModes()};
1045 
1046  int nIntCoeffs = m_ncoeffs - NumBndryCoeffs();
1047 
1048  if (outarray.size() != nIntCoeffs)
1049  {
1050  outarray = Array<OneD, unsigned int>(nIntCoeffs);
1051  }
1052 
1053  const LibUtilities::BasisType Btype[3] = {GetBasisType(0), GetBasisType(1),
1054  GetBasisType(2)};
1055 
1056  int p, q, r;
1057  int cnt = 0;
1058 
1059  int IntIdx[3][2];
1060 
1061  for (i = 0; i < 3; i++)
1062  {
1063  if (Btype[i] == LibUtilities::eModified_A)
1064  {
1065  IntIdx[i][0] = 2;
1066  IntIdx[i][1] = nummodes[i];
1067  }
1068  else
1069  {
1070  IntIdx[i][0] = 1;
1071  IntIdx[i][1] = nummodes[i] - 1;
1072  }
1073  }
1074 
1075  for (r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
1076  {
1077  for (q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
1078  {
1079  for (p = IntIdx[0][0]; p < IntIdx[0][1]; p++)
1080  {
1081  outarray[cnt++] =
1082  r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1083  }
1084  }
1085  }
1086 }
1087 
1088 /**
1089  * @param outarray Storage for computed map.
1090  */
1091 void StdHexExp::v_GetBoundaryMap(Array<OneD, unsigned int> &outarray)
1092 {
1093  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1094  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1095  "BasisType is not a boundary interior form");
1096  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1097  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1098  "BasisType is not a boundary interior form");
1099  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1100  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1101  "BasisType is not a boundary interior form");
1102 
1103  int i;
1104  int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
1105  m_base[2]->GetNumModes()};
1106 
1107  int nBndCoeffs = NumBndryCoeffs();
1108 
1109  if (outarray.size() != nBndCoeffs)
1110  {
1111  outarray = Array<OneD, unsigned int>(nBndCoeffs);
1112  }
1113 
1114  const LibUtilities::BasisType Btype[3] = {GetBasisType(0), GetBasisType(1),
1115  GetBasisType(2)};
1116 
1117  int p, q, r;
1118  int cnt = 0;
1119 
1120  int BndIdx[3][2];
1121  int IntIdx[3][2];
1122 
1123  for (i = 0; i < 3; i++)
1124  {
1125  BndIdx[i][0] = 0;
1126 
1127  if (Btype[i] == LibUtilities::eModified_A)
1128  {
1129  BndIdx[i][1] = 1;
1130  IntIdx[i][0] = 2;
1131  IntIdx[i][1] = nummodes[i];
1132  }
1133  else
1134  {
1135  BndIdx[i][1] = nummodes[i] - 1;
1136  IntIdx[i][0] = 1;
1137  IntIdx[i][1] = nummodes[i] - 1;
1138  }
1139  }
1140 
1141  for (i = 0; i < 2; i++)
1142  {
1143  r = BndIdx[2][i];
1144  for (q = 0; q < nummodes[1]; q++)
1145  {
1146  for (p = 0; p < nummodes[0]; p++)
1147  {
1148  outarray[cnt++] =
1149  r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1150  }
1151  }
1152  }
1153 
1154  for (r = IntIdx[2][0]; r < IntIdx[2][1]; r++)
1155  {
1156  for (i = 0; i < 2; i++)
1157  {
1158  q = BndIdx[1][i];
1159  for (p = 0; p < nummodes[0]; p++)
1160  {
1161  outarray[cnt++] =
1162  r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1163  }
1164  }
1165 
1166  for (q = IntIdx[1][0]; q < IntIdx[1][1]; q++)
1167  {
1168  for (i = 0; i < 2; i++)
1169  {
1170  p = BndIdx[0][i];
1171  outarray[cnt++] =
1172  r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1173  }
1174  }
1175  }
1176 
1177  sort(outarray.get(), outarray.get() + nBndCoeffs);
1178 }
1179 
1180 /**
1181  * Only for basis type Modified_A or GLL_LAGRANGE in all directions.
1182  */
1183 void StdHexExp::v_GetTraceCoeffMap(const unsigned int fid,
1184  Array<OneD, unsigned int> &maparray)
1185 {
1186  int i, j;
1187  int nummodesA = 0, nummodesB = 0;
1188 
1189  ASSERTL1(GetBasisType(0) == GetBasisType(1) &&
1190  GetBasisType(0) == GetBasisType(2),
1191  "Method only implemented if BasisType is indentical in "
1192  "all directions");
1193  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1194  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1195  "Method only implemented for Modified_A or "
1196  "GLL_Lagrange BasisType");
1197 
1198  const int nummodes0 = m_base[0]->GetNumModes();
1199  const int nummodes1 = m_base[1]->GetNumModes();
1200  const int nummodes2 = m_base[2]->GetNumModes();
1201 
1202  switch (fid)
1203  {
1204  case 0:
1205  case 5:
1206  nummodesA = nummodes0;
1207  nummodesB = nummodes1;
1208  break;
1209  case 1:
1210  case 3:
1211  nummodesA = nummodes0;
1212  nummodesB = nummodes2;
1213  break;
1214  case 2:
1215  case 4:
1216  nummodesA = nummodes1;
1217  nummodesB = nummodes2;
1218  break;
1219  default:
1220  ASSERTL0(false, "fid must be between 0 and 5");
1221  }
1222 
1223  int nFaceCoeffs = nummodesA * nummodesB;
1224 
1225  if (maparray.size() != nFaceCoeffs)
1226  {
1227  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1228  }
1229 
1230  bool modified = (GetBasisType(0) == LibUtilities::eModified_A);
1231 
1232  int offset = 0;
1233  int jump1 = 1;
1234  int jump2 = 1;
1235 
1236  switch (fid)
1237  {
1238  case 5:
1239  {
1240  if (modified)
1241  {
1242  offset = nummodes0 * nummodes1;
1243  }
1244  else
1245  {
1246  offset = (nummodes2 - 1) * nummodes0 * nummodes1;
1247  jump1 = nummodes0;
1248  }
1249  }
1250  /* Falls through. */
1251  case 0:
1252  {
1253  jump1 = nummodes0;
1254  break;
1255  }
1256  case 3:
1257  {
1258  if (modified)
1259  {
1260  offset = nummodes0;
1261  }
1262  else
1263  {
1264  offset = nummodes0 * (nummodes1 - 1);
1265  jump1 = nummodes0 * nummodes1;
1266  }
1267  }
1268  /* Falls through. */
1269  case 1:
1270  {
1271  jump1 = nummodes0 * nummodes1;
1272  break;
1273  }
1274  case 2:
1275  {
1276  if (modified)
1277  {
1278  offset = 1;
1279  }
1280  else
1281  {
1282  offset = nummodes0 - 1;
1283  jump1 = nummodes0 * nummodes1;
1284  jump2 = nummodes0;
1285  }
1286  }
1287  /* Falls through. */
1288  case 4:
1289  {
1290  jump1 = nummodes0 * nummodes1;
1291  jump2 = nummodes0;
1292  break;
1293  }
1294  default:
1295  ASSERTL0(false, "fid must be between 0 and 5");
1296  }
1297 
1298  for (i = 0; i < nummodesB; i++)
1299  {
1300  for (j = 0; j < nummodesA; j++)
1301  {
1302  maparray[i * nummodesA + j] = i * jump1 + j * jump2 + offset;
1303  }
1304  }
1305 }
1306 
1307 void StdHexExp::v_GetElmtTraceToTraceMap(const unsigned int fid,
1308  Array<OneD, unsigned int> &maparray,
1309  Array<OneD, int> &signarray,
1310  Orientation faceOrient, int P, int Q)
1311 {
1312  int i, j;
1313  int nummodesA = 0, nummodesB = 0;
1314 
1315  ASSERTL1(GetBasisType(0) == GetBasisType(1) &&
1316  GetBasisType(0) == GetBasisType(2),
1317  "Method only implemented if BasisType is indentical in "
1318  "all directions");
1319  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1320  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1321  "Method only implemented for Modified_A or "
1322  "GLL_Lagrange BasisType");
1323 
1324  const int nummodes0 = m_base[0]->GetNumModes();
1325  const int nummodes1 = m_base[1]->GetNumModes();
1326  const int nummodes2 = m_base[2]->GetNumModes();
1327 
1328  switch (fid)
1329  {
1330  case 0:
1331  case 5:
1332  nummodesA = nummodes0;
1333  nummodesB = nummodes1;
1334  break;
1335  case 1:
1336  case 3:
1337  nummodesA = nummodes0;
1338  nummodesB = nummodes2;
1339  break;
1340  case 2:
1341  case 4:
1342  nummodesA = nummodes1;
1343  nummodesB = nummodes2;
1344  break;
1345  default:
1346  ASSERTL0(false, "fid must be between 0 and 5");
1347  }
1348 
1349  if (P == -1)
1350  {
1351  P = nummodesA;
1352  Q = nummodesB;
1353  }
1354 
1355  bool modified = (GetBasisType(0) == LibUtilities::eModified_A);
1356 
1357  // check that
1358  if (modified == false)
1359  {
1360  ASSERTL1((P == nummodesA) || (Q == nummodesB),
1361  "Different trace space face dimention "
1362  "and element face dimention not possible for "
1363  "GLL-Lagrange bases");
1364  }
1365 
1366  int nFaceCoeffs = P * Q;
1367 
1368  if (maparray.size() != nFaceCoeffs)
1369  {
1370  maparray = Array<OneD, unsigned int>(nFaceCoeffs);
1371  }
1372 
1373  // fill default mapping as increasing index
1374  for (i = 0; i < nFaceCoeffs; ++i)
1375  {
1376  maparray[i] = i;
1377  }
1378 
1379  if (signarray.size() != nFaceCoeffs)
1380  {
1381  signarray = Array<OneD, int>(nFaceCoeffs, 1);
1382  }
1383  else
1384  {
1385  fill(signarray.get(), signarray.get() + nFaceCoeffs, 1);
1386  }
1387 
1388  // setup indexing to manage transpose directions
1389  Array<OneD, int> arrayindx(nFaceCoeffs);
1390  for (i = 0; i < Q; i++)
1391  {
1392  for (j = 0; j < P; j++)
1393  {
1394  if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1395  {
1396  arrayindx[i * P + j] = i * P + j;
1397  }
1398  else
1399  {
1400  arrayindx[i * P + j] = j * Q + i;
1401  }
1402  }
1403  }
1404 
1405  // zero signmap and set maparray to zero if elemental
1406  // modes are not as large as face models
1407  for (i = 0; i < nummodesB; i++)
1408  {
1409  for (j = nummodesA; j < P; j++)
1410  {
1411  signarray[arrayindx[i * P + j]] = 0.0;
1412  maparray[arrayindx[i * P + j]] = maparray[0];
1413  }
1414  }
1415 
1416  for (i = nummodesB; i < Q; i++)
1417  {
1418  for (j = 0; j < P; j++)
1419  {
1420  signarray[arrayindx[i * P + j]] = 0.0;
1421  maparray[arrayindx[i * P + j]] = maparray[0];
1422  }
1423  }
1424 
1425  // zero signmap and set maparray to zero entry if
1426  // elemental modes are not as large as face modesl
1427  for (i = 0; i < Q; i++)
1428  {
1429  // fill values into map array of trace size
1430  // for element face index
1431  for (j = 0; j < P; j++)
1432  {
1433  maparray[arrayindx[i * P + j]] = i * nummodesA + j;
1434  }
1435 
1436  // zero values if P > numModesA
1437  for (j = nummodesA; j < P; j++)
1438  {
1439  signarray[arrayindx[i * P + j]] = 0.0;
1440  maparray[arrayindx[i * P + j]] = maparray[0];
1441  }
1442  }
1443 
1444  // zero values if Q > numModesB
1445  for (i = nummodesB; i < Q; i++)
1446  {
1447  for (j = 0; j < P; j++)
1448  {
1449  signarray[arrayindx[i * P + j]] = 0.0;
1450  maparray[arrayindx[i * P + j]] = maparray[0];
1451  }
1452  }
1453 
1454  // Now reorientate indices accordign to orientation
1455  if ((faceOrient == eDir1FwdDir1_Dir2BwdDir2) ||
1456  (faceOrient == eDir1BwdDir1_Dir2BwdDir2) ||
1457  (faceOrient == eDir1BwdDir2_Dir2FwdDir1) ||
1458  (faceOrient == eDir1BwdDir2_Dir2BwdDir1))
1459  {
1460  if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1461  {
1462  if (modified)
1463  {
1464  for (int i = 3; i < Q; i += 2)
1465  {
1466  for (int j = 0; j < P; j++)
1467  {
1468  signarray[arrayindx[i * P + j]] *= -1;
1469  }
1470  }
1471 
1472  for (int i = 0; i < P; i++)
1473  {
1474  swap(maparray[i], maparray[i + P]);
1475  swap(signarray[i], signarray[i + P]);
1476  }
1477  }
1478  else
1479  {
1480  for (int i = 0; i < P; i++)
1481  {
1482  for (int j = 0; j < Q / 2; j++)
1483  {
1484  swap(maparray[i + j * P],
1485  maparray[i + P * Q - P - j * P]);
1486  swap(signarray[i + j * P],
1487  signarray[i + P * Q - P - j * P]);
1488  }
1489  }
1490  }
1491  }
1492  else
1493  {
1494  if (modified)
1495  {
1496  for (int i = 0; i < Q; i++)
1497  {
1498  for (int j = 3; j < P; j += 2)
1499  {
1500  signarray[arrayindx[i * P + j]] *= -1;
1501  }
1502  }
1503 
1504  for (int i = 0; i < Q; i++)
1505  {
1506  swap(maparray[i], maparray[i + Q]);
1507  swap(signarray[i], signarray[i + Q]);
1508  }
1509  }
1510  else
1511  {
1512  for (int i = 0; i < P; i++)
1513  {
1514  for (int j = 0; j < Q / 2; j++)
1515  {
1516  swap(maparray[i * Q + j], maparray[i * Q + Q - 1 - j]);
1517  swap(signarray[i * Q + j],
1518  signarray[i * Q + Q - 1 - j]);
1519  }
1520  }
1521  }
1522  }
1523  }
1524 
1525  if ((faceOrient == eDir1BwdDir1_Dir2FwdDir2) ||
1526  (faceOrient == eDir1BwdDir1_Dir2BwdDir2) ||
1527  (faceOrient == eDir1FwdDir2_Dir2BwdDir1) ||
1528  (faceOrient == eDir1BwdDir2_Dir2BwdDir1))
1529  {
1530  if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1531  {
1532  if (modified)
1533  {
1534  for (i = 0; i < Q; i++)
1535  {
1536  for (j = 3; j < P; j += 2)
1537  {
1538  signarray[arrayindx[i * P + j]] *= -1;
1539  }
1540  }
1541 
1542  for (i = 0; i < Q; i++)
1543  {
1544  swap(maparray[i * P], maparray[i * P + 1]);
1545  swap(signarray[i * P], signarray[i * P + 1]);
1546  }
1547  }
1548  else
1549  {
1550  for (i = 0; i < Q; i++)
1551  {
1552  for (j = 0; j < P / 2; j++)
1553  {
1554  swap(maparray[i * P + j], maparray[i * P + P - 1 - j]);
1555  swap(signarray[i * P + j],
1556  signarray[i * P + P - 1 - j]);
1557  }
1558  }
1559  }
1560  }
1561  else
1562  {
1563  if (modified)
1564  {
1565  for (i = 3; i < Q; i += 2)
1566  {
1567  for (j = 0; j < P; j++)
1568  {
1569  signarray[arrayindx[i * P + j]] *= -1;
1570  }
1571  }
1572 
1573  for (i = 0; i < P; i++)
1574  {
1575  swap(maparray[i * Q], maparray[i * Q + 1]);
1576  swap(signarray[i * Q], signarray[i * Q + 1]);
1577  }
1578  }
1579  else
1580  {
1581  for (i = 0; i < Q; i++)
1582  {
1583  for (j = 0; j < P / 2; j++)
1584  {
1585  swap(maparray[i + j * Q],
1586  maparray[i + P * Q - Q - j * Q]);
1587  swap(signarray[i + j * Q],
1588  signarray[i + P * Q - Q - j * Q]);
1589  }
1590  }
1591  }
1592  }
1593  }
1594 }
1595 
1596 /**
1597  * @param eid The edge to compute the numbering for.
1598  * @param edgeOrient Orientation of the edge.
1599  * @param maparray Storage for computed mapping array.
1600  * @param signarray ?
1601  */
1602 void StdHexExp::v_GetEdgeInteriorToElementMap(
1603  const int eid, Array<OneD, unsigned int> &maparray,
1604  Array<OneD, int> &signarray, const Orientation edgeOrient)
1605 {
1606  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1607  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1608  "BasisType is not a boundary interior form");
1609  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1610  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1611  "BasisType is not a boundary interior form");
1612  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1613  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1614  "BasisType is not a boundary interior form");
1615 
1616  ASSERTL1((eid >= 0) && (eid < 12),
1617  "local edge id must be between 0 and 11");
1618 
1619  int nEdgeIntCoeffs = GetEdgeNcoeffs(eid) - 2;
1620 
1621  if (maparray.size() != nEdgeIntCoeffs)
1622  {
1623  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1624  }
1625 
1626  if (signarray.size() != nEdgeIntCoeffs)
1627  {
1628  signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1629  }
1630  else
1631  {
1632  fill(signarray.get(), signarray.get() + nEdgeIntCoeffs, 1);
1633  }
1634 
1635  int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
1636  m_base[2]->GetNumModes()};
1637 
1638  const LibUtilities::BasisType bType[3] = {GetBasisType(0), GetBasisType(1),
1639  GetBasisType(2)};
1640 
1641  bool reverseOrdering = false;
1642  bool signChange = false;
1643 
1644  int IdxRange[3][2] = {{0, 0}, {0, 0}, {0, 0}};
1645 
1646  switch (eid)
1647  {
1648  case 0:
1649  case 1:
1650  case 2:
1651  case 3:
1652  {
1653  IdxRange[2][0] = 0;
1654  IdxRange[2][1] = 1;
1655  }
1656  break;
1657  case 8:
1658  case 9:
1659  case 10:
1660  case 11:
1661  {
1662  if (bType[2] == LibUtilities::eGLL_Lagrange)
1663  {
1664  IdxRange[2][0] = nummodes[2] - 1;
1665  IdxRange[2][1] = nummodes[2];
1666  }
1667  else
1668  {
1669  IdxRange[2][0] = 1;
1670  IdxRange[2][1] = 2;
1671  }
1672  }
1673  break;
1674  case 4:
1675  case 5:
1676  case 6:
1677  case 7:
1678  {
1679  if (bType[2] == LibUtilities::eGLL_Lagrange)
1680  {
1681  IdxRange[2][0] = 1;
1682  IdxRange[2][1] = nummodes[2] - 1;
1683 
1684  if (edgeOrient == eBackwards)
1685  {
1686  reverseOrdering = true;
1687  }
1688  }
1689  else
1690  {
1691  IdxRange[2][0] = 2;
1692  IdxRange[2][1] = nummodes[2];
1693 
1694  if (edgeOrient == eBackwards)
1695  {
1696  signChange = true;
1697  }
1698  }
1699  }
1700  break;
1701  }
1702 
1703  switch (eid)
1704  {
1705  case 0:
1706  case 4:
1707  case 5:
1708  case 8:
1709  {
1710  IdxRange[1][0] = 0;
1711  IdxRange[1][1] = 1;
1712  }
1713  break;
1714  case 2:
1715  case 6:
1716  case 7:
1717  case 10:
1718  {
1719  if (bType[1] == LibUtilities::eGLL_Lagrange)
1720  {
1721  IdxRange[1][0] = nummodes[1] - 1;
1722  IdxRange[1][1] = nummodes[1];
1723  }
1724  else
1725  {
1726  IdxRange[1][0] = 1;
1727  IdxRange[1][1] = 2;
1728  }
1729  }
1730  break;
1731  case 1:
1732  case 9:
1733  {
1734  if (bType[1] == LibUtilities::eGLL_Lagrange)
1735  {
1736  IdxRange[1][0] = 1;
1737  IdxRange[1][1] = nummodes[1] - 1;
1738 
1739  if (edgeOrient == eBackwards)
1740  {
1741  reverseOrdering = true;
1742  }
1743  }
1744  else
1745  {
1746  IdxRange[1][0] = 2;
1747  IdxRange[1][1] = nummodes[1];
1748 
1749  if (edgeOrient == eBackwards)
1750  {
1751  signChange = true;
1752  }
1753  }
1754  }
1755  break;
1756  case 3:
1757  case 11:
1758  {
1759  if (bType[1] == LibUtilities::eGLL_Lagrange)
1760  {
1761  IdxRange[1][0] = 1;
1762  IdxRange[1][1] = nummodes[1] - 1;
1763 
1764  if (edgeOrient == eForwards)
1765  {
1766  reverseOrdering = true;
1767  }
1768  }
1769  else
1770  {
1771  IdxRange[1][0] = 2;
1772  IdxRange[1][1] = nummodes[1];
1773 
1774  if (edgeOrient == eForwards)
1775  {
1776  signChange = true;
1777  }
1778  }
1779  }
1780  break;
1781  }
1782 
1783  switch (eid)
1784  {
1785  case 3:
1786  case 4:
1787  case 7:
1788  case 11:
1789  {
1790  IdxRange[0][0] = 0;
1791  IdxRange[0][1] = 1;
1792  }
1793  break;
1794  case 1:
1795  case 5:
1796  case 6:
1797  case 9:
1798  {
1799  if (bType[0] == LibUtilities::eGLL_Lagrange)
1800  {
1801  IdxRange[0][0] = nummodes[0] - 1;
1802  IdxRange[0][1] = nummodes[0];
1803  }
1804  else
1805  {
1806  IdxRange[0][0] = 1;
1807  IdxRange[0][1] = 2;
1808  }
1809  }
1810  break;
1811  case 0:
1812  case 8:
1813  {
1814  if (bType[0] == LibUtilities::eGLL_Lagrange)
1815  {
1816  IdxRange[0][0] = 1;
1817  IdxRange[0][1] = nummodes[0] - 1;
1818 
1819  if (edgeOrient == eBackwards)
1820  {
1821  reverseOrdering = true;
1822  }
1823  }
1824  else
1825  {
1826  IdxRange[0][0] = 2;
1827  IdxRange[0][1] = nummodes[0];
1828 
1829  if (edgeOrient == eBackwards)
1830  {
1831  signChange = true;
1832  }
1833  }
1834  }
1835  break;
1836  case 2:
1837  case 10:
1838  {
1839  if (bType[0] == LibUtilities::eGLL_Lagrange)
1840  {
1841  IdxRange[0][0] = 1;
1842  IdxRange[0][1] = nummodes[0] - 1;
1843 
1844  if (edgeOrient == eForwards)
1845  {
1846  reverseOrdering = true;
1847  }
1848  }
1849  else
1850  {
1851  IdxRange[0][0] = 2;
1852  IdxRange[0][1] = nummodes[0];
1853 
1854  if (edgeOrient == eForwards)
1855  {
1856  signChange = true;
1857  }
1858  }
1859  }
1860  break;
1861  }
1862 
1863  int cnt = 0;
1864 
1865  for (int r = IdxRange[2][0]; r < IdxRange[2][1]; r++)
1866  {
1867  for (int q = IdxRange[1][0]; q < IdxRange[1][1]; q++)
1868  {
1869  for (int p = IdxRange[0][0]; p < IdxRange[0][1]; p++)
1870  {
1871  maparray[cnt++] =
1872  r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
1873  }
1874  }
1875  }
1876 
1877  if (reverseOrdering)
1878  {
1879  reverse(maparray.get(), maparray.get() + nEdgeIntCoeffs);
1880  }
1881 
1882  if (signChange)
1883  {
1884  for (int p = 1; p < nEdgeIntCoeffs; p += 2)
1885  {
1886  signarray[p] = -1;
1887  }
1888  }
1889 }
1890 
1891 /**
1892  * Generate mapping describing which elemental modes lie on the
1893  * interior of a given face. Accounts for face orientation.
1894  */
1895 void StdHexExp::v_GetTraceInteriorToElementMap(
1896  const int fid, Array<OneD, unsigned int> &maparray,
1897  Array<OneD, int> &signarray, const Orientation faceOrient)
1898 {
1899  ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A ||
1900  GetBasisType(0) == LibUtilities::eGLL_Lagrange,
1901  "BasisType is not a boundary interior form");
1902  ASSERTL1(GetBasisType(1) == LibUtilities::eModified_A ||
1903  GetBasisType(1) == LibUtilities::eGLL_Lagrange,
1904  "BasisType is not a boundary interior form");
1905  ASSERTL1(GetBasisType(2) == LibUtilities::eModified_A ||
1906  GetBasisType(2) == LibUtilities::eGLL_Lagrange,
1907  "BasisType is not a boundary interior form");
1908 
1909  ASSERTL1((fid >= 0) && (fid < 6), "local face id must be between 0 and 5");
1910 
1911  int nFaceIntCoeffs = v_GetTraceIntNcoeffs(fid);
1912 
1913  if (maparray.size() != nFaceIntCoeffs)
1914  {
1915  maparray = Array<OneD, unsigned int>(nFaceIntCoeffs);
1916  }
1917 
1918  if (signarray.size() != nFaceIntCoeffs)
1919  {
1920  signarray = Array<OneD, int>(nFaceIntCoeffs, 1);
1921  }
1922  else
1923  {
1924  fill(signarray.get(), signarray.get() + nFaceIntCoeffs, 1);
1925  }
1926 
1927  int nummodes[3] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes(),
1928  m_base[2]->GetNumModes()};
1929 
1930  const LibUtilities::BasisType bType[3] = {GetBasisType(0), GetBasisType(1),
1931  GetBasisType(2)};
1932 
1933  int nummodesA = 0;
1934  int nummodesB = 0;
1935 
1936  // Determine the number of modes in face directions A & B based
1937  // on the face index given.
1938  switch (fid)
1939  {
1940  case 0:
1941  case 5:
1942  {
1943  nummodesA = nummodes[0];
1944  nummodesB = nummodes[1];
1945  }
1946  break;
1947  case 1:
1948  case 3:
1949  {
1950  nummodesA = nummodes[0];
1951  nummodesB = nummodes[2];
1952  }
1953  break;
1954  case 2:
1955  case 4:
1956  {
1957  nummodesA = nummodes[1];
1958  nummodesB = nummodes[2];
1959  }
1960  }
1961 
1962  Array<OneD, int> arrayindx(nFaceIntCoeffs);
1963 
1964  // Create a mapping array to account for transposition of the
1965  // coordinates due to face orientation.
1966  for (int i = 0; i < (nummodesB - 2); i++)
1967  {
1968  for (int j = 0; j < (nummodesA - 2); j++)
1969  {
1970  if (faceOrient < eDir1FwdDir2_Dir2FwdDir1)
1971  {
1972  arrayindx[i * (nummodesA - 2) + j] = i * (nummodesA - 2) + j;
1973  }
1974  else
1975  {
1976  arrayindx[i * (nummodesA - 2) + j] = j * (nummodesB - 2) + i;
1977  }
1978  }
1979  }
1980 
1981  int IdxRange[3][2];
1982  int Incr[3];
1983 
1984  Array<OneD, int> sign0(nummodes[0], 1);
1985  Array<OneD, int> sign1(nummodes[1], 1);
1986  Array<OneD, int> sign2(nummodes[2], 1);
1987 
1988  // Set the upper and lower bounds, and increment for the faces
1989  // involving the first coordinate direction.
1990  switch (fid)
1991  {
1992  case 0: // bottom face
1993  {
1994  IdxRange[2][0] = 0;
1995  IdxRange[2][1] = 1;
1996  Incr[2] = 1;
1997  }
1998  break;
1999  case 5: // top face
2000  {
2001  if (bType[2] == LibUtilities::eGLL_Lagrange)
2002  {
2003  IdxRange[2][0] = nummodes[2] - 1;
2004  IdxRange[2][1] = nummodes[2];
2005  Incr[2] = 1;
2006  }
2007  else
2008  {
2009  IdxRange[2][0] = 1;
2010  IdxRange[2][1] = 2;
2011  Incr[2] = 1;
2012  }
2013  }
2014  break;
2015  default: // all other faces
2016  {
2017  if (bType[2] == LibUtilities::eGLL_Lagrange)
2018  {
2019  if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 2)
2020  {
2021  IdxRange[2][0] = nummodes[2] - 2;
2022  IdxRange[2][1] = 0;
2023  Incr[2] = -1;
2024  }
2025  else
2026  {
2027  IdxRange[2][0] = 1;
2028  IdxRange[2][1] = nummodes[2] - 1;
2029  Incr[2] = 1;
2030  }
2031  }
2032  else
2033  {
2034  IdxRange[2][0] = 2;
2035  IdxRange[2][1] = nummodes[2];
2036  Incr[2] = 1;
2037 
2038  if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 2)
2039  {
2040  for (int i = 3; i < nummodes[2]; i += 2)
2041  {
2042  sign2[i] = -1;
2043  }
2044  }
2045  }
2046  }
2047  }
2048 
2049  // Set the upper and lower bounds, and increment for the faces
2050  // involving the second coordinate direction.
2051  switch (fid)
2052  {
2053  case 1:
2054  {
2055  IdxRange[1][0] = 0;
2056  IdxRange[1][1] = 1;
2057  Incr[1] = 1;
2058  }
2059  break;
2060  case 3:
2061  {
2062  if (bType[1] == LibUtilities::eGLL_Lagrange)
2063  {
2064  IdxRange[1][0] = nummodes[1] - 1;
2065  IdxRange[1][1] = nummodes[1];
2066  Incr[1] = 1;
2067  }
2068  else
2069  {
2070  IdxRange[1][0] = 1;
2071  IdxRange[1][1] = 2;
2072  Incr[1] = 1;
2073  }
2074  }
2075  break;
2076  case 0:
2077  case 5:
2078  {
2079  if (bType[1] == LibUtilities::eGLL_Lagrange)
2080  {
2081  if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 2)
2082  {
2083  IdxRange[1][0] = nummodes[1] - 2;
2084  IdxRange[1][1] = 0;
2085  Incr[1] = -1;
2086  }
2087  else
2088  {
2089  IdxRange[1][0] = 1;
2090  IdxRange[1][1] = nummodes[1] - 1;
2091  Incr[1] = 1;
2092  }
2093  }
2094  else
2095  {
2096  IdxRange[1][0] = 2;
2097  IdxRange[1][1] = nummodes[1];
2098  Incr[1] = 1;
2099 
2100  if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 2)
2101  {
2102  for (int i = 3; i < nummodes[1]; i += 2)
2103  {
2104  sign1[i] = -1;
2105  }
2106  }
2107  }
2108  }
2109  break;
2110  default: // case2: case4:
2111  {
2112  if (bType[1] == LibUtilities::eGLL_Lagrange)
2113  {
2114  if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1)
2115  {
2116  IdxRange[1][0] = nummodes[1] - 2;
2117  IdxRange[1][1] = 0;
2118  Incr[1] = -1;
2119  }
2120  else
2121  {
2122  IdxRange[1][0] = 1;
2123  IdxRange[1][1] = nummodes[1] - 1;
2124  Incr[1] = 1;
2125  }
2126  }
2127  else
2128  {
2129  IdxRange[1][0] = 2;
2130  IdxRange[1][1] = nummodes[1];
2131  Incr[1] = 1;
2132 
2133  if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1)
2134  {
2135  for (int i = 3; i < nummodes[1]; i += 2)
2136  {
2137  sign1[i] = -1;
2138  }
2139  }
2140  }
2141  }
2142  }
2143 
2144  switch (fid)
2145  {
2146  case 4:
2147  {
2148  IdxRange[0][0] = 0;
2149  IdxRange[0][1] = 1;
2150  Incr[0] = 1;
2151  }
2152  break;
2153  case 2:
2154  {
2155  if (bType[0] == LibUtilities::eGLL_Lagrange)
2156  {
2157  IdxRange[0][0] = nummodes[0] - 1;
2158  IdxRange[0][1] = nummodes[0];
2159  Incr[0] = 1;
2160  }
2161  else
2162  {
2163  IdxRange[0][0] = 1;
2164  IdxRange[0][1] = 2;
2165  Incr[0] = 1;
2166  }
2167  }
2168  break;
2169  default:
2170  {
2171  if (bType[0] == LibUtilities::eGLL_Lagrange)
2172  {
2173  if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1)
2174  {
2175  IdxRange[0][0] = nummodes[0] - 2;
2176  IdxRange[0][1] = 0;
2177  Incr[0] = -1;
2178  }
2179  else
2180  {
2181  IdxRange[0][0] = 1;
2182  IdxRange[0][1] = nummodes[0] - 1;
2183  Incr[0] = 1;
2184  }
2185  }
2186  else
2187  {
2188  IdxRange[0][0] = 2;
2189  IdxRange[0][1] = nummodes[0];
2190  Incr[0] = 1;
2191 
2192  if (((int)(faceOrient - eDir1FwdDir1_Dir2FwdDir2)) % 4 > 1)
2193  {
2194  for (int i = 3; i < nummodes[0]; i += 2)
2195  {
2196  sign0[i] = -1;
2197  }
2198  }
2199  }
2200  }
2201  }
2202 
2203  int cnt = 0;
2204 
2205  for (int r = IdxRange[2][0]; r != IdxRange[2][1]; r += Incr[2])
2206  {
2207  for (int q = IdxRange[1][0]; q != IdxRange[1][1]; q += Incr[1])
2208  {
2209  for (int p = IdxRange[0][0]; p != IdxRange[0][1]; p += Incr[0])
2210  {
2211  maparray[arrayindx[cnt]] =
2212  r * nummodes[0] * nummodes[1] + q * nummodes[0] + p;
2213  signarray[arrayindx[cnt++]] = sign0[p] * sign1[q] * sign2[r];
2214  }
2215  }
2216  }
2217 }
2218 
2219 int StdHexExp::v_GetEdgeNcoeffs(const int i) const
2220 {
2221  ASSERTL2((i >= 0) && (i <= 11), "edge id is out of range");
2222 
2223  if ((i == 0) || (i == 2) || (i == 8) || (i == 10))
2224  {
2225  return GetBasisNumModes(0);
2226  }
2227  else if ((i == 1) || (i == 3) || (i == 9) || (i == 11))
2228  {
2229  return GetBasisNumModes(1);
2230  }
2231  else
2232  {
2233  return GetBasisNumModes(2);
2234  }
2235 }
2236 
2237 DNekMatSharedPtr StdHexExp::v_GenMatrix(const StdMatrixKey &mkey)
2238 {
2239 
2240  MatrixType mtype = mkey.GetMatrixType();
2241 
2242  DNekMatSharedPtr Mat;
2243 
2244  switch (mtype)
2245  {
2246  case ePhysInterpToEquiSpaced:
2247  {
2248  int nq0 = m_base[0]->GetNumPoints();
2249  int nq1 = m_base[1]->GetNumPoints();
2250  int nq2 = m_base[2]->GetNumPoints();
2251  int nq;
2252 
2253  // take definition from key
2254  if (mkey.ConstFactorExists(eFactorConst))
2255  {
2256  nq = (int)mkey.GetConstFactor(eFactorConst);
2257  }
2258  else
2259  {
2260  nq = max(nq0, max(nq1, nq2));
2261  }
2262 
2263  int neq =
2264  LibUtilities::StdHexData::getNumberOfCoefficients(nq, nq, nq);
2265  Array<OneD, Array<OneD, NekDouble>> coords(neq);
2266  Array<OneD, NekDouble> coll(3);
2267  Array<OneD, DNekMatSharedPtr> I(3);
2268  Array<OneD, NekDouble> tmp(nq0);
2269 
2270  Mat =
2271  MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1 * nq2);
2272  int cnt = 0;
2273 
2274  for (int i = 0; i < nq; ++i)
2275  {
2276  for (int j = 0; j < nq; ++j)
2277  {
2278  for (int k = 0; k < nq; ++k, ++cnt)
2279  {
2280  coords[cnt] = Array<OneD, NekDouble>(3);
2281  coords[cnt][0] = -1.0 + 2 * k / (NekDouble)(nq - 1);
2282  coords[cnt][1] = -1.0 + 2 * j / (NekDouble)(nq - 1);
2283  coords[cnt][2] = -1.0 + 2 * i / (NekDouble)(nq - 1);
2284  }
2285  }
2286  }
2287 
2288  for (int i = 0; i < neq; ++i)
2289  {
2290  LocCoordToLocCollapsed(coords[i], coll);
2291 
2292  I[0] = m_base[0]->GetI(coll);
2293  I[1] = m_base[1]->GetI(coll + 1);
2294  I[2] = m_base[2]->GetI(coll + 2);
2295 
2296  // interpolate first coordinate direction
2297  NekDouble fac;
2298  for (int k = 0; k < nq2; ++k)
2299  {
2300  for (int j = 0; j < nq1; ++j)
2301  {
2302 
2303  fac = (I[1]->GetPtr())[j] * (I[2]->GetPtr())[k];
2304  Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
2305 
2306  Vmath::Vcopy(nq0, &tmp[0], 1,
2307  Mat->GetRawPtr() + k * nq0 * nq1 * neq +
2308  j * nq0 * neq + i,
2309  neq);
2310  }
2311  }
2312  }
2313  }
2314  break;
2315  default:
2316  {
2317  Mat = StdExpansion::CreateGeneralMatrix(mkey);
2318  }
2319  break;
2320  }
2321 
2322  return Mat;
2323 }
2324 
2325 DNekMatSharedPtr StdHexExp::v_CreateStdMatrix(const StdMatrixKey &mkey)
2326 {
2327  return v_GenMatrix(mkey);
2328 }
2329 
2330 void StdHexExp::v_MassMatrixOp(const Array<OneD, const NekDouble> &inarray,
2331  Array<OneD, NekDouble> &outarray,
2332  const StdMatrixKey &mkey)
2333 {
2334  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
2335 }
2336 
2337 void StdHexExp::v_LaplacianMatrixOp(const Array<OneD, const NekDouble> &inarray,
2338  Array<OneD, NekDouble> &outarray,
2339  const StdMatrixKey &mkey)
2340 {
2341  StdHexExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
2342 }
2343 
2344 void StdHexExp::v_LaplacianMatrixOp(const int k1, const int k2,
2345  const Array<OneD, const NekDouble> &inarray,
2346  Array<OneD, NekDouble> &outarray,
2347  const StdMatrixKey &mkey)
2348 {
2349  StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
2350 }
2351 
2352 void StdHexExp::v_WeakDerivMatrixOp(const int i,
2353  const Array<OneD, const NekDouble> &inarray,
2354  Array<OneD, NekDouble> &outarray,
2355  const StdMatrixKey &mkey)
2356 {
2357  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
2358 }
2359 
2360 void StdHexExp::v_HelmholtzMatrixOp(const Array<OneD, const NekDouble> &inarray,
2361  Array<OneD, NekDouble> &outarray,
2362  const StdMatrixKey &mkey)
2363 {
2364  StdHexExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
2365 }
2366 
2367 void StdHexExp::v_GeneralMatrixOp_MatOp(
2368  const Array<OneD, const NekDouble> &inarray,
2369  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
2370 {
2371  DNekMatSharedPtr mat = m_stdMatrixManager[mkey];
2372 
2373  if (inarray.get() == outarray.get())
2374  {
2375  Array<OneD, NekDouble> tmp(m_ncoeffs);
2376  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, tmp.get(), 1);
2377 
2378  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
2379  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
2380  }
2381  else
2382  {
2383  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
2384  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
2385  }
2386 }
2387 
2388 void StdHexExp::v_MultiplyByStdQuadratureMetric(
2389  const Array<OneD, const NekDouble> &inarray,
2390  Array<OneD, NekDouble> &outarray)
2391 {
2392  int nquad0 = m_base[0]->GetNumPoints();
2393  int nquad1 = m_base[1]->GetNumPoints();
2394  int nquad2 = m_base[2]->GetNumPoints();
2395  int nq01 = nquad0 * nquad1;
2396  int nq12 = nquad1 * nquad2;
2397 
2398  const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
2399  const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
2400  const Array<OneD, const NekDouble> &w2 = m_base[2]->GetW();
2401 
2402  for (int i = 0; i < nq12; ++i)
2403  {
2404  Vmath::Vmul(nquad0, inarray.get() + i * nquad0, 1, w0.get(), 1,
2405  outarray.get() + i * nquad0, 1);
2406  }
2407 
2408  for (int i = 0; i < nq12; ++i)
2409  {
2410  Vmath::Smul(nquad0, w1[i % nquad1], outarray.get() + i * nquad0, 1,
2411  outarray.get() + i * nquad0, 1);
2412  }
2413 
2414  for (int i = 0; i < nquad2; ++i)
2415  {
2416  Vmath::Smul(nq01, w2[i], outarray.get() + i * nq01, 1,
2417  outarray.get() + i * nq01, 1);
2418  }
2419 }
2420 
2421 void StdHexExp::v_SVVLaplacianFilter(Array<OneD, NekDouble> &array,
2422  const StdMatrixKey &mkey)
2423 {
2424  // Generate an orthonogal expansion
2425  int qa = m_base[0]->GetNumPoints();
2426  int qb = m_base[1]->GetNumPoints();
2427  int qc = m_base[2]->GetNumPoints();
2428  int nmodes_a = m_base[0]->GetNumModes();
2429  int nmodes_b = m_base[1]->GetNumModes();
2430  int nmodes_c = m_base[2]->GetNumModes();
2431  // Declare orthogonal basis.
2432  LibUtilities::PointsKey pa(qa, m_base[0]->GetPointsType());
2433  LibUtilities::PointsKey pb(qb, m_base[1]->GetPointsType());
2434  LibUtilities::PointsKey pc(qc, m_base[2]->GetPointsType());
2435 
2436  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodes_a, pa);
2437  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodes_b, pb);
2438  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A, nmodes_c, pc);
2439  StdHexExp OrthoExp(Ba, Bb, Bc);
2440 
2441  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2442  int cnt = 0;
2443 
2444  // project onto modal space.
2445  OrthoExp.FwdTrans(array, orthocoeffs);
2446 
2447  if (mkey.ConstFactorExists(eFactorSVVPowerKerDiffCoeff))
2448  {
2449  // Rodrigo's power kernel
2450  NekDouble cutoff = mkey.GetConstFactor(eFactorSVVCutoffRatio);
2451  NekDouble SvvDiffCoeff =
2452  mkey.GetConstFactor(eFactorSVVPowerKerDiffCoeff) *
2453  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2454 
2455  for (int i = 0; i < nmodes_a; ++i)
2456  {
2457  for (int j = 0; j < nmodes_b; ++j)
2458  {
2459  NekDouble fac1 = std::max(
2460  pow((1.0 * i) / (nmodes_a - 1), cutoff * nmodes_a),
2461  pow((1.0 * j) / (nmodes_b - 1), cutoff * nmodes_b));
2462 
2463  for (int k = 0; k < nmodes_c; ++k)
2464  {
2465  NekDouble fac =
2466  std::max(fac1, pow((1.0 * k) / (nmodes_c - 1),
2467  cutoff * nmodes_c));
2468 
2469  orthocoeffs[cnt] *= SvvDiffCoeff * fac;
2470  cnt++;
2471  }
2472  }
2473  }
2474  }
2475  else if (mkey.ConstFactorExists(
2476  eFactorSVVDGKerDiffCoeff)) // Rodrigo/Mansoor's DG Kernel
2477  {
2478  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff) *
2479  mkey.GetConstFactor(eFactorSVVDiffCoeff);
2480 
2481  int max_abc = max(nmodes_a - kSVVDGFiltermodesmin,
2482  nmodes_b - kSVVDGFiltermodesmin);
2483  max_abc = max(max_abc, nmodes_c - kSVVDGFiltermodesmin);
2484  // clamp max_abc
2485  max_abc = max(max_abc, 0);
2486  max_abc = min(max_abc, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
2487 
2488  for (int i = 0; i < nmodes_a; ++i)
2489  {
2490  for (int j = 0; j < nmodes_b; ++j)
2491  {
2492  int maxij = max(i, j);
2493 
2494  for (int k = 0; k < nmodes_c; ++k)
2495  {
2496  int maxijk = max(maxij, k);
2497  maxijk = min(maxijk, kSVVDGFiltermodesmax - 1);
2498 
2499  orthocoeffs[cnt] *=
2500  SvvDiffCoeff * kSVVDGFilter[max_abc][maxijk];
2501  cnt++;
2502  }
2503  }
2504  }
2505  }
2506  else
2507  {
2508 
2509  int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio) *
2510  min(nmodes_a, nmodes_b));
2511  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
2512  // Filter just trilinear space
2513  int nmodes = max(nmodes_a, nmodes_b);
2514  nmodes = max(nmodes, nmodes_c);
2515 
2516  Array<OneD, NekDouble> fac(nmodes, 1.0);
2517  for (int j = cutoff; j < nmodes; ++j)
2518  {
2519  fac[j] = fabs((j - nmodes) / ((NekDouble)(j - cutoff + 1.0)));
2520  fac[j] *= fac[j]; // added this line to conform with equation
2521  }
2522 
2523  for (int i = 0; i < nmodes_a; ++i)
2524  {
2525  for (int j = 0; j < nmodes_b; ++j)
2526  {
2527  for (int k = 0; k < nmodes_c; ++k)
2528  {
2529  if ((i >= cutoff) || (j >= cutoff) || (k >= cutoff))
2530  {
2531  orthocoeffs[i * nmodes_a * nmodes_b + j * nmodes_c +
2532  k] *=
2533  (SvvDiffCoeff * exp(-(fac[i] + fac[j] + fac[k])));
2534  }
2535  else
2536  {
2537  orthocoeffs[i * nmodes_a * nmodes_b + j * nmodes_c +
2538  k] *= 0.0;
2539  }
2540  }
2541  }
2542  }
2543  }
2544 
2545  // backward transform to physical space
2546  OrthoExp.BwdTrans(orthocoeffs, array);
2547 }
2548 
2549 void StdHexExp::v_ExponentialFilter(Array<OneD, NekDouble> &array,
2550  const NekDouble alpha,
2551  const NekDouble exponent,
2552  const NekDouble cutoff)
2553 {
2554  // Generate an orthogonal expansion
2555  int qa = m_base[0]->GetNumPoints();
2556  int qb = m_base[1]->GetNumPoints();
2557  int qc = m_base[2]->GetNumPoints();
2558  int nmodesA = m_base[0]->GetNumModes();
2559  int nmodesB = m_base[1]->GetNumModes();
2560  int nmodesC = m_base[2]->GetNumModes();
2561  int P = nmodesA - 1;
2562  int Q = nmodesB - 1;
2563  int R = nmodesC - 1;
2564 
2565  // Declare orthogonal basis.
2566  LibUtilities::PointsKey pa(qa, m_base[0]->GetPointsType());
2567  LibUtilities::PointsKey pb(qb, m_base[1]->GetPointsType());
2568  LibUtilities::PointsKey pc(qc, m_base[2]->GetPointsType());
2569 
2570  LibUtilities::BasisKey Ba(LibUtilities::eOrtho_A, nmodesA, pa);
2571  LibUtilities::BasisKey Bb(LibUtilities::eOrtho_A, nmodesB, pb);
2572  LibUtilities::BasisKey Bc(LibUtilities::eOrtho_A, nmodesC, pc);
2573  StdHexExp OrthoExp(Ba, Bb, Bc);
2574 
2575  // Cutoff
2576  int Pcut = cutoff * P;
2577  int Qcut = cutoff * Q;
2578  int Rcut = cutoff * R;
2579 
2580  // Project onto orthogonal space.
2581  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
2582  OrthoExp.FwdTrans(array, orthocoeffs);
2583 
2584  //
2585  NekDouble fac, fac1, fac2, fac3;
2586  int index = 0;
2587  for (int i = 0; i < nmodesA; ++i)
2588  {
2589  for (int j = 0; j < nmodesB; ++j)
2590  {
2591  for (int k = 0; k < nmodesC; ++k, ++index)
2592  {
2593  // to filter out only the "high-modes"
2594  if (i > Pcut || j > Qcut || k > Rcut)
2595  {
2596  fac1 = (NekDouble)(i - Pcut) / ((NekDouble)(P - Pcut));
2597  fac2 = (NekDouble)(j - Qcut) / ((NekDouble)(Q - Qcut));
2598  fac3 = (NekDouble)(k - Rcut) / ((NekDouble)(R - Rcut));
2599  fac = max(max(fac1, fac2), fac3);
2600  fac = pow(fac, exponent);
2601  orthocoeffs[index] *= exp(-alpha * fac);
2602  }
2603  }
2604  }
2605  }
2606 
2607  // backward transform to physical space
2608  OrthoExp.BwdTrans(orthocoeffs, array);
2609 }
2610 
2611 void StdHexExp::v_GetSimplexEquiSpacedConnectivity(Array<OneD, int> &conn,
2612  bool standard)
2613 {
2614  boost::ignore_unused(standard);
2615 
2616  int np0 = m_base[0]->GetNumPoints();
2617  int np1 = m_base[1]->GetNumPoints();
2618  int np2 = m_base[2]->GetNumPoints();
2619  int np = max(np0, max(np1, np2));
2620 
2621  conn = Array<OneD, int>(6 * (np - 1) * (np - 1) * (np - 1));
2622 
2623  int row = 0;
2624  int rowp1 = 0;
2625  int cnt = 0;
2626  int plane = 0;
2627  for (int i = 0; i < np - 1; ++i)
2628  {
2629  for (int j = 0; j < np - 1; ++j)
2630  {
2631  rowp1 += np;
2632  for (int k = 0; k < np - 1; ++k)
2633  {
2634  conn[cnt++] = plane + row + k;
2635  conn[cnt++] = plane + row + k + 1;
2636  conn[cnt++] = plane + rowp1 + k;
2637 
2638  conn[cnt++] = plane + rowp1 + k + 1;
2639  conn[cnt++] = plane + rowp1 + k;
2640  conn[cnt++] = plane + row + k + 1;
2641  }
2642  row += np;
2643  }
2644  plane += np * np;
2645  }
2646 }
2647 } // namespace StdRegions
2648 } // namespace Nektar
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
ASSERTL1
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
ASSERTL2
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272
Nektar::LibUtilities::BasisKey
Describes the specification for a Basis.
Definition: Basis.h:50
Nektar::LibUtilities::PointsKey
Defines a specification for a set of points.
Definition: Points.h:59
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:168
Nektar::StdRegions::StdExpansion3D::v_LaplacianMatrixOp_MatFree
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: StdExpansion3D.cpp:286
Nektar::StdRegions::StdExpansion3D::BwdTrans_SumFacKernel
void BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdExpansion3D.cpp:109
Nektar::StdRegions::StdExpansion3D::IProductWRTBase_SumFacKernel
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdExpansion3D.cpp:121
Nektar::StdRegions::StdExpansion3D::GetEdgeNcoeffs
int GetEdgeNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th edge.
Definition: StdExpansion3D.h:140
Nektar::StdRegions::StdExpansion3D::v_HelmholtzMatrixOp_MatFree
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: StdExpansion3D.cpp:329
Nektar::StdRegions::StdExpansion3D::PhysTensorDeriv
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d1, Array< OneD, NekDouble > &outarray_d2, Array< OneD, NekDouble > &outarray_d3)
Calculate the 3D derivative in the local tensor/collapsed coordinate at the physical points.
Definition: StdExpansion3D.cpp:67
Nektar::StdRegions::StdExpansion
The base class for all shapes.
Definition: StdExpansion.h:71
Nektar::StdRegions::StdExpansion::GetNcoeffs
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
Nektar::StdRegions::StdExpansion::GetTotPoints
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
Nektar::StdRegions::StdExpansion::WeakDerivMatrixOp_MatFree
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.cpp:818
Nektar::StdRegions::StdExpansion::GetBasisType
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:163
Nektar::StdRegions::StdExpansion::GetStdMatrix
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:611
Nektar::StdRegions::StdExpansion::LocCoordToLocCollapsed
void LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Convert local cartesian coordinate xi into local collapsed coordinates eta.
Definition: StdExpansion.h:974
Nektar::StdRegions::StdExpansion::BwdTrans_SumFac
void BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:1179
Nektar::StdRegions::StdExpansion::IProductWRTDerivBase_SumFac
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:1185
Nektar::StdRegions::StdExpansion::CreateGeneralMatrix
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
Definition: StdExpansion.cpp:253
Nektar::StdRegions::StdExpansion::IProductWRTBase
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
Definition: StdExpansion.h:536
Nektar::StdRegions::StdExpansion::LaplacianMatrixOp_MatFree
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:1217
Nektar::StdRegions::StdExpansion::DetShapeType
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:375
Nektar::StdRegions::StdExpansion::BwdTrans
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Backward transformation from coefficient space to physical space.
Definition: StdExpansion.h:432
Nektar::StdRegions::StdExpansion::GetPointsType
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:213
Nektar::StdRegions::StdExpansion::FwdTrans
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
Definition: StdExpansion.h:1663
Nektar::StdRegions::StdExpansion::GetNumPoints
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:226
Nektar::StdRegions::StdExpansion::MultiplyByQuadratureMetric
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:731
Nektar::StdRegions::StdExpansion::m_stdMatrixManager
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
Definition: StdExpansion.h:1150
Nektar::StdRegions::StdExpansion::GetBasisNumModes
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:176
Nektar::StdRegions::StdExpansion::PhysDeriv
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:850
Nektar::StdRegions::StdExpansion::m_base
Array< OneD, LibUtilities::BasisSharedPtr > m_base
Definition: StdExpansion.h:1145
Nektar::StdRegions::StdExpansion::MassMatrixOp_MatFree
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.cpp:666
Nektar::StdRegions::StdHexExp
Class representing a hexehedral element in reference space.
Definition: StdHexExp.h:48
Nektar::StdRegions::StdHexExp::v_IProductWRTBase_MatOp
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:352
Nektar::StdRegions::StdHexExp::v_ExponentialFilter
virtual void v_ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
Definition: StdHexExp.cpp:2549
Nektar::StdRegions::StdHexExp::v_PhysDeriv
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Differentiation Methods.
Definition: StdHexExp.cpp:87
Nektar::StdRegions::StdHexExp::v_NumBndryCoeffs
virtual int v_NumBndryCoeffs() const
Definition: StdHexExp.cpp:670
Nektar::StdRegions::StdHexExp::v_HelmholtzMatrixOp
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2360
Nektar::StdRegions::StdHexExp::v_GetTraceNumModes
virtual void v_GetTraceNumModes(const int fid, int &numModes0, int &numModes1, Orientation faceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdHexExp.cpp:854
Nektar::StdRegions::StdHexExp::v_SVVLaplacianFilter
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2421
Nektar::StdRegions::StdHexExp::v_PhysEvaluateBasis
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdHexExp.cpp:629
Nektar::StdRegions::StdHexExp::v_GetTraceBasisKey
virtual const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int k) const
Definition: StdHexExp.cpp:798
Nektar::StdRegions::StdHexExp::v_GetTraceNumPoints
virtual int v_GetTraceNumPoints(const int i) const
Definition: StdHexExp.cpp:750
Nektar::StdRegions::StdHexExp::v_LocCoordToLocCollapsed
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdHexExp.cpp:565
Nektar::StdRegions::StdHexExp::v_FillMode
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:584
Nektar::StdRegions::StdHexExp::v_GetTraceNcoeffs
virtual int v_GetTraceNcoeffs(const int i) const
Definition: StdHexExp.cpp:709
Nektar::StdRegions::StdHexExp::v_BwdTrans_SumFac
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:195
Nektar::StdRegions::StdHexExp::v_GetEdgeNcoeffs
virtual int v_GetEdgeNcoeffs(const int i) const
Definition: StdHexExp.cpp:2219
Nektar::StdRegions::StdHexExp::v_GetBoundaryMap
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdHexExp.cpp:1091
Nektar::StdRegions::StdHexExp::v_IProductWRTBase_SumFac
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multbyweights=true)
Definition: StdHexExp.cpp:367
Nektar::StdRegions::StdHexExp::v_GenMatrix
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2237
Nektar::StdRegions::StdHexExp::v_GetTraceIntNcoeffs
virtual int v_GetTraceIntNcoeffs(const int i) const
Definition: StdHexExp.cpp:726
Nektar::StdRegions::StdHexExp::v_StdPhysDeriv
virtual void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
Definition: StdHexExp.cpp:133
Nektar::StdRegions::StdHexExp::v_GeneralMatrixOp_MatOp
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2367
Nektar::StdRegions::StdHexExp::v_IProductWRTDerivBase_SumFac
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:519
Nektar::StdRegions::StdHexExp::v_GetVertexMap
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdHexExp.cpp:902
Nektar::StdRegions::StdHexExp::v_GetNverts
virtual int v_GetNverts() const
Definition: StdHexExp.cpp:650
Nektar::StdRegions::StdHexExp::v_GetTraceInteriorToElementMap
virtual void v_GetTraceInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdHexExp.cpp:1895
Nektar::StdRegions::StdHexExp::v_LocCollapsedToLocCoord
virtual void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
Definition: StdHexExp.cpp:573
Nektar::StdRegions::StdHexExp::v_GetTraceCoeffMap
virtual void v_GetTraceCoeffMap(const unsigned int fid, Array< OneD, unsigned int > &maparray)
Definition: StdHexExp.cpp:1183
Nektar::StdRegions::StdHexExp::v_GetTracePointsKey
virtual LibUtilities::PointsKey v_GetTracePointsKey(const int i, const int j) const
Definition: StdHexExp.cpp:768
Nektar::StdRegions::StdHexExp::v_MassMatrixOp
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2330
Nektar::StdRegions::StdHexExp::v_IProductWRTBase_SumFacKernel
virtual void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdHexExp.cpp:401
Nektar::StdRegions::StdHexExp::v_GetTotalTraceIntNcoeffs
virtual int v_GetTotalTraceIntNcoeffs() const
Definition: StdHexExp.cpp:743
Nektar::StdRegions::StdHexExp::v_GetNedges
virtual int v_GetNedges() const
Definition: StdHexExp.cpp:655
Nektar::StdRegions::StdHexExp::v_GetCoords
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
Definition: StdHexExp.cpp:826
Nektar::StdRegions::StdHexExp::v_LaplacianMatrixOp
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2337
Nektar::StdRegions::StdHexExp::v_IsBoundaryInteriorExpansion
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdHexExp.cpp:70
Nektar::StdRegions::StdHexExp::v_BwdTrans
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:168
Nektar::StdRegions::StdHexExp::v_CreateStdMatrix
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2325
Nektar::StdRegions::StdHexExp::v_MultiplyByStdQuadratureMetric
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:2388
Nektar::StdRegions::StdHexExp::v_GetInteriorMap
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdHexExp.cpp:1030
Nektar::StdRegions::StdHexExp::v_NumDGBndryCoeffs
virtual int v_NumDGBndryCoeffs() const
Definition: StdHexExp.cpp:690
Nektar::StdRegions::StdHexExp::v_CalcNumberOfCoefficients
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdHexExp.cpp:788
Nektar::StdRegions::StdHexExp::v_GetNtraces
virtual int v_GetNtraces() const
Definition: StdHexExp.cpp:660
Nektar::StdRegions::StdHexExp::v_FwdTrans
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:277
Nektar::StdRegions::StdHexExp::v_GetSimplexEquiSpacedConnectivity
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: StdHexExp.cpp:2611
Nektar::StdRegions::StdHexExp::v_BwdTrans_SumFacKernel
virtual void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
Definition: StdHexExp.cpp:220
Nektar::StdRegions::StdHexExp::v_DetShapeType
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdHexExp.cpp:665
Nektar::StdRegions::StdHexExp::v_WeakDerivMatrixOp
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdHexExp.cpp:2352
Nektar::StdRegions::StdHexExp::v_IProductWRTDerivBase_MatOp
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:489
Nektar::StdRegions::StdHexExp::v_GetElmtTraceToTraceMap
virtual void v_GetElmtTraceToTraceMap(const unsigned int fid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation faceOrient, int P, int Q)
Definition: StdHexExp.cpp:1307
Nektar::StdRegions::StdHexExp::v_GetEdgeInteriorToElementMap
virtual void v_GetEdgeInteriorToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation traceOrient=eDir1FwdDir1_Dir2FwdDir2)
Definition: StdHexExp.cpp:1602
Nektar::StdRegions::StdHexExp::v_IProductWRTDerivBase
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:482
Nektar::StdRegions::StdHexExp::v_IProductWRTBase
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdHexExp.cpp:335
Nektar::StdRegions::StdMatrixKey::GetMatrixType
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
Nektar::StdRegions::StdMatrixKey::GetConstFactor
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:126
Nektar::StdRegions::StdMatrixKey::ConstFactorExists
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:135
Blas::Dgemv
static void Dgemv(const char &trans, const int &m, const int &n, const double &alpha, const double *a, const int &lda, const double *x, const int &incx, const double &beta, double *y, const int &incy)
BLAS level 2: Matrix vector multiply y = A x where A[m x n].
Definition: Blas.hpp:246
Blas::Dscal
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:168
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:368
CellMLToNektar.cellml_metadata.p
p
Definition: cellml_metadata.py:350
Nektar::LibUtilities::StdHexData::getNumberOfCoefficients
int getNumberOfCoefficients(int Na, int Nb, int Nc)
Definition: ShapeType.hpp:156
Nektar::LibUtilities::eModified_B
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:51
Nektar::LibUtilities::eOrtho_A
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:44
Nektar::LibUtilities::eModified_C
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:52
Nektar::LibUtilities::eGLL_Lagrange
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:58
Nektar::LibUtilities::eOrtho_C
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:48
Nektar::LibUtilities::eOrtho_B
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
Nektar::NekConstants::kNekZeroTol
static const NekDouble kNekZeroTol
Definition: NektarUnivConsts.hpp:51
Nektar::StdRegions::kSVVDGFiltermodesmin
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:352
Nektar::StdRegions::kSVVDGFiltermodesmax
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:353
Nektar::StdRegions::EvaluateQuadFaceBasisKey
LibUtilities::BasisKey EvaluateQuadFaceBasisKey(const int facedir, const LibUtilities::BasisType faceDirBasisType, const int numpoints, const int nummodes)
Definition: StdExpansion3D.cpp:443
Nektar::StdRegions::kSVVDGFilter
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:355
Nektar
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
Nektar::NullNekDouble1DArray
static Array< OneD, NekDouble > NullNekDouble1DArray
Definition: SharedArray.hpp:883
Nektar::DNekMatSharedPtr
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmul
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.cpp:209
Vmath::Smul
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.cpp:248
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
Generated on Tue Aug 23 2022 07:27:15 for Nektar++ by doxygen 1.9.1