Nektar++
StdTriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File StdTriExp.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: Triangle routines built upon StdExpansion2D
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
39 #include <StdRegions/StdSegExp.h> // for StdSegExp, etc
40 #include <StdRegions/StdTriExp.h>
41 
42 using namespace std;
43 
44 namespace Nektar
45 {
46 namespace StdRegions
47 {
48 
49 StdTriExp::StdTriExp()
50 {
51 }
52 
53 StdTriExp::StdTriExp(const LibUtilities::BasisKey &Ba,
54  const LibUtilities::BasisKey &Bb)
55  : StdExpansion(LibUtilities::StdTriData::getNumberOfCoefficients(
56  Ba.GetNumModes(), Bb.GetNumModes()),
57  2, Ba, Bb),
58  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
59  Ba.GetNumModes(), Bb.GetNumModes()),
60  Ba, Bb)
61 {
62  ASSERTL0(Ba.GetNumModes() <= Bb.GetNumModes(),
63  "order in 'a' direction is higher than order "
64  "in 'b' direction");
65 }
66 
68 {
69 }
70 
72 {
73 }
74 
75 //-------------------------------
76 // Integration Methods
77 //-------------------------------
79 {
80  int i;
81  int nquad1 = m_base[1]->GetNumPoints();
82  Array<OneD, NekDouble> w1_tmp(nquad1);
83 
84  Array<OneD, const NekDouble> w0 = m_base[0]->GetW();
85  Array<OneD, const NekDouble> w1 = m_base[1]->GetW();
86  Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
87 
88  switch (m_base[1]->GetPointsType())
89  {
90  case LibUtilities::eGaussRadauMAlpha1Beta0: // (0,1) Jacobi Inner
91  // product
92  {
93  Vmath::Smul(nquad1, 0.5, w1, 1, w1_tmp, 1);
94  break;
95  }
96  default:
97  {
98  // include jacobian factor on whatever coordinates are defined.
99  for (i = 0; i < nquad1; ++i)
100  {
101  w1_tmp[i] = 0.5 * (1 - z1[i]) * w1[i];
102  }
103  break;
104  }
105  }
106 
107  return StdExpansion2D::Integral(inarray, w0, w1_tmp);
108 }
109 
110 //-----------------------------
111 // Differentiation Methods
112 //-----------------------------
113 
114 /**
115  * \brief Calculate the derivative of the physical points.
116  *
117  * \f$ \frac{\partial u}{\partial x_1} = \left .
118  * \frac{2.0}{1-\eta_2} \frac{\partial u}{\partial d\eta_1}
119  * \right |_{\eta_2}\f$
120  *
121  * \f$ \frac{\partial u}{\partial x_2} = \left .
122  * \frac{1+\eta_1}{1-\eta_2} \frac{\partial u}{\partial d\eta_1}
123  * \right |_{\eta_2} + \left . \frac{\partial u}{\partial d\eta_2}
124  * \right |_{\eta_1} \f$
125  */
127  Array<OneD, NekDouble> &out_d0,
128  Array<OneD, NekDouble> &out_d1,
129  Array<OneD, NekDouble> &out_d2)
130 {
131  boost::ignore_unused(out_d2);
132 
133  int i;
134  int nquad0 = m_base[0]->GetNumPoints();
135  int nquad1 = m_base[1]->GetNumPoints();
136  Array<OneD, NekDouble> wsp(std::max(nquad0, nquad1));
137 
138  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
139  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
140 
141  // set up geometric factor: 2/(1-z1)
142  Vmath::Sadd(nquad1, -1.0, z1, 1, wsp, 1);
143  Vmath::Sdiv(nquad1, -2.0, wsp, 1, wsp, 1);
144 
145  if (out_d0.size() > 0)
146  {
147  PhysTensorDeriv(inarray, out_d0, out_d1);
148 
149  for (i = 0; i < nquad1; ++i)
150  {
151  Blas::Dscal(nquad0, wsp[i], &out_d0[0] + i * nquad0, 1);
152  }
153 
154  // if no d1 required do not need to calculate both deriv
155  if (out_d1.size() > 0)
156  {
157  // set up geometric factor: (1+z0)/2
158  Vmath::Sadd(nquad0, 1.0, z0, 1, wsp, 1);
159  Vmath::Smul(nquad0, 0.5, wsp, 1, wsp, 1);
160 
161  for (i = 0; i < nquad1; ++i)
162  {
163  Vmath::Vvtvp(nquad0, &wsp[0], 1, &out_d0[0] + i * nquad0, 1,
164  &out_d1[0] + i * nquad0, 1,
165  &out_d1[0] + i * nquad0, 1);
166  }
167  }
168  }
169  else if (out_d1.size() > 0)
170  {
171  Array<OneD, NekDouble> diff0(nquad0 * nquad1);
172  PhysTensorDeriv(inarray, diff0, out_d1);
173 
174  for (i = 0; i < nquad1; ++i)
175  {
176  Blas::Dscal(nquad0, wsp[i], &diff0[0] + i * nquad0, 1);
177  }
178 
179  Vmath::Sadd(nquad0, 1.0, z0, 1, wsp, 1);
180  Vmath::Smul(nquad0, 0.5, wsp, 1, wsp, 1);
181 
182  for (i = 0; i < nquad1; ++i)
183  {
184  Vmath::Vvtvp(nquad0, &wsp[0], 1, &diff0[0] + i * nquad0, 1,
185  &out_d1[0] + i * nquad0, 1, &out_d1[0] + i * nquad0,
186  1);
187  }
188  }
189 }
190 
191 void StdTriExp::v_PhysDeriv(const int dir,
192  const Array<OneD, const NekDouble> &inarray,
193  Array<OneD, NekDouble> &outarray)
194 {
195  switch (dir)
196  {
197  case 0:
198  {
199  v_PhysDeriv(inarray, outarray, NullNekDouble1DArray);
200  break;
201  }
202  case 1:
203  {
204  v_PhysDeriv(inarray, NullNekDouble1DArray, outarray);
205  break;
206  }
207  default:
208  {
209  ASSERTL1(false, "input dir is out of range");
210  break;
211  }
212  }
213 }
214 
216  Array<OneD, NekDouble> &out_d0,
217  Array<OneD, NekDouble> &out_d1,
218  Array<OneD, NekDouble> &out_d2)
219 {
220  boost::ignore_unused(out_d2);
221  StdTriExp::v_PhysDeriv(inarray, out_d0, out_d1);
222 }
223 
224 void StdTriExp::v_StdPhysDeriv(const int dir,
225  const Array<OneD, const NekDouble> &inarray,
226  Array<OneD, NekDouble> &outarray)
227 {
228  StdTriExp::v_PhysDeriv(dir, inarray, outarray);
229 }
230 
231 //---------------------------------------
232 // Transforms
233 //---------------------------------------
234 
235 /**
236  * \brief Backward tranform for triangular elements
237  *
238  * @note 'q' (base[1]) runs fastest in this element.
239  */
241  Array<OneD, NekDouble> &outarray)
242 {
243  v_BwdTrans_SumFac(inarray, outarray);
244 }
245 
247  Array<OneD, NekDouble> &outarray)
248 {
250  m_base[1]->GetNumModes());
251 
252  BwdTrans_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(), inarray,
253  outarray, wsp);
254 }
255 
257  const Array<OneD, const NekDouble> &base0,
258  const Array<OneD, const NekDouble> &base1,
259  const Array<OneD, const NekDouble> &inarray,
261  bool doCheckCollDir0, bool doCheckCollDir1)
262 {
263  boost::ignore_unused(doCheckCollDir0, doCheckCollDir1);
264 
265  int i;
266  int mode;
267  int nquad0 = m_base[0]->GetNumPoints();
268  int nquad1 = m_base[1]->GetNumPoints();
269  int nmodes0 = m_base[0]->GetNumModes();
270  int nmodes1 = m_base[1]->GetNumModes();
271 
272  ASSERTL1(wsp.size() >= nquad0 * nmodes1,
273  "Workspace size is not sufficient");
276  "Basis[1] is not of general tensor type");
277 
278  for (i = mode = 0; i < nmodes0; ++i)
279  {
280  Blas::Dgemv('N', nquad1, nmodes1 - i, 1.0, base1.get() + mode * nquad1,
281  nquad1, &inarray[0] + mode, 1, 0.0, &wsp[0] + i * nquad1,
282  1);
283  mode += nmodes1 - i;
284  }
285 
286  // fix for modified basis by splitting top vertex mode
288  {
289  Blas::Daxpy(nquad1, inarray[1], base1.get() + nquad1, 1,
290  &wsp[0] + nquad1, 1);
291  }
292 
293  Blas::Dgemm('N', 'T', nquad0, nquad1, nmodes0, 1.0, base0.get(), nquad0,
294  &wsp[0], nquad1, 0.0, &outarray[0], nquad0);
295 }
296 
298  Array<OneD, NekDouble> &outarray)
299 {
300  v_IProductWRTBase(inarray, outarray);
301 
302  // get Mass matrix inverse
303  StdMatrixKey masskey(eInvMass, DetShapeType(), *this);
304  DNekMatSharedPtr matsys = GetStdMatrix(masskey);
305 
306  // copy inarray in case inarray == outarray
307  NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
308  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
309 
310  out = (*matsys) * in;
311 }
312 
314  const Array<OneD, const NekDouble> &inarray,
315  Array<OneD, NekDouble> &outarray)
316 {
317  int i, j;
318  int npoints[2] = {m_base[0]->GetNumPoints(), m_base[1]->GetNumPoints()};
319  int nmodes[2] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes()};
320 
321  fill(outarray.get(), outarray.get() + m_ncoeffs, 0.0);
322 
323  Array<OneD, NekDouble> physEdge[3];
324  Array<OneD, NekDouble> coeffEdge[3];
325  for (i = 0; i < 3; i++)
326  {
327  physEdge[i] = Array<OneD, NekDouble>(npoints[i != 0]);
328  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i != 0]);
329  }
330 
331  for (i = 0; i < npoints[0]; i++)
332  {
333  physEdge[0][i] = inarray[i];
334  }
335 
336  for (i = 0; i < npoints[1]; i++)
337  {
338  physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
339  physEdge[2][i] =
340  inarray[(npoints[1] - 1) * npoints[0] - i * npoints[0]];
341  }
342 
343  StdSegExpSharedPtr segexp[2] = {
345  m_base[0]->GetBasisKey()),
347  m_base[1]->GetBasisKey())};
348 
349  Array<OneD, unsigned int> mapArray;
350  Array<OneD, int> signArray;
351  NekDouble sign;
352 
353  for (i = 0; i < 3; i++)
354  {
355  segexp[i != 0]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
356 
357  GetTraceToElementMap(i, mapArray, signArray);
358  for (j = 0; j < nmodes[i != 0]; j++)
359  {
360  sign = (NekDouble)signArray[j];
361  outarray[mapArray[j]] = sign * coeffEdge[i][j];
362  }
363  }
364 
367 
368  StdMatrixKey masskey(eMass, DetShapeType(), *this);
369  MassMatrixOp(outarray, tmp0, masskey);
370  v_IProductWRTBase(inarray, tmp1);
371 
372  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
373 
374  // get Mass matrix inverse (only of interior DOF)
375  // use block (1,1) of the static condensed system
376  // note: this block alreay contains the inverse matrix
377  DNekMatSharedPtr matsys =
378  (m_stdStaticCondMatrixManager[masskey])->GetBlock(1, 1);
379 
380  int nBoundaryDofs = v_NumBndryCoeffs();
381  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
382 
383  Array<OneD, NekDouble> rhs(nInteriorDofs);
384  Array<OneD, NekDouble> result(nInteriorDofs);
385 
386  v_GetInteriorMap(mapArray);
387 
388  for (i = 0; i < nInteriorDofs; i++)
389  {
390  rhs[i] = tmp1[mapArray[i]];
391  }
392 
393  Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, 1.0, &(matsys->GetPtr())[0],
394  nInteriorDofs, rhs.get(), 1, 0.0, result.get(), 1);
395 
396  for (i = 0; i < nInteriorDofs; i++)
397  {
398  outarray[mapArray[i]] = result[i];
399  }
400 }
401 
402 //---------------------------------------
403 // Inner product functions
404 //---------------------------------------
405 
406 /**
407  * \brief Calculate the inner product of inarray with respect to the
408  * basis B=base0[p]*base1[pq] and put into outarray.
409  *
410  * \f$
411  * \begin{array}{rcl}
412  * I_{pq} = (\phi^A_q \phi^B_{pq}, u) &=&
413  * \sum_{i=0}^{nq_0}\sum_{j=0}^{nq_1}
414  * \phi^A_p(\eta_{0,i})\phi^B_{pq}(\eta_{1,j}) w^0_i w^1_j
415  * u(\xi_{0,i} \xi_{1,j}) \\
416  * & = & \sum_{i=0}^{nq_0} \phi^A_p(\eta_{0,i})
417  * \sum_{j=0}^{nq_1} \phi^B_{pq}(\eta_{1,j}) \tilde{u}_{i,j}
418  * \end{array}
419  * \f$
420  *
421  * where
422  *
423  * \f$ \tilde{u}_{i,j} = w^0_i w^1_j u(\xi_{0,i},\xi_{1,j}) \f$
424  *
425  * which can be implemented as
426  *
427  * \f$ f_{pj} = \sum_{i=0}^{nq_0} \phi^A_p(\eta_{0,i})
428  * \tilde{u}_{i,j}
429  * \rightarrow {\bf B_1 U} \f$
430  * \f$ I_{pq} = \sum_{j=0}^{nq_1} \phi^B_{pq}(\eta_{1,j}) f_{pj}
431  * \rightarrow {\bf B_2[p*skip] f[skip]} \f$
432  *
433  * \b Recall: \f$ \eta_{1} = \frac{2(1+\xi_1)}{(1-\xi_2)}-1, \,
434  * \eta_2 = \xi_2\f$
435  *
436  * \b Note: For the orthgonality of this expansion to be realised the
437  * 'q' ordering must run fastest in contrast to the Quad and Hex
438  * ordering where 'p' index runs fastest to be consistent with the
439  * quadrature ordering.
440  *
441  * In the triangular space the i (i.e. \f$\eta_1\f$ direction)
442  * ordering still runs fastest by convention.
443  */
445  Array<OneD, NekDouble> &outarray)
446 {
447  StdTriExp::v_IProductWRTBase_SumFac(inarray, outarray);
448 }
449 
451  const Array<OneD, const NekDouble> &inarray,
452  Array<OneD, NekDouble> &outarray)
453 {
454  int nq = GetTotPoints();
455  StdMatrixKey iprodmatkey(eIProductWRTBase, DetShapeType(), *this);
456  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
457 
458  Blas::Dgemv('N', m_ncoeffs, nq, 1.0, iprodmat->GetPtr().get(), m_ncoeffs,
459  inarray.get(), 1, 0.0, outarray.get(), 1);
460 }
461 
463  const Array<OneD, const NekDouble> &inarray,
464  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
465 {
466  int nquad0 = m_base[0]->GetNumPoints();
467  int nquad1 = m_base[1]->GetNumPoints();
468  int order0 = m_base[0]->GetNumModes();
469 
470  if (multiplybyweights)
471  {
472  Array<OneD, NekDouble> tmp(nquad0 * nquad1 + nquad1 * order0);
473  Array<OneD, NekDouble> wsp(tmp + nquad0 * nquad1);
474 
475  // multiply by integration constants
476  MultiplyByQuadratureMetric(inarray, tmp);
477  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
478  m_base[1]->GetBdata(), tmp, outarray, wsp);
479  }
480  else
481  {
482  Array<OneD, NekDouble> wsp(nquad1 * order0);
483  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
484  m_base[1]->GetBdata(), inarray, outarray,
485  wsp);
486  }
487 }
488 
490  const Array<OneD, const NekDouble> &base0,
491  const Array<OneD, const NekDouble> &base1,
492  const Array<OneD, const NekDouble> &inarray,
494  bool doCheckCollDir0, bool doCheckCollDir1)
495 {
496  boost::ignore_unused(doCheckCollDir0, doCheckCollDir1);
497 
498  int i;
499  int mode;
500  int nquad0 = m_base[0]->GetNumPoints();
501  int nquad1 = m_base[1]->GetNumPoints();
502  int nmodes0 = m_base[0]->GetNumModes();
503  int nmodes1 = m_base[1]->GetNumModes();
504 
505  ASSERTL1(wsp.size() >= nquad1 * nmodes0,
506  "Workspace size is not sufficient");
507 
508  Blas::Dgemm('T', 'N', nquad1, nmodes0, nquad0, 1.0, inarray.get(), nquad0,
509  base0.get(), nquad0, 0.0, wsp.get(), nquad1);
510 
511  // Inner product with respect to 'b' direction
512  for (mode = i = 0; i < nmodes0; ++i)
513  {
514  Blas::Dgemv('T', nquad1, nmodes1 - i, 1.0, base1.get() + mode * nquad1,
515  nquad1, wsp.get() + i * nquad1, 1, 0.0,
516  outarray.get() + mode, 1);
517  mode += nmodes1 - i;
518  }
519 
520  // fix for modified basis by splitting top vertex mode
522  {
523  outarray[1] +=
524  Blas::Ddot(nquad1, base1.get() + nquad1, 1, wsp.get() + nquad1, 1);
525  }
526 }
527 
529  const int dir, const Array<OneD, const NekDouble> &inarray,
530  Array<OneD, NekDouble> &outarray)
531 {
532  StdTriExp::v_IProductWRTDerivBase_SumFac(dir, inarray, outarray);
533 }
534 
536  const int dir, const Array<OneD, const NekDouble> &inarray,
537  Array<OneD, NekDouble> &outarray)
538 {
539  int nq = GetTotPoints();
541 
542  switch (dir)
543  {
544  case 0:
545  {
546  mtype = eIProductWRTDerivBase0;
547  break;
548  }
549  case 1:
550  {
551  mtype = eIProductWRTDerivBase1;
552  break;
553  }
554  default:
555  {
556  ASSERTL1(false, "input dir is out of range");
557  break;
558  }
559  }
560 
561  StdMatrixKey iprodmatkey(mtype, DetShapeType(), *this);
562  DNekMatSharedPtr iprodmat = GetStdMatrix(iprodmatkey);
563 
564  Blas::Dgemv('N', m_ncoeffs, nq, 1.0, iprodmat->GetPtr().get(), m_ncoeffs,
565  inarray.get(), 1, 0.0, outarray.get(), 1);
566 }
567 
569  const int dir, const Array<OneD, const NekDouble> &inarray,
570  Array<OneD, NekDouble> &outarray)
571 {
572  int i;
573  int nquad0 = m_base[0]->GetNumPoints();
574  int nquad1 = m_base[1]->GetNumPoints();
575  int nqtot = nquad0 * nquad1;
576  int nmodes0 = m_base[0]->GetNumModes();
577  int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
578 
579  Array<OneD, NekDouble> gfac0(2 * wspsize);
580  Array<OneD, NekDouble> tmp0(gfac0 + wspsize);
581 
582  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
583 
584  // set up geometric factor: 2/(1-z1)
585  for (i = 0; i < nquad1; ++i)
586  {
587  gfac0[i] = 2.0 / (1 - z1[i]);
588  }
589 
590  for (i = 0; i < nquad1; ++i)
591  {
592  Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
593  &tmp0[0] + i * nquad0, 1);
594  }
595 
596  MultiplyByQuadratureMetric(tmp0, tmp0);
597 
598  switch (dir)
599  {
600  case 0:
601  {
602  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
603  m_base[1]->GetBdata(), tmp0, outarray,
604  gfac0);
605  break;
606  }
607  case 1:
608  {
610  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
611 
612  for (i = 0; i < nquad0; ++i)
613  {
614  gfac0[i] = 0.5 * (1 + z0[i]);
615  }
616 
617  for (i = 0; i < nquad1; ++i)
618  {
619  Vmath::Vmul(nquad0, &gfac0[0], 1, &tmp0[0] + i * nquad0, 1,
620  &tmp0[0] + i * nquad0, 1);
621  }
622 
623  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(),
624  m_base[1]->GetBdata(), tmp0, tmp3,
625  gfac0);
626 
627  MultiplyByQuadratureMetric(inarray, tmp0);
628  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
629  m_base[1]->GetDbdata(), tmp0, outarray,
630  gfac0);
631  Vmath::Vadd(m_ncoeffs, &tmp3[0], 1, &outarray[0], 1, &outarray[0],
632  1);
633  break;
634  }
635  default:
636  {
637  ASSERTL1(false, "input dir is out of range");
638  break;
639  }
640  }
641 }
642 
643 //---------------------------------------
644 // Evaluation functions
645 //---------------------------------------
646 
649 {
650  NekDouble d1 = 1. - xi[1];
651  if (fabs(d1) < NekConstants::kNekZeroTol)
652  {
653  if (d1 >= 0.)
654  {
656  }
657  else
658  {
660  }
661  }
662  eta[0] = 2. * (1. + xi[0]) / d1 - 1.0;
663  eta[1] = xi[1];
664 }
665 
668 {
669  xi[0] = (1.0 + eta[0]) * (1.0 - eta[1]) * 0.5 - 1.0;
670  xi[1] = eta[1];
671 }
672 
673 void StdTriExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
674 {
675  int i, m;
676  int nquad0 = m_base[0]->GetNumPoints();
677  int nquad1 = m_base[1]->GetNumPoints();
678  int order0 = m_base[0]->GetNumModes();
679  int order1 = m_base[1]->GetNumModes();
680  int mode0 = 0;
681  Array<OneD, const NekDouble> base0 = m_base[0]->GetBdata();
682  Array<OneD, const NekDouble> base1 = m_base[1]->GetBdata();
683 
684  ASSERTL2(mode <= m_ncoeffs, "calling argument mode is larger than "
685  "total expansion order");
686 
687  m = order1;
688  for (i = 0; i < order0; ++i, m += order1 - i)
689  {
690  if (m > mode)
691  {
692  mode0 = i;
693  break;
694  }
695  }
696 
697  // deal with top vertex mode in modified basis
698  if (mode == 1 && m_base[0]->GetBasisType() == LibUtilities::eModified_A)
699  {
700  Vmath::Fill(nquad0 * nquad1, 1.0, outarray, 1);
701  }
702  else
703  {
704  for (i = 0; i < nquad1; ++i)
705  {
706  Vmath::Vcopy(nquad0, (NekDouble *)(base0.get() + mode0 * nquad0), 1,
707  &outarray[0] + i * nquad0, 1);
708  }
709  }
710 
711  for (i = 0; i < nquad0; ++i)
712  {
713  Vmath::Vmul(nquad1, (NekDouble *)(base1.get() + mode * nquad1), 1,
714  &outarray[0] + i, nquad0, &outarray[0] + i, nquad0);
715  }
716 }
717 
719  const Array<OneD, const NekDouble> &coords, int mode)
720 {
721  Array<OneD, NekDouble> coll(2);
722  LocCoordToLocCollapsed(coords, coll);
723 
724  // From mode we need to determine mode0 and mode1 in the (p,q)
725  // direction. mode1 can be directly inferred from mode.
726  const int nm1 = m_base[1]->GetNumModes();
727  const double c = 1 + 2 * nm1;
728  const int mode0 = floor(0.5 * (c - sqrt(c * c - 8 * mode)));
729 
730  if (mode == 1 && m_base[0]->GetBasisType() == LibUtilities::eModified_A)
731  {
732  // Account for collapsed vertex.
733  return StdExpansion::BaryEvaluateBasis<1>(coll[1], 1);
734  }
735  else
736  {
737  return StdExpansion::BaryEvaluateBasis<0>(coll[0], mode0) *
738  StdExpansion::BaryEvaluateBasis<1>(coll[1], mode);
739  }
740 }
741 
743 {
744  return 3;
745 }
746 
748 {
749  return 3;
750 }
751 
753 {
755 }
756 
758 {
760  "BasisType is not a boundary interior form");
762  "BasisType is not a boundary interior form");
763 
764  return 3 + (GetBasisNumModes(0) - 2) + 2 * (GetBasisNumModes(1) - 2);
765 }
766 
768 {
770  "BasisType is not a boundary interior form");
772  "BasisType is not a boundary interior form");
773 
774  return GetBasisNumModes(0) + 2 * GetBasisNumModes(1);
775 }
776 
777 int StdTriExp::v_GetTraceNcoeffs(const int i) const
778 {
779  ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
780 
781  if (i == 0)
782  {
783  return GetBasisNumModes(0);
784  }
785  else
786  {
787  return GetBasisNumModes(1);
788  }
789 }
790 
791 int StdTriExp::v_GetTraceNumPoints(const int i) const
792 {
793  ASSERTL2((i >= 0) && (i <= 2), "edge id is out of range");
794 
795  if (i == 0)
796  {
797  return GetNumPoints(0);
798  }
799  else
800  {
801  return GetNumPoints(1);
802  }
803 }
804 
806  const std::vector<unsigned int> &nummodes, int &modes_offset)
807 {
809  nummodes[modes_offset], nummodes[modes_offset + 1]);
810  modes_offset += 2;
811 
812  return nmodes;
813 }
814 
816  Array<OneD, NekDouble> &coords_1,
817  Array<OneD, NekDouble> &coords_2)
818 {
819  boost::ignore_unused(coords_2);
820 
821  Array<OneD, const NekDouble> z0 = m_base[0]->GetZ();
822  Array<OneD, const NekDouble> z1 = m_base[1]->GetZ();
823  int nq0 = GetNumPoints(0);
824  int nq1 = GetNumPoints(1);
825  int i, j;
826 
827  for (i = 0; i < nq1; ++i)
828  {
829  for (j = 0; j < nq0; ++j)
830  {
831  coords_0[i * nq0 + j] = (1 + z0[j]) * (1 - z1[i]) / 2.0 - 1.0;
832  }
833  Vmath::Fill(nq0, z1[i], &coords_1[0] + i * nq0, 1);
834  }
835 }
836 
838 {
839  return m_base[0]->GetBasisType() == LibUtilities::eModified_A &&
840  m_base[1]->GetBasisType() == LibUtilities::eModified_B;
841 }
842 
844  const int j) const
845 {
846  boost::ignore_unused(j);
847  ASSERTL2(i >= 0 && i <= 2, "edge id is out of range");
848 
849  if (i == 0)
850  {
851  return GetBasis(0)->GetBasisKey();
852  }
853  else
854  {
855  switch (m_base[1]->GetBasisType())
856  {
858  {
859  switch (m_base[1]->GetPointsType())
860  {
861  case LibUtilities::eGaussRadauMAlpha1Beta0:
862  {
864  m_base[1]
865  ->GetBasisKey()
866  .GetPointsKey()
867  .GetNumPoints() +
868  1,
871  m_base[1]->GetNumModes(),
872  pkey);
873  break;
874  }
875 
876  break;
877  // Currently this does not increase the points by
878  // 1 since when using this quadrature we are
879  // presuming it is already been increased by one
880  // when comopared to
881  // GaussRadauMAlpha1Beta0. Currently used in the
882  // GJP option
883  //
884  // Note have put down it back to numpoints +1 to
885  // test for use on tri faces and GJP.
886 
888  {
890  m_base[1]
891  ->GetBasisKey()
892  .GetPointsKey()
893  .GetNumPoints() +
894  1,
897  m_base[1]->GetNumModes(),
898  pkey);
899  break;
900  }
901 
903  {
905  m_base[1]
906  ->GetBasisKey()
907  .GetPointsKey()
908  .GetNumPoints() +
909  1,
912  m_base[1]->GetNumModes(),
913  pkey);
914  break;
915  }
916 
917  break;
918  default:
920  "Unexpected points distribution " +
922  [m_base[1]->GetPointsType()] +
923  " in StdTriExp::v_GetTraceBasisKey");
924  break;
925  }
926  break;
927  }
928  default:
930  "Information not available to set edge key");
931  break;
932  }
933  }
935 }
936 
937 //--------------------------
938 // Mappings
939 //--------------------------
940 
941 int StdTriExp::v_GetVertexMap(const int localVertexId, bool useCoeffPacking)
942 {
945  "Mapping not defined for this type of basis");
946 
947  int localDOF = 0;
948  if (useCoeffPacking == true)
949  {
950  switch (localVertexId)
951  {
952  case 0:
953  {
954  localDOF = 0;
955  break;
956  }
957  case 1:
958  {
959  localDOF = 1;
960  break;
961  }
962  case 2:
963  {
964  localDOF = m_base[1]->GetNumModes();
965  break;
966  }
967  default:
968  {
969  ASSERTL0(false, "eid must be between 0 and 2");
970  break;
971  }
972  }
973  }
974  else // follow book format for vertex indexing.
975  {
976  switch (localVertexId)
977  {
978  case 0:
979  {
980  localDOF = 0;
981  break;
982  }
983  case 1:
984  {
985  localDOF = m_base[1]->GetNumModes();
986  break;
987  }
988  case 2:
989  {
990  localDOF = 1;
991  break;
992  }
993  default:
994  {
995  ASSERTL0(false, "eid must be between 0 and 2");
996  break;
997  }
998  }
999  }
1000 
1001  return localDOF;
1002 }
1003 
1005 {
1008  "Expansion not of a proper type");
1009 
1010  int i, j;
1011  int cnt = 0;
1012  int nummodes0, nummodes1;
1013  int startvalue;
1014  if (outarray.size() != GetNcoeffs() - NumBndryCoeffs())
1015  {
1017  }
1018 
1019  nummodes0 = m_base[0]->GetNumModes();
1020  nummodes1 = m_base[1]->GetNumModes();
1021 
1022  startvalue = 2 * nummodes1;
1023 
1024  for (i = 0; i < nummodes0 - 2; i++)
1025  {
1026  for (j = 0; j < nummodes1 - 3 - i; j++)
1027  {
1028  outarray[cnt++] = startvalue + j;
1029  }
1030  startvalue += nummodes1 - 2 - i;
1031  }
1032 }
1033 
1035 {
1038  "Expansion not of expected type");
1039  int i;
1040  int cnt;
1041  int nummodes0, nummodes1;
1042  int value;
1043 
1044  if (outarray.size() != NumBndryCoeffs())
1045  {
1047  }
1048 
1049  nummodes0 = m_base[0]->GetNumModes();
1050  nummodes1 = m_base[1]->GetNumModes();
1051 
1052  value = 2 * nummodes1 - 1;
1053  for (i = 0; i < value; i++)
1054  {
1055  outarray[i] = i;
1056  }
1057  cnt = value;
1058 
1059  for (i = 0; i < nummodes0 - 2; i++)
1060  {
1061  outarray[cnt++] = value;
1062  value += nummodes1 - 2 - i;
1063  }
1064 }
1065 
1066 void StdTriExp::v_GetTraceCoeffMap(const unsigned int eid,
1067  Array<OneD, unsigned int> &maparray)
1068 {
1069 
1072  "Mapping not defined for this type of basis");
1073 
1074  ASSERTL1(eid < 3, "eid must be between 0 and 2");
1075 
1076  int i;
1077  int order0 = m_base[0]->GetNumModes();
1078  int order1 = m_base[1]->GetNumModes();
1079  int numModes = (eid == 0) ? order0 : order1;
1080 
1081  if (maparray.size() != numModes)
1082  {
1083  maparray = Array<OneD, unsigned int>(numModes);
1084  }
1085 
1086  switch (eid)
1087  {
1088  case 0:
1089  {
1090  int cnt = 0;
1091  for (i = 0; i < numModes; cnt += order1 - i, ++i)
1092  {
1093  maparray[i] = cnt;
1094  }
1095  break;
1096  }
1097  case 1:
1098  {
1099  maparray[0] = order1;
1100  maparray[1] = 1;
1101  for (i = 2; i < numModes; i++)
1102  {
1103  maparray[i] = order1 - 1 + i;
1104  }
1105  break;
1106  }
1107  case 2:
1108  {
1109  for (i = 0; i < numModes; i++)
1110  {
1111  maparray[i] = i;
1112  }
1113  break;
1114  }
1115  default:
1116  ASSERTL0(false, "eid must be between 0 and 2");
1117  break;
1118  }
1119 }
1120 
1122  const int eid, Array<OneD, unsigned int> &maparray,
1123  Array<OneD, int> &signarray, const Orientation edgeOrient)
1124 {
1127  "Mapping not defined for this type of basis");
1128  int i;
1129  const int nummodes1 = m_base[1]->GetNumModes();
1130  const int nEdgeIntCoeffs = GetTraceNcoeffs(eid) - 2;
1131 
1132  if (maparray.size() != nEdgeIntCoeffs)
1133  {
1134  maparray = Array<OneD, unsigned int>(nEdgeIntCoeffs);
1135  }
1136 
1137  if (signarray.size() != nEdgeIntCoeffs)
1138  {
1139  signarray = Array<OneD, int>(nEdgeIntCoeffs, 1);
1140  }
1141  else
1142  {
1143  fill(signarray.get(), signarray.get() + nEdgeIntCoeffs, 1);
1144  }
1145 
1146  switch (eid)
1147  {
1148  case 0:
1149  {
1150  int cnt = 2 * nummodes1 - 1;
1151  for (i = 0; i < nEdgeIntCoeffs; cnt += nummodes1 - 2 - i, ++i)
1152  {
1153  maparray[i] = cnt;
1154  }
1155  break;
1156  }
1157  case 1:
1158  {
1159  for (i = 0; i < nEdgeIntCoeffs; i++)
1160  {
1161  maparray[i] = nummodes1 + 1 + i;
1162  }
1163  break;
1164  }
1165  case 2:
1166  {
1167  for (i = 0; i < nEdgeIntCoeffs; i++)
1168  {
1169  maparray[i] = 2 + i;
1170  }
1171  break;
1172  }
1173  default:
1174  {
1175  ASSERTL0(false, "eid must be between 0 and 2");
1176  break;
1177  }
1178  }
1179 
1180  if (edgeOrient == eBackwards)
1181  {
1182  for (i = 1; i < nEdgeIntCoeffs; i += 2)
1183  {
1184  signarray[i] = -1;
1185  }
1186  }
1187 }
1188 
1189 //---------------------------------------
1190 // Wrapper functions
1191 //---------------------------------------
1192 
1194 {
1195 
1196  MatrixType mtype = mkey.GetMatrixType();
1197 
1198  DNekMatSharedPtr Mat;
1199 
1200  switch (mtype)
1201  {
1203  {
1204  int nq0, nq1, nq;
1205 
1206  nq0 = m_base[0]->GetNumPoints();
1207  nq1 = m_base[1]->GetNumPoints();
1208 
1209  // take definition from key
1210  if (mkey.ConstFactorExists(eFactorConst))
1211  {
1212  nq = (int)mkey.GetConstFactor(eFactorConst);
1213  }
1214  else
1215  {
1216  nq = max(nq0, nq1);
1217  }
1218 
1221  Array<OneD, NekDouble> coll(2);
1223  Array<OneD, NekDouble> tmp(nq0);
1224 
1225  Mat = MemoryManager<DNekMat>::AllocateSharedPtr(neq, nq0 * nq1);
1226  int cnt = 0;
1227 
1228  for (int i = 0; i < nq; ++i)
1229  {
1230  for (int j = 0; j < nq - i; ++j, ++cnt)
1231  {
1232  coords[cnt] = Array<OneD, NekDouble>(2);
1233  coords[cnt][0] = -1.0 + 2 * j / (NekDouble)(nq - 1);
1234  coords[cnt][1] = -1.0 + 2 * i / (NekDouble)(nq - 1);
1235  }
1236  }
1237 
1238  for (int i = 0; i < neq; ++i)
1239  {
1240  LocCoordToLocCollapsed(coords[i], coll);
1241 
1242  I[0] = m_base[0]->GetI(coll);
1243  I[1] = m_base[1]->GetI(coll + 1);
1244 
1245  // interpolate first coordinate direction
1246  for (int j = 0; j < nq1; ++j)
1247  {
1248  NekDouble fac = (I[1]->GetPtr())[j];
1249  Vmath::Smul(nq0, fac, I[0]->GetPtr(), 1, tmp, 1);
1250 
1251  Vmath::Vcopy(nq0, &tmp[0], 1,
1252  Mat->GetRawPtr() + j * nq0 * neq + i, neq);
1253  }
1254  }
1255  break;
1256  }
1257  default:
1258  {
1260  break;
1261  }
1262  }
1263 
1264  return Mat;
1265 }
1266 
1268 {
1269  return v_GenMatrix(mkey);
1270 }
1271 
1272 //---------------------------------------
1273 // Operator evaluation functions
1274 //---------------------------------------
1275 
1277  Array<OneD, NekDouble> &outarray,
1278  const StdMatrixKey &mkey)
1279 {
1280  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1281 }
1282 
1284  Array<OneD, NekDouble> &outarray,
1285  const StdMatrixKey &mkey)
1286 {
1287  StdTriExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1288 }
1289 
1290 void StdTriExp::v_LaplacianMatrixOp(const int k1, const int k2,
1291  const Array<OneD, const NekDouble> &inarray,
1292  Array<OneD, NekDouble> &outarray,
1293  const StdMatrixKey &mkey)
1294 {
1295  StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1296 }
1297 
1299  const Array<OneD, const NekDouble> &inarray,
1300  Array<OneD, NekDouble> &outarray,
1301  const StdMatrixKey &mkey)
1302 {
1303  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1304 }
1305 
1307  Array<OneD, NekDouble> &outarray,
1308  const StdMatrixKey &mkey)
1309 {
1310  StdTriExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1311 }
1312 
1314  const StdMatrixKey &mkey)
1315 {
1316  int qa = m_base[0]->GetNumPoints();
1317  int qb = m_base[1]->GetNumPoints();
1318  int nmodes_a = m_base[0]->GetNumModes();
1319  int nmodes_b = m_base[1]->GetNumModes();
1320 
1321  // Declare orthogonal basis.
1324 
1327  StdTriExp OrthoExp(Ba, Bb);
1328 
1329  Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
1330 
1331  // project onto physical space.
1332  OrthoExp.FwdTrans(array, orthocoeffs);
1333 
1334  if (mkey.ConstFactorExists(
1335  eFactorSVVPowerKerDiffCoeff)) // Rodrigo's power kern
1336  {
1338  NekDouble SvvDiffCoeff =
1341 
1342  int cnt = 0;
1343  for (int j = 0; j < nmodes_a; ++j)
1344  {
1345  for (int k = 0; k < nmodes_b - j; ++k, ++cnt)
1346  {
1347  NekDouble fac = std::max(
1348  pow((1.0 * j) / (nmodes_a - 1), cutoff * nmodes_a),
1349  pow((1.0 * k) / (nmodes_b - 1), cutoff * nmodes_b));
1350 
1351  orthocoeffs[cnt] *= (SvvDiffCoeff * fac);
1352  }
1353  }
1354  }
1355  else if (mkey.ConstFactorExists(
1356  eFactorSVVDGKerDiffCoeff)) // Rodrigo/mansoor's DG kernel
1357  {
1358  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDGKerDiffCoeff) *
1360  int max_ab = max(nmodes_a - kSVVDGFiltermodesmin,
1361  nmodes_b - kSVVDGFiltermodesmin);
1362  max_ab = max(max_ab, 0);
1363  max_ab = min(max_ab, kSVVDGFiltermodesmax - kSVVDGFiltermodesmin);
1364 
1365  int cnt = 0;
1366  for (int j = 0; j < nmodes_a; ++j)
1367  {
1368  for (int k = 0; k < nmodes_b - j; ++k, ++cnt)
1369  {
1370  int maxjk = max(j, k);
1371  maxjk = min(maxjk, kSVVDGFiltermodesmax - 1);
1372 
1373  orthocoeffs[cnt] *= SvvDiffCoeff * kSVVDGFilter[max_ab][maxjk];
1374  }
1375  }
1376  }
1377  else
1378  {
1379  NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
1380 
1381  int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio) *
1382  min(nmodes_a, nmodes_b));
1383 
1384  NekDouble epsilon = 1.0;
1385  int nmodes = min(nmodes_a, nmodes_b);
1386 
1387  int cnt = 0;
1388 
1389  // apply SVV filter (JEL)
1390  for (int j = 0; j < nmodes_a; ++j)
1391  {
1392  for (int k = 0; k < nmodes_b - j; ++k)
1393  {
1394  if (j + k >= cutoff)
1395  {
1396  orthocoeffs[cnt] *=
1397  (SvvDiffCoeff *
1398  exp(-(j + k - nmodes) * (j + k - nmodes) /
1399  ((NekDouble)((j + k - cutoff + epsilon) *
1400  (j + k - cutoff + epsilon)))));
1401  }
1402  else
1403  {
1404  orthocoeffs[cnt] *= 0.0;
1405  }
1406  cnt++;
1407  }
1408  }
1409  }
1410 
1411  // backward transform to physical space
1412  OrthoExp.BwdTrans(orthocoeffs, array);
1413 }
1414 
1416  const Array<OneD, const NekDouble> &inarray,
1417  Array<OneD, NekDouble> &outarray)
1418 {
1419  int n_coeffs = inarray.size();
1420  int nquad0 = m_base[0]->GetNumPoints();
1421  int nquad1 = m_base[1]->GetNumPoints();
1422  Array<OneD, NekDouble> coeff(n_coeffs);
1423  Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
1426  int nqtot = nquad0 * nquad1;
1427  Array<OneD, NekDouble> phys_tmp(nqtot, 0.0);
1428 
1429  int nmodes0 = m_base[0]->GetNumModes();
1430  int nmodes1 = m_base[1]->GetNumModes();
1431  int numMin2 = nmodes0;
1432  int i;
1433 
1434  const LibUtilities::PointsKey Pkey0(nmodes0,
1436  const LibUtilities::PointsKey Pkey1(nmodes1,
1438 
1439  LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1440  LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1441 
1442  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1443  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B, nmodes1, Pkey1);
1444 
1445  StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1447 
1450  bortho0, bortho1);
1451 
1452  m_TriExp->BwdTrans(inarray, phys_tmp);
1453  m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
1454 
1455  for (i = 0; i < n_coeffs; i++)
1456  {
1457  if (i == numMin)
1458  {
1459  coeff[i] = 0.0;
1460  numMin += numMin2 - 1;
1461  numMin2 -= 1.0;
1462  }
1463  }
1464 
1465  m_OrthoTriExp->BwdTrans(coeff, phys_tmp);
1466  m_TriExp->FwdTrans(phys_tmp, outarray);
1467 }
1468 
1470  const Array<OneD, const NekDouble> &inarray,
1471  Array<OneD, NekDouble> &outarray, const StdMatrixKey &mkey)
1472 {
1474 
1475  if (inarray.get() == outarray.get())
1476  {
1478  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, tmp.get(), 1);
1479 
1480  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1481  m_ncoeffs, tmp.get(), 1, 0.0, outarray.get(), 1);
1482  }
1483  else
1484  {
1485  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, 1.0, mat->GetPtr().get(),
1486  m_ncoeffs, inarray.get(), 1, 0.0, outarray.get(), 1);
1487  }
1488 }
1489 
1490 //---------------------------------------
1491 // Private helper functions
1492 //---------------------------------------
1493 
1495  const Array<OneD, const NekDouble> &inarray,
1496  Array<OneD, NekDouble> &outarray)
1497 {
1498  int i;
1499  int nquad0 = m_base[0]->GetNumPoints();
1500  int nquad1 = m_base[1]->GetNumPoints();
1501 
1502  const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
1503  const Array<OneD, const NekDouble> &w1 = m_base[1]->GetW();
1504  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1505 
1506  // multiply by integration constants
1507  for (i = 0; i < nquad1; ++i)
1508  {
1509  Vmath::Vmul(nquad0, inarray.get() + i * nquad0, 1, w0.get(), 1,
1510  outarray.get() + i * nquad0, 1);
1511  }
1512 
1513  switch (m_base[1]->GetPointsType())
1514  {
1515  // (1,0) Jacobi Inner product
1516  case LibUtilities::eGaussRadauMAlpha1Beta0:
1517  for (i = 0; i < nquad1; ++i)
1518  {
1519  Blas::Dscal(nquad0, 0.5 * w1[i], outarray.get() + i * nquad0,
1520  1);
1521  }
1522  break;
1523  // Legendre inner product
1524  default:
1525  for (i = 0; i < nquad1; ++i)
1526  {
1527  Blas::Dscal(nquad0, 0.5 * (1 - z1[i]) * w1[i],
1528  outarray.get() + i * nquad0, 1);
1529  }
1530  break;
1531  }
1532 }
1533 
1535  bool standard)
1536 {
1537  boost::ignore_unused(standard);
1538 
1539  int np1 = m_base[0]->GetNumPoints();
1540  int np2 = m_base[1]->GetNumPoints();
1541  int np = max(np1, np2);
1542 
1543  conn = Array<OneD, int>(3 * (np - 1) * (np - 1));
1544 
1545  int row = 0;
1546  int rowp1 = 0;
1547  int cnt = 0;
1548  for (int i = 0; i < np - 1; ++i)
1549  {
1550  rowp1 += np - i;
1551  for (int j = 0; j < np - i - 2; ++j)
1552  {
1553  conn[cnt++] = row + j;
1554  conn[cnt++] = row + j + 1;
1555  conn[cnt++] = rowp1 + j;
1556 
1557  conn[cnt++] = rowp1 + j + 1;
1558  conn[cnt++] = rowp1 + j;
1559  conn[cnt++] = row + j + 1;
1560  }
1561 
1562  conn[cnt++] = row + np - i - 2;
1563  conn[cnt++] = row + np - i - 1;
1564  conn[cnt++] = rowp1 + np - i - 2;
1565 
1566  row += np - i;
1567  }
1568 }
1569 } // namespace StdRegions
1570 } // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:215
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:209
#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
int GetNumModes() const
Returns the order of the basis.
Definition: Basis.h:83
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 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 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
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)
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 void v_GetCoords(Array< OneD, NekDouble > &coords_x, Array< OneD, NekDouble > &coords_y, Array< OneD, NekDouble > &coords_z)
Definition: StdTriExp.cpp:815
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: StdTriExp.cpp:256
virtual void v_GetTraceCoeffMap(const unsigned int traceid, Array< OneD, unsigned int > &maparray)
Definition: StdTriExp.cpp:1066
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:313
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1469
virtual void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true)
Definition: StdTriExp.cpp:1534
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Definition: StdTriExp.cpp:78
virtual int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false)
Definition: StdTriExp.cpp:941
virtual int v_GetTraceNcoeffs(const int i) const
Definition: StdTriExp.cpp:777
virtual int v_NumBndryCoeffs() const
Definition: StdTriExp.cpp:757
virtual DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1193
virtual int v_GetTraceNumPoints(const int i) const
Definition: StdTriExp.cpp:791
virtual void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:246
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1298
virtual void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi)
Definition: StdTriExp.cpp:666
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: StdTriExp.cpp:215
virtual void v_GetTraceInteriorToElementMap(const int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, const Orientation edgeOrient=eForwards)
Definition: StdTriExp.cpp:1121
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1306
virtual void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray)
Definition: StdTriExp.cpp:1034
virtual int v_GetNtraces() const
Definition: StdTriExp.cpp:747
virtual void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta)
Definition: StdTriExp.cpp:647
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:450
virtual void v_GetInteriorMap(Array< OneD, unsigned int > &outarray)
Definition: StdTriExp.cpp:1004
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1313
virtual int v_NumDGBndryCoeffs() const
Definition: StdTriExp.cpp:767
virtual int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset)
Definition: StdTriExp.cpp:805
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[p]*base1[pq] and put into ou...
Definition: StdTriExp.cpp:444
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Backward tranform for triangular elements.
Definition: StdTriExp.cpp:240
virtual int v_GetNverts() const
Definition: StdTriExp.cpp:742
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: StdTriExp.cpp:462
virtual bool v_IsBoundaryInteriorExpansion()
Definition: StdTriExp.cpp:837
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:535
virtual void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:1494
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:1415
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: StdTriExp.cpp:297
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: StdTriExp.cpp:126
virtual LibUtilities::ShapeType v_DetShapeType() const
Definition: StdTriExp.cpp:752
virtual void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:673
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdTriExp.cpp:718
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1283
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: StdTriExp.cpp:489
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:568
virtual const LibUtilities::BasisKey v_GetTraceBasisKey(const int i, const int j) const
Definition: StdTriExp.cpp:843
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1276
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdTriExp.cpp:528
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey)
Definition: StdTriExp.cpp:1267
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 Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
Definition: Blas.hpp:168
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
static void Daxpy(const int &n, const double &alpha, const double *x, const int &incx, const double *y, const int &incy)
BLAS level 1: y = alpha x plus y.
Definition: Blas.hpp:154
int getNumberOfCoefficients(int Na, int Nb)
Definition: ShapeType.hpp:112
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
const std::string kPointsTypeStr[]
Definition: Foundations.hpp:54
@ eGaussRadauMLegendre
1D Gauss-Radau-Legendre quadrature points, pinned at x=-1
Definition: PointsType.h:49
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:75
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:51
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:44
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:46
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
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::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition: StdTriExp.h:235
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 Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.cpp:574
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
Definition: Vmath.cpp:359
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 Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/y.
Definition: Vmath.cpp:324
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha - x.
Definition: Vmath.cpp:384
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
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:291