Nektar++
StdQuadExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdQuadExp.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: Quadrilateral routines built upon StdExpansion2D
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
39 #include <StdRegions/StdQuadExp.h>
40 #include <StdRegions/StdSegExp.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46 namespace StdRegions
47 {
48 
49 StdQuadExp::StdQuadExp()
50 {
51 }
52 
53 /** \brief Constructor using BasisKey class for quadrature
54  * points and order definition
55  */
56 StdQuadExp::StdQuadExp(const LibUtilities::BasisKey &Ba,
57  const LibUtilities::BasisKey &Bb)
58  : StdExpansion(Ba.GetNumModes() * Bb.GetNumModes(), 2, Ba, Bb),
59  StdExpansion2D(Ba.GetNumModes() * Bb.GetNumModes(), Ba, Bb)
60 {
61 }
62 
63 /** \brief Copy Constructor */
65 {
66 }
67 
68 /** \brief Destructor */
70 {
71 }
72 
73 /////////////////////////
74 // Integration Methods //
75 /////////////////////////
76 
78 {
79  Array<OneD, const NekDouble> w0 = m_base[0]->GetW();
80  Array<OneD, const NekDouble> w1 = m_base[1]->GetW();
81 
82  return StdExpansion2D::Integral(inarray, w0, w1);
83 }
84 
85 /////////////////////////////
86 // Differentiation Methods //
87 /////////////////////////////
88 
89 /** \brief Calculate the derivative of the physical points
90  *
91  * For quadrilateral region can use the Tensor_Deriv function
92  * defined under StdExpansion.
93  */
94 
96  Array<OneD, NekDouble> &out_d0,
97  Array<OneD, NekDouble> &out_d1,
98  Array<OneD, NekDouble> &out_d2)
99 {
100  boost::ignore_unused(out_d2);
101  PhysTensorDeriv(inarray, out_d0, out_d1);
102 }
103 
104 void StdQuadExp::v_PhysDeriv(const int dir,
105  const Array<OneD, const NekDouble> &inarray,
106  Array<OneD, NekDouble> &outarray)
107 {
108  switch (dir)
109  {
110  case 0:
111  {
112  PhysTensorDeriv(inarray, outarray, NullNekDouble1DArray);
113  }
114  break;
115  case 1:
116  {
117  PhysTensorDeriv(inarray, NullNekDouble1DArray, outarray);
118  }
119  break;
120  default:
121  {
122  ASSERTL1(false, "input dir is out of range");
123  }
124  break;
125  }
126 }
127 
129  Array<OneD, NekDouble> &out_d0,
130  Array<OneD, NekDouble> &out_d1,
131  Array<OneD, NekDouble> &out_d2)
132 {
133  boost::ignore_unused(out_d2);
134  // PhysTensorDeriv(inarray, out_d0, out_d1);
135  StdQuadExp::v_PhysDeriv(inarray, out_d0, out_d1);
136 }
137 
138 void StdQuadExp::v_StdPhysDeriv(const int dir,
139  const Array<OneD, const NekDouble> &inarray,
140  Array<OneD, NekDouble> &outarray)
141 {
142  // PhysTensorDeriv(inarray, outarray);
143  StdQuadExp::v_PhysDeriv(dir, inarray, outarray);
144 }
145 
146 ////////////////
147 // Transforms //
148 ////////////////
149 
151  Array<OneD, NekDouble> &outarray)
152 {
153  if (m_base[0]->Collocation() && m_base[1]->Collocation())
154  {
156  inarray, 1, outarray, 1);
157  }
158  else
159  {
160  StdQuadExp::v_BwdTrans_SumFac(inarray, outarray);
161  }
162 }
163 
165  Array<OneD, NekDouble> &outarray)
166 {
168  m_base[1]->GetNumModes());
169 
170  BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(), inarray,
171  outarray, wsp, true, true);
172 }
173 
174 // The arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify
175 // whether to check if the basis has collocation properties (i.e. for the
176 // classical spectral element basis, In this case the 1D 'B' matrix is equal to
177 // the identity matrix which can be exploited to speed up the calculations).
178 // However, as this routine also allows to pass the matrix 'DB' (derivative of
179 // the basis), the collocation property cannot always be used. Therefor follow
180 // this rule: if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
181 // base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
182 // base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
183 // base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
185  const Array<OneD, const NekDouble> &base0,
186  const Array<OneD, const NekDouble> &base1,
187  const Array<OneD, const NekDouble> &inarray,
189  bool doCheckCollDir0, bool doCheckCollDir1)
190 {
191  int nquad0 = m_base[0]->GetNumPoints();
192  int nquad1 = m_base[1]->GetNumPoints();
193  int nmodes0 = m_base[0]->GetNumModes();
194  int nmodes1 = m_base[1]->GetNumModes();
195 
196  bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
197  bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
198 
199  if (colldir0 && colldir1)
200  {
201  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, outarray.get(), 1);
202  }
203  else if (colldir0)
204  {
205  Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes1, 1.0, &inarray[0], nquad0,
206  base1.get(), nquad1, 0.0, &outarray[0], nquad0);
207  }
208  else if (colldir1)
209  {
210  Blas::Dgemm('N', 'N', nquad0, nmodes1, nmodes0, 1.0, base0.get(),
211  nquad0, &inarray[0], nmodes0, 0.0, &outarray[0], nquad0);
212  }
213  else
214  {
215  ASSERTL1(wsp.size() >= nquad0 * nmodes1,
216  "Workspace size is not sufficient");
217 
218  // Those two calls correpsond to the operation
219  // out = B0*in*Transpose(B1);
220  Blas::Dgemm('N', 'N', nquad0, nmodes1, nmodes0, 1.0, base0.get(),
221  nquad0, &inarray[0], nmodes0, 0.0, &wsp[0], nquad0);
222  Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes1, 1.0, &wsp[0], nquad0,
223  base1.get(), nquad1, 0.0, &outarray[0], nquad0);
224  }
225 }
226 
228  Array<OneD, NekDouble> &outarray)
229 {
230  if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()))
231  {
232  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
233  }
234  else
235  {
236  StdQuadExp::v_IProductWRTBase(inarray, outarray);
237 
238  // get Mass matrix inverse
239  StdMatrixKey masskey(eInvMass, DetShapeType(), *this);
240  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
241 
242  // copy inarray in case inarray == outarray
243  NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
244  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
245 
246  out = (*matsys) * in;
247  }
248 }
249 
251  const Array<OneD, const NekDouble> &inarray,
252  Array<OneD, NekDouble> &outarray)
253 {
254  if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()))
255  {
256  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
257  }
258  else
259  {
260  int i, j;
261  int npoints[2] = {m_base[0]->GetNumPoints(), m_base[1]->GetNumPoints()};
262  int nmodes[2] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes()};
263 
264  fill(outarray.get(), outarray.get() + m_ncoeffs, 0.0);
265 
266  Array<OneD, NekDouble> physEdge[4];
267  Array<OneD, NekDouble> coeffEdge[4];
268  for (i = 0; i < 4; i++)
269  {
270  physEdge[i] = Array<OneD, NekDouble>(npoints[i % 2]);
271  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i % 2]);
272  }
273 
274  for (i = 0; i < npoints[0]; i++)
275  {
276  physEdge[0][i] = inarray[i];
277  physEdge[2][i] = inarray[npoints[0] * npoints[1] - 1 - i];
278  }
279 
280  for (i = 0; i < npoints[1]; i++)
281  {
282  physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
283  physEdge[3][i] =
284  inarray[(npoints[1] - 1) * npoints[0] - i * npoints[0]];
285  }
286 
287  StdSegExpSharedPtr segexp[2] = {
289  m_base[0]->GetBasisKey()),
291  m_base[1]->GetBasisKey())};
292 
293  Array<OneD, unsigned int> mapArray;
294  Array<OneD, int> signArray;
295  NekDouble sign;
296 
297  for (i = 0; i < 4; i++)
298  {
299  segexp[i % 2]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
300 
301  GetTraceToElementMap(i, mapArray, signArray);
302  for (j = 0; j < nmodes[i % 2]; j++)
303  {
304  sign = (NekDouble)signArray[j];
305  outarray[mapArray[j]] = sign * coeffEdge[i][j];
306  }
307  }
308 
311 
312  StdMatrixKey masskey(eMass, DetShapeType(), *this);
313  MassMatrixOp(outarray, tmp0, masskey);
314  IProductWRTBase(inarray, tmp1);
315 
316  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
317 
318  // get Mass matrix inverse (only of interior DOF)
319  // use block (1,1) of the static condensed system
320  // note: this block alreay contains the inverse matrix
321  DNekMatSharedPtr matsys =
322  (m_stdStaticCondMatrixManager[masskey])->GetBlock(1, 1);
323 
324  int nBoundaryDofs = NumBndryCoeffs();
325  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
326 
327  Array<OneD, NekDouble> rhs(nInteriorDofs);
328  Array<OneD, NekDouble> result(nInteriorDofs);
329 
330  GetInteriorMap(mapArray);
331 
332  for (i = 0; i < nInteriorDofs; i++)
333  {
334  rhs[i] = tmp1[mapArray[i]];
335  }
336 
337  Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, 1.0,
338  &(matsys->GetPtr())[0], nInteriorDofs, rhs.get(), 1, 0.0,
339  result.get(), 1);
340 
341  for (i = 0; i < nInteriorDofs; i++)
342  {
343  outarray[mapArray[i]] = result[i];
344  }
345  }
346 }
347 
348 /////////////////////////////
349 // Inner Product Functions //
350 /////////////////////////////
351 
352 /** \brief Calculate the inner product of inarray with respect to
353  * the basis B=base0*base1 and put into outarray
354  *
355  * \f$
356  * \begin{array}{rcl}
357  * I_{pq} = (\phi_p \phi_q, u) & = & \sum_{i=0}^{nq_0}
358  * \sum_{j=0}^{nq_1}
359  * \phi_p(\xi_{0,i}) \phi_q(\xi_{1,j}) w^0_i w^1_j u(\xi_{0,i}
360  * \xi_{1,j}) \\
361  * & = & \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i})
362  * \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j}) \tilde{u}_{i,j}
363  * \end{array}
364  * \f$
365  *
366  * where
367  *
368  * \f$ \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
369  *
370  * which can be implemented as
371  *
372  * \f$ f_{qi} = \sum_{j=0}^{nq_1} \phi_q(\xi_{1,j})
373  * \tilde{u}_{i,j} = {\bf B_1 U} \f$
374  * \f$ I_{pq} = \sum_{i=0}^{nq_0} \phi_p(\xi_{0,i}) f_{qi} =
375  * {\bf B_0 F} \f$
376  */
378  Array<OneD, NekDouble> &outarray)
379 {
380  if (m_base[0]->Collocation() && m_base[1]->Collocation())
381  {
382  MultiplyByQuadratureMetric(inarray, outarray);
383  }
384  else
385  {
386  StdQuadExp::v_IProductWRTBase_SumFac(inarray, outarray);
387  }
388 }
389 
391  const Array<OneD, const NekDouble> &inarray,
392  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
393 {
394  int nquad0 = m_base[0]->GetNumPoints();
395  int nquad1 = m_base[1]->GetNumPoints();
396  int order0 = m_base[0]->GetNumModes();
397 
398  if (multiplybyweights)
399  {
400  Array<OneD, NekDouble> tmp(nquad0 * nquad1 + nquad1 * order0);
401  Array<OneD, NekDouble> wsp(tmp + nquad0 * nquad1);
402 
403  // multiply by integration constants
404  MultiplyByQuadratureMetric(inarray, tmp);
405  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
406  m_base[1]->GetBdata(), tmp, outarray, wsp,
407  true, true);
408  }
409  else
410  {
411  Array<OneD, NekDouble> wsp(nquad1 * order0);
412  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
413  m_base[1]->GetBdata(), inarray, outarray,
414  wsp, true, true);
415  }
416 }
417 
419  const Array<OneD, const NekDouble> &inarray,
420  Array<OneD, NekDouble> &outarray)
421 {
422  int nq = GetTotPoints();
423  StdMatrixKey iprodmatkey(eIProductWRTBase, DetShapeType(), *this);
424  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
425 
426  Blas::Dgemv('N', m_ncoeffs, nq, 1.0, iprodmat->GetPtr().get(), m_ncoeffs,
427  inarray.get(), 1, 0.0, outarray.get(), 1);
428 }
429 
431  const int dir, const Array<OneD, const NekDouble> &inarray,
432  Array<OneD, NekDouble> &outarray)
433 {
434  v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
435 }
436 
438  const int dir, const Array<OneD, const NekDouble> &inarray,
439  Array<OneD, NekDouble> &outarray)
440 {
441  ASSERTL0((dir == 0) || (dir == 1), "input dir is out of range");
442 
443  int nquad0 = m_base[0]->GetNumPoints();
444  int nquad1 = m_base[1]->GetNumPoints();
445  int nqtot = nquad0 * nquad1;
446  int order0 = m_base[0]->GetNumModes();
447 
448  Array<OneD, NekDouble> tmp(nqtot + nquad1 * order0);
449  Array<OneD, NekDouble> wsp(tmp + nqtot);
450 
451  // multiply by integration constants
452  MultiplyByQuadratureMetric(inarray, tmp);
453 
454  if (dir) // dir == 1
455  {
456  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
457  m_base[1]->GetDbdata(), tmp, outarray, wsp,
458  true, false);
459  }
460  else // dir == 0
461  {
462  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
463  m_base[1]->GetBdata(), tmp, outarray, wsp,
464  false, true);
465  }
466 }
467 
469  const int dir, const Array<OneD, const NekDouble> &inarray,
470  Array<OneD, NekDouble> &outarray)
471 {
472  ASSERTL0((dir == 0) || (dir == 1), "input dir is out of range");
473 
474  int nq = GetTotPoints();
475  MatrixType mtype;
476 
477  if (dir) // dir == 1
478  {
479  mtype = eIProductWRTDerivBase1;
480  }
481  else // dir == 0
482  {
483  mtype = eIProductWRTDerivBase0;
484  }
485 
486  StdMatrixKey iprodmatkey(mtype, DetShapeType(), *this);
487  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
488 
489  Blas::Dgemv('N', m_ncoeffs, nq, 1.0, iprodmat->GetPtr().get(), m_ncoeffs,
490  inarray.get(), 1, 0.0, outarray.get(), 1);
491 }
492 
493 // the arguments doCheckCollDir0 and doCheckCollDir1 allow you to specify
494 // whether to check if the basis has collocation properties (i.e. for the
495 // classical spectral element basis, In this case the 1D 'B' matrix is equal to
496 // the identity matrix which can be exploited to speed up the calculations).
497 // However, as this routine also allows to pass the matrix 'DB' (derivative of
498 // the basis), the collocation property cannot always be used. Therefor follow
499 // this rule: if base0 == m_base[0]->GetBdata() --> set doCheckCollDir0 == true;
500 // base1 == m_base[1]->GetBdata() --> set doCheckCollDir1 == true;
501 // base0 == m_base[0]->GetDbdata() --> set doCheckCollDir0 == false;
502 // base1 == m_base[1]->GetDbdata() --> set doCheckCollDir1 == false;
504  const Array<OneD, const NekDouble> &base0,
505  const Array<OneD, const NekDouble> &base1,
506  const Array<OneD, const NekDouble> &inarray,
508  bool doCheckCollDir0, bool doCheckCollDir1)
509 {
510  int nquad0 = m_base[0]->GetNumPoints();
511  int nquad1 = m_base[1]->GetNumPoints();
512  int nmodes0 = m_base[0]->GetNumModes();
513  int nmodes1 = m_base[1]->GetNumModes();
514 
515  bool colldir0 = doCheckCollDir0 ? (m_base[0]->Collocation()) : false;
516  bool colldir1 = doCheckCollDir1 ? (m_base[1]->Collocation()) : false;
517 
518  if (colldir0 && colldir1)
519  {
520  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, outarray.get(), 1);
521  }
522  else if (colldir0)
523  {
524  Blas::Dgemm('N', 'N', nmodes0, nmodes1, nquad1, 1.0, inarray.get(),
525  nmodes0, base1.get(), nquad1, 0.0, outarray.get(), nmodes0);
526  }
527  else if (colldir1)
528  {
529  Blas::Dgemm('T', 'N', nmodes0, nquad1, nquad0, 1.0, base0.get(), nquad0,
530  inarray.get(), nquad0, 0.0, outarray.get(), nmodes0);
531  }
532  else
533  {
534  ASSERTL1(wsp.size() >= nquad1 * nmodes0,
535  "Workspace size is not sufficient");
536 
537 #if 1
538  Blas::Dgemm('T', 'N', nmodes0, nquad1, nquad0, 1.0, base0.get(), nquad0,
539  inarray.get(), nquad0, 0.0, wsp.get(), nmodes0);
540 
541 #else
542  for (int i = 0; i < nmodes0; ++i)
543  {
544  for (int j = 0; j < nquad1; ++j)
545  {
546  wsp[j * nmodes0 + i] =
547  Blas::Ddot(nquad0, base0.get() + i * nquad0, 1,
548  inarray.get() + j * nquad0, 1);
549  }
550  }
551 #endif
552  Blas::Dgemm('N', 'N', nmodes0, nmodes1, nquad1, 1.0, wsp.get(), nmodes0,
553  base1.get(), nquad1, 0.0, outarray.get(), nmodes0);
554  }
555 }
556 
557 //////////////////////////
558 // Evaluation functions //
559 //////////////////////////
560 
563 {
564  eta[0] = xi[0];
565  eta[1] = xi[1];
566 }
567 
570 {
571  xi[0] = eta[0];
572  xi[1] = eta[1];
573 }
574 
575 /** \brief Fill outarray with mode \a mode of expansion
576  *
577  * Note for quadrilateral expansions _base[0] (i.e. p) modes run
578  * fastest
579  */
580 
581 void StdQuadExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
582 {
583  int i;
584  int nquad0 = m_base[0]->GetNumPoints();
585  int nquad1 = m_base[1]->GetNumPoints();
586  Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
587  Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
588  int btmp0 = m_base[0]->GetNumModes();
589  int mode0 = mode % btmp0;
590  int mode1 = mode / btmp0;
591 
592  ASSERTL2(mode1 == (int)floor((1.0 * mode) / btmp0),
593  "Integer Truncation not Equiv to Floor");
594 
595  ASSERTL2(m_ncoeffs > mode,
596  "calling argument mode is larger than total expansion order");
597 
598  for (i = 0; i < nquad1; ++i)
599  {
600  Vmath::Vcopy(nquad0, (NekDouble *)(base0.get() + mode0 * nquad0), 1,
601  &outarray[0] + i * nquad0, 1);
602  }
603 
604  for (i = 0; i < nquad0; ++i)
605  {
606  Vmath::Vmul(nquad1, (NekDouble *)(base1.get() + mode1 * nquad1), 1,
607  &outarray[0] + i, nquad0, &outarray[0] + i, nquad0);
608  }
609 }
610 
611 //////////////////////
612 // Helper functions //
613 //////////////////////
614 
616 {
617  return 4;
618 }
619 
621 {
622  return 4;
623 }
624 
625 int StdQuadExp::v_GetTraceNcoeffs(const int i) const
626 {
627  ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
628 
629  if ((i == 0) || (i == 2))
630  {
631  return GetBasisNumModes(0);
632  }
633  else
634  {
635  return GetBasisNumModes(1);
636  }
637 }
638 
639 int StdQuadExp::v_GetTraceNumPoints(const int i) const
640 {
641  ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
642 
643  if ((i == 0) || (i == 2))
644  {
645  return GetNumPoints(0);
646  }
647  else
648  {
649  return GetNumPoints(1);
650  }
651 }
652 
654  const int j) const
655 {
656  boost::ignore_unused(j);
657  ASSERTL2((i >= 0) && (i <= 3), "edge id is out of range");
658 
659  if ((i == 0) || (i == 2))
660  {
661  return GetBasis(0)->GetBasisKey();
662  }
663  else
664  {
665  return GetBasis(1)->GetBasisKey();
666  }
667 }
668 
670 {
672 }
673 
675 {
679  "BasisType is not a boundary interior form");
683  "BasisType is not a boundary interior form");
684 
685  return 4 + 2 * (GetBasisNumModes(0) - 2) + 2 * (GetBasisNumModes(1) - 2);
686 }
687 
689 {
693  "BasisType is not a boundary interior form");
697  "BasisType is not a boundary interior form");
698 
699  return 2 * GetBasisNumModes(0) + 2 * GetBasisNumModes(1);
700 }
701 
703  const std::vector<unsigned int> &nummodes, int &modes_offset)
704 {
705  int nmodes = nummodes[modes_offset] * nummodes[modes_offset + 1];
706  modes_offset += 2;
707 
708  return nmodes;
709 }
710 
712 {
713  bool returnval = false;
714 
717  {
720  {
721  returnval = true;
722  }
723  }
724 
725  return returnval;
726 }
727 
729  Array<OneD, NekDouble> &coords_1,
730  Array<OneD, NekDouble> &coords_2)
731 {
732  boost::ignore_unused(coords_2);
733  Array<OneD, const NekDouble> z0 = m_base[0]->GetZ();
734  Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
735  int nq0 = GetNumPoints(0);
736  int nq1 = GetNumPoints(1);
737  int i;
738 
739  for (i = 0; i < nq1; ++i)
740  {
741  Blas::Dcopy(nq0, z0.get(), 1, &coords_0[0] + i * nq0, 1);
742  Vmath::Fill(nq0, z1[i], &coords_1[0] + i * nq0, 1);
743  }
744 }
745 
746 /**
747  * @brief This function evaluates the basis function mode @p mode at a
748  * point @p coords of the domain.
749  *
750  * This function uses barycentric interpolation with the tensor
751  * product separation of the basis function to improve performance.
752  *
753  * @param coord The coordinate inside the standard region.
754  * @param mode The mode number to be evaluated.
755  *
756  * @return The value of the basis function @p mode at @p coords.
757  */
759  const Array<OneD, const NekDouble> &coords, int mode)
760 {
761  ASSERTL2(coords[0] > -1 - NekConstants::kNekZeroTol, "coord[0] < -1");
762  ASSERTL2(coords[0] < 1 + NekConstants::kNekZeroTol, "coord[0] > 1");
763  ASSERTL2(coords[1] > -1 - NekConstants::kNekZeroTol, "coord[1] < -1");
764  ASSERTL2(coords[1] < 1 + NekConstants::kNekZeroTol, "coord[1] > 1");
765 
766  const int nm0 = m_base[0]->GetNumModes();
767  const int nm1 = m_base[1]->GetNumModes();
768 
769  return StdExpansion::BaryEvaluateBasis<0>(coords[0], mode % nm1) *
770  StdExpansion::BaryEvaluateBasis<1>(coords[1], mode / nm0);
771 }
772 
773 //////////////
774 // Mappings //
775 //////////////
776 
778 {
779  int i;
780  int cnt = 0;
781  int nummodes0, nummodes1;
782  int value1 = 0, value2 = 0;
783  if (outarray.size() != NumBndryCoeffs())
784  {
786  }
787 
788  nummodes0 = m_base[0]->GetNumModes();
789  nummodes1 = m_base[1]->GetNumModes();
790 
791  const LibUtilities::BasisType Btype0 = GetBasisType(0);
792  const LibUtilities::BasisType Btype1 = GetBasisType(1);
793 
794  switch (Btype1)
795  {
798  value1 = nummodes0;
799  break;
801  value1 = 2 * nummodes0;
802  break;
803  default:
804  ASSERTL0(0, "Mapping array is not defined for this expansion");
805  break;
806  }
807 
808  for (i = 0; i < value1; i++)
809  {
810  outarray[i] = i;
811  }
812  cnt = value1;
813 
814  switch (Btype0)
815  {
818  value2 = value1 + nummodes0 - 1;
819  break;
821  value2 = value1 + 1;
822  break;
823  default:
824  ASSERTL0(0, "Mapping array is not defined for this expansion");
825  break;
826  }
827 
828  for (i = 0; i < nummodes1 - 2; i++)
829  {
830  outarray[cnt++] = value1 + i * nummodes0;
831  outarray[cnt++] = value2 + i * nummodes0;
832  }
833 
834  if (Btype1 == LibUtilities::eGLL_Lagrange ||
836  {
837  for (i = nummodes0 * (nummodes1 - 1); i < GetNcoeffs(); i++)
838  {
839  outarray[cnt++] = i;
840  }
841  }
842 }
843 
845 {
846  int i, j;
847  int cnt = 0;
848  int nummodes0, nummodes1;
849  int startvalue = 0;
850  if (outarray.size() != GetNcoeffs() - NumBndryCoeffs())
851  {
853  }
854 
855  nummodes0 = m_base[0]->GetNumModes();
856  nummodes1 = m_base[1]->GetNumModes();
857 
858  const LibUtilities::BasisType Btype0 = GetBasisType(0);
859  const LibUtilities::BasisType Btype1 = GetBasisType(1);
860 
861  switch (Btype1)
862  {
864  startvalue = nummodes0;
865  break;
867  startvalue = 2 * nummodes0;
868  break;
869  default:
870  ASSERTL0(0, "Mapping array is not defined for this expansion");
871  break;
872  }
873 
874  switch (Btype0)
875  {
877  startvalue++;
878  break;
880  startvalue += 2;
881  break;
882  default:
883  ASSERTL0(0, "Mapping array is not defined for this expansion");
884  break;
885  }
886 
887  for (i = 0; i < nummodes1 - 2; i++)
888  {
889  for (j = 0; j < nummodes0 - 2; j++)
890  {
891  outarray[cnt++] = startvalue + j;
892  }
893  startvalue += nummodes0;
894  }
895 }
896 
897 int StdQuadExp::v_GetVertexMap(int localVertexId, bool useCoeffPacking)
898 {
899  int localDOF = 0;
900 
901  if (useCoeffPacking == true)
902  {
903  switch (localVertexId)
904  {
905  case 0:
906  {
907  localDOF = 0;
908  }
909  break;
910  case 1:
911  {
913  {
914  localDOF = m_base[0]->GetNumModes() - 1;
915  }
916  else
917  {
918  localDOF = 1;
919  }
920  }
921  break;
922  case 2:
923  {
925  {
926  localDOF = m_base[0]->GetNumModes() *
927  (m_base[1]->GetNumModes() - 1);
928  }
929  else
930  {
931  localDOF = m_base[0]->GetNumModes();
932  }
933  }
934  break;
935  case 3:
936  {
938  {
939  localDOF =
940  m_base[0]->GetNumModes() * m_base[1]->GetNumModes() - 1;
941  }
942  else
943  {
944  localDOF = m_base[0]->GetNumModes() + 1;
945  }
946  }
947  break;
948  default:
949  ASSERTL0(false, "eid must be between 0 and 3");
950  break;
951  }
952  }
953  else
954  {
955  switch (localVertexId)
956  {
957  case 0:
958  {
959  localDOF = 0;
960  }
961  break;
962  case 1:
963  {
965  {
966  localDOF = m_base[0]->GetNumModes() - 1;
967  }
968  else
969  {
970  localDOF = 1;
971  }
972  }
973  break;
974  case 2:
975  {
977  {
978  localDOF =
979  m_base[0]->GetNumModes() * m_base[1]->GetNumModes() - 1;
980  }
981  else
982  {
983  localDOF = m_base[0]->GetNumModes() + 1;
984  }
985  }
986  break;
987  case 3:
988  {
990  {
991  localDOF = m_base[0]->GetNumModes() *
992  (m_base[1]->GetNumModes() - 1);
993  }
994  else
995  {
996  localDOF = m_base[0]->GetNumModes();
997  }
998  }
999  break;
1000  default:
1001  ASSERTL0(false, "eid must be between 0 and 3");
1002  break;
1003  }
1004  }
1005  return localDOF;
1006 }
1007 
1008 /** \brief Get the map of the coefficient location to teh
1009  * local trace coefficients
1010  */
1011 
1012 void StdQuadExp::v_GetTraceCoeffMap(const unsigned int traceid,
1013  Array<OneD, unsigned int> &maparray)
1014 {
1015  ASSERTL1(traceid < 4, "traceid must be between 0 and 3");
1016 
1017  unsigned int i;
1018  unsigned int order0 = m_base[0]->GetNumModes();
1019  unsigned int order1 = m_base[1]->GetNumModes();
1020  unsigned int numModes = (traceid % 2) ? order1 : order0;
1021 
1022  if (maparray.size() != numModes)
1023  {
1024  maparray = Array<OneD, unsigned int>(numModes);
1025  }
1026 
1027  const LibUtilities::BasisType bType = GetBasisType(traceid % 2);
1028 
1029  if (bType == LibUtilities::eModified_A)
1030  {
1031  switch (traceid)
1032  {
1033  case 0:
1034  {
1035  for (i = 0; i < numModes; i++)
1036  {
1037  maparray[i] = i;
1038  }
1039  }
1040  break;
1041  case 1:
1042  {
1043  for (i = 0; i < numModes; i++)
1044  {
1045  maparray[i] = i * order0 + 1;
1046  }
1047  }
1048  break;
1049  case 2:
1050  {
1051  for (i = 0; i < numModes; i++)
1052  {
1053  maparray[i] = order0 + i;
1054  }
1055  }
1056  break;
1057  case 3:
1058  {
1059  for (i = 0; i < numModes; i++)
1060  {
1061  maparray[i] = i * order0;
1062  }
1063  }
1064  break;
1065  default:
1066  break;
1067  }
1068  }
1069  else if (bType == LibUtilities::eGLL_Lagrange ||
1071  {
1072  switch (traceid)
1073  {
1074  case 0:
1075  {
1076  for (i = 0; i < numModes; i++)
1077  {
1078  maparray[i] = i;
1079  }
1080  }
1081  break;
1082  case 1:
1083  {
1084  for (i = 0; i < numModes; i++)
1085  {
1086  maparray[i] = (i + 1) * order0 - 1;
1087  }
1088  }
1089  break;
1090  case 2:
1091  {
1092  for (i = 0; i < numModes; i++)
1093  {
1094  maparray[i] = order0 * (order1 - 1) + i;
1095  }
1096  }
1097  break;
1098  case 3:
1099  {
1100  for (i = 0; i < numModes; i++)
1101  {
1102  maparray[i] = order0 * i;
1103  }
1104  }
1105  break;
1106  default:
1107  break;
1108  }
1109  }
1110  else
1111  {
1112  ASSERTL0(false, "Mapping not defined for this type of basis");
1113  }
1114 }
1115 
1117  const int eid, Array<OneD, unsigned int> &maparray,
1118  Array<OneD, int> &signarray, const Orientation edgeOrient)
1119 {
1120  int i;
1121  const int nummodes0 = m_base[0]->GetNumModes();
1122  const int nummodes1 = m_base[1]->GetNumModes();
1123  const int nEdgeIntCoeffs = GetTraceNcoeffs(eid) - 2;
1124  const LibUtilities::BasisType bType = GetBasisType(eid % 2);
1125 
1126  if (maparray.size() != nEdgeIntCoeffs)
1127  {
1128  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1129  }
1130 
1131  if (signarray.size() != nEdgeIntCoeffs)
1132  {
1133  signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1134  }
1135  else
1136  {
1137  fill(signarray.get(), signarray.get() + nEdgeIntCoeffs, 1);
1138  }
1139 
1140  if (bType == LibUtilities::eModified_A)
1141  {
1142  switch (eid)
1143  {
1144  case 0:
1145  {
1146  for (i = 0; i < nEdgeIntCoeffs; i++)
1147  {
1148  maparray[i] = i + 2;
1149  }
1150  }
1151  break;
1152  case 1:
1153  {
1154  for (i = 0; i < nEdgeIntCoeffs; i++)
1155  {
1156  maparray[i] = (i + 2) * nummodes0 + 1;
1157  }
1158  }
1159  break;
1160  case 2:
1161  {
1162  for (i = 0; i < nEdgeIntCoeffs; i++)
1163  {
1164  maparray[i] = nummodes0 + i + 2;
1165  }
1166  }
1167  break;
1168  case 3:
1169  {
1170  for (i = 0; i < nEdgeIntCoeffs; i++)
1171  {
1172  maparray[i] = (i + 2) * nummodes0;
1173  }
1174  }
1175  break;
1176  default:
1177  ASSERTL0(false, "eid must be between 0 and 3");
1178  break;
1179  }
1180 
1181  if (edgeOrient == eBackwards)
1182  {
1183  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1184  {
1185  signarray[i] = -1;
1186  }
1187  }
1188  }
1189  else if (bType == LibUtilities::eGLL_Lagrange)
1190  {
1191  switch (eid)
1192  {
1193  case 0:
1194  {
1195  for (i = 0; i < nEdgeIntCoeffs; i++)
1196  {
1197  maparray[i] = i + 1;
1198  }
1199  }
1200  break;
1201  case 1:
1202  {
1203  for (i = 0; i < nEdgeIntCoeffs; i++)
1204  {
1205  maparray[i] = (i + 2) * nummodes0 - 1;
1206  }
1207  }
1208  break;
1209  case 2:
1210  {
1211  for (i = 0; i < nEdgeIntCoeffs; i++)
1212  {
1213  maparray[i] = nummodes0 * (nummodes1 - 1) + i + 1;
1214  }
1215  }
1216  break;
1217  case 3:
1218  {
1219  for (i = 0; i < nEdgeIntCoeffs; i++)
1220  {
1221  maparray[i] = nummodes0 * (i + 1);
1222  }
1223  }
1224  break;
1225  default:
1226  ASSERTL0(false, "eid must be between 0 and 3");
1227  break;
1228  }
1229  if (edgeOrient == eBackwards)
1230  {
1231  reverse(maparray.get(), maparray.get() + nEdgeIntCoeffs);
1232  }
1233  }
1234  else
1235  {
1236  ASSERTL0(false, "Mapping not defined for this type of basis");
1237  }
1238 }
1239 
1240 ///////////////////////
1241 // Wrapper Functions //
1242 ///////////////////////
1243 
1245 {
1246  int i;
1247  int order0 = GetBasisNumModes(0);
1248  int order1 = GetBasisNumModes(1);
1249  MatrixType mtype = mkey.GetMatrixType();
1250 
1251  DNekMatSharedPtr Mat;
1252 
1253  switch (mtype)
1254  {
1256  {
1257  int nq0 = m_base[0]->GetNumPoints();
1258  int nq1 = m_base[1]->GetNumPoints();
1259  int nq;
1260 
1261  // take definition from key
1262  if (mkey.ConstFactorExists(eFactorConst))
1263  {
1264  nq = (int)mkey.GetConstFactor(eFactorConst);
1265  }
1266  else
1267  {
1268  nq = max(nq0, nq1);
1269  }
1270 
1271  int neq =
1274  Array<OneD, NekDouble> coll(2);
1276  Array<OneD, NekDouble> tmp(nq0);
1277 
1278  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1);
1279  int cnt = 0;
1280 
1281  for (int i = 0; i < nq; ++i)
1282  {
1283  for (int j = 0; j < nq; ++j, ++cnt)
1284  {
1285  coords[cnt] = Array<OneD, NekDouble>(2);
1286  coords[cnt][0] = -1.0 + 2 * j / (NekDouble)(nq - 1);
1287  coords[cnt][1] = -1.0 + 2 * i / (NekDouble)(nq - 1);
1288  }
1289  }
1290 
1291  for (int i = 0; i < neq; ++i)
1292  {
1293  LocCoordToLocCollapsed(coords[i], coll);
1294 
1295  I[0] = m_base[0]->GetI(coll);
1296  I[1] = m_base[1]->GetI(coll + 1);
1297 
1298  // interpolate first coordinate direction
1299  for (int j = 0; j < nq1; ++j)
1300  {
1301  NekDouble fac = (I[1]->GetPtr())[j];
1302  Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1303 
1304  Vmath::Vcopy(nq0, &tmp[0], 1,
1305  Mat->GetRawPtr() + j * nq0 * neq + i, neq);
1306  }
1307  }
1308  break;
1309  }
1310  case eMass:
1311  {
1313  // For Fourier basis set the imaginary component of mean mode
1314  // to have a unit diagonal component in mass matrix
1316  {
1317  for (i = 0; i < order1; ++i)
1318  {
1319  (*Mat)(order0 * i + 1, i * order0 + 1) = 1.0;
1320  }
1321  }
1322 
1324  {
1325  for (i = 0; i < order0; ++i)
1326  {
1327  (*Mat)(order0 + i, order0 + i) = 1.0;
1328  }
1329  }
1330  break;
1331  }
1332  case eFwdTrans:
1333  {
1334  Mat =
1336  StdMatrixKey iprodkey(eIProductWRTBase, DetShapeType(), *this);
1337  DNekMat &Iprod = *GetStdMatrix(iprodkey);
1338  StdMatrixKey imasskey(eInvMass, DetShapeType(), *this);
1339  DNekMat &Imass = *GetStdMatrix(imasskey);
1340 
1341  (*Mat) = Imass * Iprod;
1342  break;
1343  }
1344  case eGaussDG:
1345  {
1346  ConstFactorMap factors = mkey.GetConstFactors();
1347 
1348  int edge = (int)factors[StdRegions::eFactorGaussEdge];
1349  int dir = (edge + 1) % 2;
1350  int nCoeffs = m_base[dir]->GetNumModes();
1351 
1352  const LibUtilities::PointsKey BS_p(
1355  nCoeffs, BS_p);
1356 
1357  Array<OneD, NekDouble> coords(1, 0.0);
1358  coords[0] = (edge == 0 || edge == 3) ? -1.0 : 1.0;
1359 
1362  DNekMatSharedPtr m_Ix = basis->GetI(coords);
1363 
1364  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(1.0, nCoeffs);
1365  Vmath::Vcopy(nCoeffs, m_Ix->GetPtr(), 1, Mat->GetPtr(), 1);
1366  break;
1367  }
1368  default:
1369  {
1371  break;
1372  }
1373  }
1374 
1375  return Mat;
1376 }
1377 
1379 {
1380  return GenMatrix(mkey);
1381 }
1382 
1383 ///////////////////////////////////
1384 // Operator evaluation functions //
1385 ///////////////////////////////////
1386 
1388  const Array<OneD, const NekDouble> &inarray,
1389  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1390 {
1391 
1393 
1394  if (inarray.get() == outarray.get())
1395  {
1397  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, tmp.get(), 1);
1398 
1399  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1400  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1401  }
1402  else
1403  {
1404  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1405  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1406  }
1407 }
1408 
1410  const StdMatrixKey &mkey)
1411 {
1412  // Generate an orthonogal expansion
1413  int qa = m_base[0]->GetNumPoints();
1414  int qb = m_base[1]->GetNumPoints();
1415  int nmodes_a = m_base[0]->GetNumModes();
1416  int nmodes_b = m_base[1]->GetNumModes();
1417  int nmodes = min(nmodes_a, nmodes_b);
1418  // Declare orthogonal basis.
1421 
1424  StdQuadExp OrthoExp(Ba, Bb);
1425 
1426  // SVV parameters loaded from the .xml case file
1427  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1428 
1429  // project onto modal space.
1430  OrthoExp.FwdTrans(array, orthocoeffs);
1431 
1432  if (mkey.ConstFactorExists(
1433  eFactorSVVPowerKerDiffCoeff)) // Rodrigo's power kernel
1434  {
1436  NekDouble SvvDiffCoeff =
1439 
1440  for (int j = 0; j < nmodes_a; ++j)
1441  {
1442  for (int k = 0; k < nmodes_b; ++k)
1443  {
1444  // linear space but makes high modes very negative
1445  NekDouble fac = std::max(
1446  pow((1.0 * j) / (nmodes_a - 1), cutoff * nmodes_a),
1447  pow((1.0 * k) / (nmodes_b - 1), cutoff * nmodes_b));
1448 
1449  orthocoeffs[j * nmodes_b + k] *= SvvDiffCoeff * fac;
1450  }
1451  }
1452  }
1453  else if (mkey.ConstFactorExists(
1454  eFactorSVVDGKerDiffCoeff)) // Rodrigo/mansoor's DG kernel
1455  {
1456  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff) *
1458  int max_ab = max(nmodes_a - kSVVDGFiltermodesmin,
1459  nmodes_b - kSVVDGFiltermodesmin);
1460  max_ab = max(max_ab, 0);
1461  max_ab = min(max_ab, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
1462 
1463  for (int j = 0; j < nmodes_a; ++j)
1464  {
1465  for (int k = 0; k < nmodes_b; ++k)
1466  {
1467  int maxjk = max(j, k);
1468  maxjk = min(maxjk, kSVVDGFiltermodesmax - 1);
1469 
1470  orthocoeffs[j * nmodes_b + k] *=
1471  SvvDiffCoeff * kSVVDGFilter[max_ab][maxjk];
1472  }
1473  }
1474  }
1475  else
1476  {
1477  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
1478  // Exponential Kernel implementation
1479  int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio) *
1480  min(nmodes_a, nmodes_b));
1481 
1482  // counters for scanning through orthocoeffs array
1483  int cnt = 0;
1484 
1485  //------"New" Version August 22nd '13--------------------
1486  for (int j = 0; j < nmodes_a; ++j)
1487  {
1488  for (int k = 0; k < nmodes_b; ++k)
1489  {
1490  if (j + k >= cutoff) // to filter out only the "high-modes"
1491  {
1492  orthocoeffs[j * nmodes_b + k] *=
1493  (SvvDiffCoeff *
1494  exp(-(j + k - nmodes) * (j + k - nmodes) /
1495  ((NekDouble)((j + k - cutoff + 1) *
1496  (j + k - cutoff + 1)))));
1497  }
1498  else
1499  {
1500  orthocoeffs[j * nmodes_b + k] *= 0.0;
1501  }
1502  cnt++;
1503  }
1504  }
1505  }
1506 
1507  // backward transform to physical space
1508  OrthoExp.BwdTrans(orthocoeffs, array);
1509 }
1510 
1512  const NekDouble alpha,
1513  const NekDouble exponent,
1514  const NekDouble cutoff)
1515 {
1516  // Generate an orthogonal expansion
1517  int qa = m_base[0]->GetNumPoints();
1518  int qb = m_base[1]->GetNumPoints();
1519  int nmodesA = m_base[0]->GetNumModes();
1520  int nmodesB = m_base[1]->GetNumModes();
1521  int P = nmodesA - 1;
1522  int Q = nmodesB - 1;
1523 
1524  // Declare orthogonal basis.
1527 
1530  StdQuadExp OrthoExp(Ba, Bb);
1531 
1532  // Cutoff
1533  int Pcut = cutoff * P;
1534  int Qcut = cutoff * Q;
1535 
1536  // Project onto orthogonal space.
1537  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1538  OrthoExp.FwdTrans(array, orthocoeffs);
1539 
1540  //
1541  NekDouble fac, fac1, fac2;
1542  for (int i = 0; i < nmodesA; ++i)
1543  {
1544  for (int j = 0; j < nmodesB; ++j)
1545  {
1546  // to filter out only the "high-modes"
1547  if (i > Pcut || j > Qcut)
1548  {
1549  fac1 = (NekDouble)(i - Pcut) / ((NekDouble)(P - Pcut));
1550  fac2 = (NekDouble)(j - Qcut) / ((NekDouble)(Q - Qcut));
1551  fac = max(fac1, fac2);
1552  fac = pow(fac, exponent);
1553  orthocoeffs[i * nmodesB + j] *= exp(-alpha * fac);
1554  }
1555  }
1556  }
1557 
1558  // backward transform to physical space
1559  OrthoExp.BwdTrans(orthocoeffs, array);
1560 }
1561 
1563  int numMin, const Array<OneD, const NekDouble> &inarray,
1564  Array<OneD, NekDouble> &outarray)
1565 {
1566  int n_coeffs = inarray.size();
1567 
1568  Array<OneD, NekDouble> coeff(n_coeffs);
1569  Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
1572 
1573  int nmodes0 = m_base[0]->GetNumModes();
1574  int nmodes1 = m_base[1]->GetNumModes();
1575  int numMax = nmodes0;
1576 
1577  Vmath::Vcopy(n_coeffs, inarray, 1, coeff_tmp, 1);
1578 
1579  const LibUtilities::PointsKey Pkey0(nmodes0,
1581  const LibUtilities::PointsKey Pkey1(nmodes1,
1583 
1584  LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1585  LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1586 
1587  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1588  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_A, nmodes1, Pkey1);
1589 
1590  LibUtilities::InterpCoeff2D(b0, b1, coeff_tmp, bortho0, bortho1, coeff);
1591 
1592  Vmath::Zero(n_coeffs, coeff_tmp, 1);
1593 
1594  int cnt = 0;
1595  for (int i = 0; i < numMin + 1; ++i)
1596  {
1597  Vmath::Vcopy(numMin, tmp = coeff + cnt, 1, tmp2 = coeff_tmp + cnt, 1);
1598 
1599  cnt = i * numMax;
1600  }
1601 
1602  LibUtilities::InterpCoeff2D(bortho0, bortho1, coeff_tmp, b0, b1, outarray);
1603 }
1604 
1606  Array<OneD, NekDouble> &outarray,
1607  const StdMatrixKey &mkey)
1608 {
1609  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1610 }
1611 
1613  const Array<OneD, const NekDouble> &inarray,
1614  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1615 {
1616  StdQuadExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1617 }
1618 
1620  const int k1, const int k2, const Array<OneD, const NekDouble> &inarray,
1621  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1622 {
1623  StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1624 }
1625 
1627  const int i, const Array<OneD, const NekDouble> &inarray,
1628  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1629 {
1630  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1631 }
1632 
1634  const Array<OneD, const NekDouble> &inarray,
1635  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1636 {
1637  StdQuadExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1638 }
1639 
1640 // up to here
1642  const Array<OneD, const NekDouble> &inarray,
1643  Array<OneD, NekDouble> &outarray)
1644 {
1645  int i;
1646  int nquad0 = m_base[0]->GetNumPoints();
1647  int nquad1 = m_base[1]->GetNumPoints();
1648 
1649  const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
1650  const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
1651 
1652  // multiply by integration constants
1653  for (i = 0; i < nquad1; ++i)
1654  {
1655  Vmath::Vmul(nquad0, inarray.get() + i * nquad0, 1, w0.get(), 1,
1656  outarray.get() + i * nquad0, 1);
1657  }
1658 
1659  for (i = 0; i < nquad0; ++i)
1660  {
1661  Vmath::Vmul(nquad1, outarray.get() + i, nquad0, w1.get(), 1,
1662  outarray.get() + i, nquad0);
1663  }
1664 }
1665 
1667  bool standard)
1668 {
1669  boost::ignore_unused(standard);
1670 
1671  int np1 = m_base[0]->GetNumPoints();
1672  int np2 = m_base[1]->GetNumPoints();
1673  int np = max(np1, np2);
1674 
1675  conn = Array<OneD, int>(6 * (np - 1) * (np - 1));
1676 
1677  int row = 0;
1678  int rowp1 = 0;
1679  int cnt = 0;
1680  for (int i = 0; i < np - 1; ++i)
1681  {
1682  rowp1 += np;
1683  for (int j = 0; j < np - 1; ++j)
1684  {
1685  conn[cnt++] = row + j;
1686  conn[cnt++] = row + j + 1;
1687  conn[cnt++] = rowp1 + j;
1688 
1689  conn[cnt++] = rowp1 + j + 1;
1690  conn[cnt++] = rowp1 + j;
1691  conn[cnt++] = row + j + 1;
1692  }
1693  row += np;
1694  }
1695 }
1696 
1697 } // namespace StdRegions
1698 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:249
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:272
#define sign(a, b)
return the sign(b)*a
Definition: Polylib.cpp:15
Describes the specification for a Basis.
Definition: Basis.h:50
Defines a specification for a set of points.
Definition: Points.h:59
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void PhysTensorDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray_d0, Array< OneD, NekDouble > &outarray_d1)
Calculate the 2D derivative in the local tensor/collapsed coordinate at the physical points.
NekDouble Integral(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &w0, const Array< OneD, const NekDouble > &w1)
void BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
virtual void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
virtual void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0=true, bool doCheckCollDir1=true)
The base class for all shapes.
Definition: StdExpansion.h:71
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:130
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
const LibUtilities::BasisSharedPtr & GetBasis(int dir) const
This function gets the shared point to basis in the dir direction.
Definition: StdExpansion.h:118
void WeakDerivMatrixOp_MatFree(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:163
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:611
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
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:760
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
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
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void GetTraceToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation traceOrient=eForwards, int P=-1, int Q=-1)
Definition: StdExpansion.h:692
LibUtilities::ShapeType DetShapeType() const
This function returns the shape of the expansion domain.
Definition: StdExpansion.h:375
void GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdExpansion.h:682
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
int GetTraceNcoeffs(const int i) const
This function returns the number of expansion coefficients belonging to the i-th trace.
Definition: StdExpansion.h:269
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:845
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:213
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
LibUtilities::NekManager< StdMatrixKey, DNekBlkMat, StdMatrixKey::opLess > m_stdStaticCondMatrixManager
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:226
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:731
LibUtilities::NekManager< StdMatrixKey, DNekMat, StdMatrixKey::opLess > m_stdMatrixManager
int GetBasisNumModes(const int dir) const
This function returns the number of expansion modes in the dir direction.
Definition: StdExpansion.h:176
Array< OneD, LibUtilities::BasisSharedPtr > m_base
void MassMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:140
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:126
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:135
virtual int v_NumBndryCoeffs() const
Definition: StdQuadExp.cpp:674
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:430
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Transform a given function from physical quadrature space to coefficient space.
Definition: StdQuadExp.cpp:227
virtual const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int j) const
Definition: StdQuadExp.cpp:653
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:164
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:418
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=NullNekDouble1DArray)
Definition: StdQuadExp.cpp:128
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Calculate the inner product of inarray with respect to the basis B=base0*base1 and put into outarray.
Definition: StdQuadExp.cpp:377
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdQuadExp.cpp:777
void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual int v_NumDGBndryCoeffs() const
Definition: StdQuadExp.cpp:688
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdQuadExp.cpp:390
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Definition: StdQuadExp.cpp:77
virtual void v_ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
virtual void v_GetTraceCoeffMap(const unsigned int traceid, Array< OneD, unsigned int > &maparray)
Get the map of the coefficient location to teh local trace coefficients.
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdQuadExp.cpp:669
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdQuadExp.cpp:702
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
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=NullNekDouble1DArray)
Calculate the derivative of the physical points.
Definition: StdQuadExp.cpp:95
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:250
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdQuadExp.cpp:711
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdQuadExp.cpp:561
virtual int v_GetNverts() const
Definition: StdQuadExp.cpp:615
virtual int v_GetNtraces() const
Definition: StdQuadExp.cpp:620
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &array)
Fill outarray with mode mode of expansion.
Definition: StdQuadExp.cpp:581
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdQuadExp.cpp:897
virtual NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode)
This function evaluates the basis function mode mode at a point coords of the domain.
Definition: StdQuadExp.cpp:758
virtual void v_IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1)
Definition: StdQuadExp.cpp:503
virtual void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
Definition: StdQuadExp.cpp:568
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdQuadExp.cpp:844
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_0, Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2)
Definition: StdQuadExp.cpp:728
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:468
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:437
virtual int v_GetTraceNcoeffs(const int i) const
Definition: StdQuadExp.cpp:625
virtual int v_GetTraceNumPoints(const int i) const
Definition: StdQuadExp.cpp:639
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
virtual void v_BwdTrans_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1)
Definition: StdQuadExp.cpp:184
virtual void v_GetTraceInteriorToElementMap(const int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation edgeOrient=eForwards)
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdQuadExp.cpp:150
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
static void Dcopy(const int &n, const double *x, const int &incx, double *y, const int &incy)
BLAS level 1: Copy x to y.
Definition: Blas.hpp:147
static double Ddot(const int &n, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: output = .
Definition: Blas.hpp:182
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
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:136
BasisManagerT & BasisManager(void)
std::shared_ptr< Basis > BasisSharedPtr
void InterpCoeff2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
Definition: InterpCoeff.cpp:71
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:59
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:58
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
@ eFourier
Fourier Expansion .
Definition: BasisType.h:57
static const NekDouble kNekZeroTol
const int kSVVDGFiltermodesmin
Definition: StdRegions.hpp:352
const int kSVVDGFiltermodesmax
Definition: StdRegions.hpp:353
const NekDouble kSVVDGFilter[9][11]
Definition: StdRegions.hpp:355
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
std::shared_ptr< StdSegExp > StdSegExpSharedPtr
Definition: StdSegExp.h:46
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
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
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
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.cpp:492
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1255
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
Definition: Vmath.cpp:419