Nektar++
TriExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File TriExp.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 // The above copyright notice and this permission notice shall be included
20 // in all copies or substantial portions of the Software.
21 //
22 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 // DEALINGS IN THE SOFTWARE.
29 //
30 // Description: Expasion for triangular elements.
31 //
32 ///////////////////////////////////////////////////////////////////////////////
33 
34 #include <boost/core/ignore_unused.hpp>
35 
39 #include <LocalRegions/SegExp.h>
40 #include <LocalRegions/TriExp.h>
42 
43 using namespace std;
44 
45 namespace Nektar
46 {
47 namespace LocalRegions
48 {
49 TriExp::TriExp(const LibUtilities::BasisKey &Ba,
50  const LibUtilities::BasisKey &Bb,
52  : StdExpansion(LibUtilities::StdTriData::getNumberOfCoefficients(
53  Ba.GetNumModes(), (Bb.GetNumModes())),
54  2, Ba, Bb),
55  StdExpansion2D(LibUtilities::StdTriData::getNumberOfCoefficients(
56  Ba.GetNumModes(), (Bb.GetNumModes())),
57  Ba, Bb),
58  StdTriExp(Ba, Bb), Expansion(geom), Expansion2D(geom),
59  m_matrixManager(
60  std::bind(&Expansion2D::CreateMatrix, this, std::placeholders::_1),
61  std::string("TriExpMatrix")),
62  m_staticCondMatrixManager(std::bind(&Expansion::CreateStaticCondMatrix,
63  this, std::placeholders::_1),
64  std::string("TriExpStaticCondMatrix"))
65 {
66 }
67 
69  : StdExpansion(T), StdExpansion2D(T), StdTriExp(T), Expansion(T),
70  Expansion2D(T), m_matrixManager(T.m_matrixManager),
71  m_staticCondMatrixManager(T.m_staticCondMatrixManager)
72 {
73 }
74 
76 {
77 }
78 
80 {
81  int nquad0 = m_base[0]->GetNumPoints();
82  int nquad1 = m_base[1]->GetNumPoints();
84  NekDouble ival;
85  Array<OneD, NekDouble> tmp(nquad0 * nquad1);
86 
87  // multiply inarray with Jacobian
88  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
89  {
90  Vmath::Vmul(nquad0 * nquad1, jac, 1, inarray, 1, tmp, 1);
91  }
92  else
93  {
94  Vmath::Smul(nquad0 * nquad1, jac[0], inarray, 1, tmp, 1);
95  }
96 
97  // call StdQuadExp version;
98  ival = StdTriExp::v_Integral(tmp);
99  return ival;
100 }
101 
103  Array<OneD, NekDouble> &out_d0,
104  Array<OneD, NekDouble> &out_d1,
105  Array<OneD, NekDouble> &out_d2)
106 {
107  int nquad0 = m_base[0]->GetNumPoints();
108  int nquad1 = m_base[1]->GetNumPoints();
109  int nqtot = nquad0 * nquad1;
110  const Array<TwoD, const NekDouble> &df =
111  m_metricinfo->GetDerivFactors(GetPointsKeys());
112 
113  Array<OneD, NekDouble> diff0(2 * nqtot);
114  Array<OneD, NekDouble> diff1(diff0 + nqtot);
115 
116  StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
117 
118  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
119  {
120  if (out_d0.size())
121  {
122  Vmath::Vmul(nqtot, df[0], 1, diff0, 1, out_d0, 1);
123  Vmath::Vvtvp(nqtot, df[1], 1, diff1, 1, out_d0, 1, out_d0, 1);
124  }
125 
126  if (out_d1.size())
127  {
128  Vmath::Vmul(nqtot, df[2], 1, diff0, 1, out_d1, 1);
129  Vmath::Vvtvp(nqtot, df[3], 1, diff1, 1, out_d1, 1, out_d1, 1);
130  }
131 
132  if (out_d2.size())
133  {
134  Vmath::Vmul(nqtot, df[4], 1, diff0, 1, out_d2, 1);
135  Vmath::Vvtvp(nqtot, df[5], 1, diff1, 1, out_d2, 1, out_d2, 1);
136  }
137  }
138  else // regular geometry
139  {
140  if (out_d0.size())
141  {
142  Vmath::Smul(nqtot, df[0][0], diff0, 1, out_d0, 1);
143  Blas::Daxpy(nqtot, df[1][0], diff1, 1, out_d0, 1);
144  }
145 
146  if (out_d1.size())
147  {
148  Vmath::Smul(nqtot, df[2][0], diff0, 1, out_d1, 1);
149  Blas::Daxpy(nqtot, df[3][0], diff1, 1, out_d1, 1);
150  }
151 
152  if (out_d2.size())
153  {
154  Vmath::Smul(nqtot, df[4][0], diff0, 1, out_d2, 1);
155  Blas::Daxpy(nqtot, df[5][0], diff1, 1, out_d2, 1);
156  }
157  }
158 }
159 
160 void TriExp::v_PhysDeriv(const int dir,
161  const Array<OneD, const NekDouble> &inarray,
162  Array<OneD, NekDouble> &outarray)
163 {
164  switch (dir)
165  {
166  case 0:
167  {
168  PhysDeriv(inarray, outarray, NullNekDouble1DArray,
170  }
171  break;
172  case 1:
173  {
174  PhysDeriv(inarray, NullNekDouble1DArray, outarray,
176  }
177  break;
178  case 2:
179  {
181  outarray);
182  }
183  break;
184  default:
185  {
186  ASSERTL1(false, "input dir is out of range");
187  }
188  break;
189  }
190 }
191 
193  const Array<OneD, const NekDouble> &inarray,
195 {
196  if (!out.size())
197  {
198  return;
199  }
200 
201  int nquad0 = m_base[0]->GetNumPoints();
202  int nquad1 = m_base[1]->GetNumPoints();
203  int nqtot = nquad0 * nquad1;
204 
205  const Array<TwoD, const NekDouble> &df =
206  m_metricinfo->GetDerivFactors(GetPointsKeys());
207 
208  Array<OneD, NekDouble> diff0(2 * nqtot);
209  Array<OneD, NekDouble> diff1(diff0 + nqtot);
210 
211  // diff0 = du/d_xi, diff1 = du/d_eta
212  StdTriExp::v_PhysDeriv(inarray, diff0, diff1);
213 
214  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
215  {
217 
218  // D^v_xi = v_x*d_xi/dx + v_y*d_xi/dy + v_z*d_xi/dz
219  // D^v_eta = v_x*d_eta/dx + v_y*d_eta/dy + v_z*d_eta/dz
220  for (int i = 0; i < 2; ++i)
221  {
222  tangmat[i] = Array<OneD, NekDouble>(nqtot, 0.0);
223  for (int k = 0; k < (m_geom->GetCoordim()); ++k)
224  {
225  Vmath::Vvtvp(nqtot, &df[2 * k + i][0], 1, &direction[k * nqtot],
226  1, &tangmat[i][0], 1, &tangmat[i][0], 1);
227  }
228  }
229 
230  /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
231  Vmath::Vmul(nqtot, &tangmat[0][0], 1, &diff0[0], 1, &out[0], 1);
232  Vmath::Vvtvp(nqtot, &tangmat[1][0], 1, &diff1[0], 1, &out[0], 1,
233  &out[0], 1);
234  }
235  else
236  {
238 
239  for (int i = 0; i < 2; ++i)
240  {
241  tangmat[i] = Array<OneD, NekDouble>(nqtot, 0.0);
242  for (int k = 0; k < (m_geom->GetCoordim()); ++k)
243  {
244  Vmath::Svtvp(nqtot, df[2 * k + i][0], &direction[k * nqtot], 1,
245  &tangmat[i][0], 1, &tangmat[i][0], 1);
246  }
247  }
248 
249  /// D_v = D^v_xi * du/d_xi + D^v_eta * du/d_eta
250  Vmath::Vmul(nqtot, &tangmat[0][0], 1, &diff0[0], 1, &out[0], 1);
251 
252  Vmath::Vvtvp(nqtot, &tangmat[1][0], 1, &diff1[0], 1, &out[0], 1,
253  &out[0], 1);
254  }
255 }
256 
258  Array<OneD, NekDouble> &outarray)
259 {
260  IProductWRTBase(inarray, outarray);
261 
262  // get Mass matrix inverse
263  MatrixKey masskey(StdRegions::eInvMass, DetShapeType(), *this);
264  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
265 
266  // copy inarray in case inarray == outarray
267  NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
268  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
269 
270  out = (*matsys) * in;
271 }
272 
274  const Array<OneD, const NekDouble> &inarray,
275  Array<OneD, NekDouble> &outarray)
276 {
277  int i, j;
278  int npoints[2] = {m_base[0]->GetNumPoints(), m_base[1]->GetNumPoints()};
279  int nmodes[2] = {m_base[0]->GetNumModes(), m_base[1]->GetNumModes()};
280 
281  fill(outarray.get(), outarray.get() + m_ncoeffs, 0.0);
282 
283  if (nmodes[0] == 1 && nmodes[1] == 1)
284  {
285  outarray[0] = inarray[0];
286  return;
287  }
288 
289  Array<OneD, NekDouble> physEdge[3];
290  Array<OneD, NekDouble> coeffEdge[3];
291  for (i = 0; i < 3; i++)
292  {
293  // define physEdge and add 1 so can interpolate grl10 points if
294  // necessary
295  physEdge[i] = Array<OneD, NekDouble>(max(npoints[i != 0], npoints[0]));
296  coeffEdge[i] = Array<OneD, NekDouble>(nmodes[i != 0]);
297  }
298 
299  for (i = 0; i < npoints[0]; i++)
300  {
301  physEdge[0][i] = inarray[i];
302  }
303 
304  // extract data in cartesian directions
305  for (i = 0; i < npoints[1]; i++)
306  {
307  physEdge[1][i] = inarray[npoints[0] - 1 + i * npoints[0]];
308  physEdge[2][i] = inarray[i * npoints[0]];
309  }
310 
311  SegExpSharedPtr segexp[3];
313  m_base[0]->GetBasisKey(), GetGeom2D()->GetEdge(0));
314 
316  {
317  for (i = 1; i < 3; i++)
318  {
320  m_base[i != 0]->GetBasisKey(), GetGeom2D()->GetEdge(i));
321  }
322  }
323  else // interploate using edge 0 GLL distribution
324  {
325  for (i = 1; i < 3; i++)
326  {
328  m_base[0]->GetBasisKey(), GetGeom2D()->GetEdge(i));
329 
330  LibUtilities::Interp1D(m_base[1]->GetPointsKey(), physEdge[i],
331  m_base[0]->GetPointsKey(), physEdge[i]);
332  }
333  npoints[1] = npoints[0];
334  }
335 
336  Array<OneD, unsigned int> mapArray;
337  Array<OneD, int> signArray;
338  NekDouble sign;
339  // define an orientation to get EdgeToElmtMapping from Cartesian data
340  StdRegions::Orientation orient[3] = {
342 
343  for (i = 0; i < 3; i++)
344  {
345  segexp[i]->FwdTransBndConstrained(physEdge[i], coeffEdge[i]);
346 
347  // this orient goes with the one above and so could
348  // probably set both to eForwards
349  GetTraceToElementMap(i, mapArray, signArray, orient[i]);
350  for (j = 0; j < nmodes[i != 0]; j++)
351  {
352  sign = (NekDouble)signArray[j];
353  outarray[mapArray[j]] = sign * coeffEdge[i][j];
354  }
355  }
356 
357  int nBoundaryDofs = NumBndryCoeffs();
358  int nInteriorDofs = m_ncoeffs - nBoundaryDofs;
359 
360  if (nInteriorDofs > 0)
361  {
364 
366  *this);
367  MassMatrixOp(outarray, tmp0, stdmasskey);
368  IProductWRTBase(inarray, tmp1);
369 
370  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
371 
372  // get Mass matrix inverse (only of interior DOF)
373  // use block (1,1) of the static condensed system
374  // note: this block alreay contains the inverse matrix
375  MatrixKey masskey(StdRegions::eMass, DetShapeType(), *this);
376  DNekScalMatSharedPtr matsys =
377  (m_staticCondMatrixManager[masskey])->GetBlock(1, 1);
378 
379  Array<OneD, NekDouble> rhs(nInteriorDofs);
380  Array<OneD, NekDouble> result(nInteriorDofs);
381 
382  GetInteriorMap(mapArray);
383 
384  for (i = 0; i < nInteriorDofs; i++)
385  {
386  rhs[i] = tmp1[mapArray[i]];
387  }
388 
389  Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, matsys->Scale(),
390  &((matsys->GetOwnedMatrix())->GetPtr())[0], nInteriorDofs,
391  rhs.get(), 1, 0.0, result.get(), 1);
392 
393  for (i = 0; i < nInteriorDofs; i++)
394  {
395  outarray[mapArray[i]] = result[i];
396  }
397  }
398 }
399 
401  Array<OneD, NekDouble> &outarray)
402 {
403  IProductWRTBase_SumFac(inarray, outarray);
404 }
405 
407  const Array<OneD, const NekDouble> &inarray,
408  Array<OneD, NekDouble> &outarray)
409 {
410  IProductWRTDerivBase_SumFac(dir, inarray, outarray);
411 }
412 
414  const Array<OneD, const NekDouble> &inarray,
415  Array<OneD, NekDouble> &outarray, bool multiplybyweights)
416 {
417  int nquad0 = m_base[0]->GetNumPoints();
418  int nquad1 = m_base[1]->GetNumPoints();
419  int order0 = m_base[0]->GetNumModes();
420 
421  if (multiplybyweights)
422  {
423  Array<OneD, NekDouble> tmp(nquad0 * nquad1 + nquad1 * order0);
424  Array<OneD, NekDouble> wsp(tmp + nquad0 * nquad1);
425 
426  MultiplyByQuadratureMetric(inarray, tmp);
427  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
428  m_base[1]->GetBdata(), tmp, outarray, wsp);
429  }
430  else
431  {
432  Array<OneD, NekDouble> wsp(+nquad1 * order0);
433 
434  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(),
435  m_base[1]->GetBdata(), inarray, outarray,
436  wsp);
437  }
438 }
439 
441  const Array<OneD, const NekDouble> &inarray,
442  Array<OneD, NekDouble> &outarray)
443 {
444  int nq = GetTotPoints();
445  MatrixKey iprodmatkey(StdRegions::eIProductWRTBase, DetShapeType(), *this);
446  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
447 
448  Blas::Dgemv('N', m_ncoeffs, nq, iprodmat->Scale(),
449  (iprodmat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
450  inarray.get(), 1, 0.0, outarray.get(), 1);
451 }
452 
454  const int dir, const Array<OneD, const NekDouble> &inarray,
455  Array<OneD, NekDouble> &outarray)
456 {
457  int nquad0 = m_base[0]->GetNumPoints();
458  int nquad1 = m_base[1]->GetNumPoints();
459  int nqtot = nquad0 * nquad1;
460  int nmodes0 = m_base[0]->GetNumModes();
461  int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
462 
463  Array<OneD, NekDouble> tmp0(4 * wspsize);
464  Array<OneD, NekDouble> tmp1(tmp0 + wspsize);
465  Array<OneD, NekDouble> tmp2(tmp0 + 2 * wspsize);
466  Array<OneD, NekDouble> tmp3(tmp0 + 3 * wspsize);
467 
469  tmp2D[0] = tmp1;
470  tmp2D[1] = tmp2;
471 
472  TriExp::v_AlignVectorToCollapsedDir(dir, inarray, tmp2D);
473 
474  MultiplyByQuadratureMetric(tmp1, tmp1);
475  MultiplyByQuadratureMetric(tmp2, tmp2);
476 
477  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
478  tmp1, tmp3, tmp0);
479  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
480  tmp2, outarray, tmp0);
481  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
482 }
483 
485  const int dir, const Array<OneD, const NekDouble> &inarray,
486  Array<OneD, Array<OneD, NekDouble>> &outarray)
487 {
488  ASSERTL1((dir == 0) || (dir == 1) || (dir == 2), "Invalid direction.");
489  ASSERTL1((dir == 2) ? (m_geom->GetCoordim() == 3) : true,
490  "Invalid direction.");
491 
492  int nquad0 = m_base[0]->GetNumPoints();
493  int nquad1 = m_base[1]->GetNumPoints();
494  int nqtot = nquad0 * nquad1;
495  int nmodes0 = m_base[0]->GetNumModes();
496  int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
497 
498  const Array<TwoD, const NekDouble> &df =
499  m_metricinfo->GetDerivFactors(GetPointsKeys());
500 
501  Array<OneD, NekDouble> tmp0(wspsize);
502  Array<OneD, NekDouble> tmp3(wspsize);
503  Array<OneD, NekDouble> gfac0(wspsize);
504  Array<OneD, NekDouble> gfac1(wspsize);
505 
506  Array<OneD, NekDouble> tmp1 = outarray[0];
507  Array<OneD, NekDouble> tmp2 = outarray[1];
508 
509  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
510  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
511 
512  // set up geometric factor: 2/(1-z1)
513  for (int i = 0; i < nquad1; ++i)
514  {
515  gfac0[i] = 2.0 / (1 - z1[i]);
516  }
517  for (int i = 0; i < nquad0; ++i)
518  {
519  gfac1[i] = 0.5 * (1 + z0[i]);
520  }
521 
522  for (int i = 0; i < nquad1; ++i)
523  {
524  Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
525  &tmp0[0] + i * nquad0, 1);
526  }
527 
528  for (int i = 0; i < nquad1; ++i)
529  {
530  Vmath::Vmul(nquad0, &gfac1[0], 1, &tmp0[0] + i * nquad0, 1,
531  &tmp1[0] + i * nquad0, 1);
532  }
533 
534  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
535  {
536  Vmath::Vmul(nqtot, &df[2 * dir][0], 1, &tmp0[0], 1, &tmp0[0], 1);
537  Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, &tmp1[0], 1, &tmp1[0], 1);
538  Vmath::Vmul(nqtot, &df[2 * dir + 1][0], 1, &inarray[0], 1, &tmp2[0], 1);
539  }
540  else
541  {
542  Vmath::Smul(nqtot, df[2 * dir][0], tmp0, 1, tmp0, 1);
543  Vmath::Smul(nqtot, df[2 * dir + 1][0], tmp1, 1, tmp1, 1);
544  Vmath::Smul(nqtot, df[2 * dir + 1][0], inarray, 1, tmp2, 1);
545  }
546  Vmath::Vadd(nqtot, tmp0, 1, tmp1, 1, tmp1, 1);
547 }
548 
550  const int dir, const Array<OneD, const NekDouble> &inarray,
551  Array<OneD, NekDouble> &outarray)
552 {
553  int nq = GetTotPoints();
555 
556  switch (dir)
557  {
558  case 0:
559  {
561  }
562  break;
563  case 1:
564  {
566  }
567  break;
568  case 2:
569  {
571  }
572  break;
573  default:
574  {
575  ASSERTL1(false, "input dir is out of range");
576  }
577  break;
578  }
579 
580  MatrixKey iprodmatkey(mtype, DetShapeType(), *this);
581  DNekScalMatSharedPtr iprodmat = m_matrixManager[iprodmatkey];
582 
583  Blas::Dgemv('N', m_ncoeffs, nq, iprodmat->Scale(),
584  (iprodmat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
585  inarray.get(), 1, 0.0, outarray.get(), 1);
586 }
587 
589  const Array<OneD, const NekDouble> &direction,
590  const Array<OneD, const NekDouble> &inarray,
591  Array<OneD, NekDouble> &outarray)
592 {
593  IProductWRTDirectionalDerivBase_SumFac(direction, inarray, outarray);
594 }
595 
596 /**
597  * @brief Directinoal Derivative in the modal space in the dir
598  * direction of varcoeffs.
599  */
601  const Array<OneD, const NekDouble> &direction,
602  const Array<OneD, const NekDouble> &inarray,
603  Array<OneD, NekDouble> &outarray)
604 {
605  int i;
606  int shapedim = 2;
607  int nquad0 = m_base[0]->GetNumPoints();
608  int nquad1 = m_base[1]->GetNumPoints();
609  int nqtot = nquad0 * nquad1;
610  int nmodes0 = m_base[0]->GetNumModes();
611  int wspsize = max(max(nqtot, m_ncoeffs), nquad1 * nmodes0);
612 
613  const Array<TwoD, const NekDouble> &df =
614  m_metricinfo->GetDerivFactors(GetPointsKeys());
615 
616  Array<OneD, NekDouble> tmp0(6 * wspsize);
617  Array<OneD, NekDouble> tmp1(tmp0 + wspsize);
618  Array<OneD, NekDouble> tmp2(tmp0 + 2 * wspsize);
619  Array<OneD, NekDouble> tmp3(tmp0 + 3 * wspsize);
620  Array<OneD, NekDouble> gfac0(tmp0 + 4 * wspsize);
621  Array<OneD, NekDouble> gfac1(tmp0 + 5 * wspsize);
622 
623  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
624  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
625 
626  // set up geometric factor: 2/(1-z1)
627  for (i = 0; i < nquad1; ++i)
628  {
629  gfac0[i] = 2.0 / (1 - z1[i]);
630  }
631  for (i = 0; i < nquad0; ++i)
632  {
633  gfac1[i] = 0.5 * (1 + z0[i]);
634  }
635  for (i = 0; i < nquad1; ++i)
636  {
637  Vmath::Smul(nquad0, gfac0[i], &inarray[0] + i * nquad0, 1,
638  &tmp0[0] + i * nquad0, 1);
639  }
640  for (i = 0; i < nquad1; ++i)
641  {
642  Vmath::Vmul(nquad0, &gfac1[0], 1, &tmp0[0] + i * nquad0, 1,
643  &tmp1[0] + i * nquad0, 1);
644  }
645 
646  // Compute gmat \cdot e^j
647  Array<OneD, Array<OneD, NekDouble>> dfdir(shapedim);
648  Expansion::ComputeGmatcdotMF(df, direction, dfdir);
649 
650  Vmath::Vmul(nqtot, &dfdir[0][0], 1, &tmp0[0], 1, &tmp0[0], 1);
651  Vmath::Vmul(nqtot, &dfdir[1][0], 1, &tmp1[0], 1, &tmp1[0], 1);
652  Vmath::Vmul(nqtot, &dfdir[1][0], 1, &inarray[0], 1, &tmp2[0], 1);
653 
654  Vmath::Vadd(nqtot, &tmp0[0], 1, &tmp1[0], 1, &tmp1[0], 1);
655 
656  MultiplyByQuadratureMetric(tmp1, tmp1);
657  MultiplyByQuadratureMetric(tmp2, tmp2);
658 
659  IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
660  tmp1, tmp3, tmp0);
661  IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
662  tmp2, outarray, tmp0);
663  Vmath::Vadd(m_ncoeffs, tmp3, 1, outarray, 1, outarray, 1);
664 }
665 
669  Array<OneD, NekDouble> &outarray)
670 {
671  int nq = m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints();
672  Array<OneD, NekDouble> Fn(nq);
673 
675  GetLeftAdjacentElementExp()->GetTraceNormal(
677 
678  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
679  {
680  Vmath::Vvtvvtp(nq, &normals[0][0], 1, &Fx[0], 1, &normals[1][0], 1,
681  &Fy[0], 1, &Fn[0], 1);
682  Vmath::Vvtvp(nq, &normals[2][0], 1, &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
683  }
684  else
685  {
686  Vmath::Svtsvtp(nq, normals[0][0], &Fx[0], 1, normals[1][0], &Fy[0], 1,
687  &Fn[0], 1);
688  Vmath::Svtvp(nq, normals[2][0], &Fz[0], 1, &Fn[0], 1, &Fn[0], 1);
689  }
690 
691  IProductWRTBase(Fn, outarray);
692 }
693 
695  const Array<OneD, const Array<OneD, NekDouble>> &Fvec,
696  Array<OneD, NekDouble> &outarray)
697 {
698  NormVectorIProductWRTBase(Fvec[0], Fvec[1], Fvec[2], outarray);
699 }
700 
702 {
703 
705  m_base[0]->GetBasisKey(), m_base[1]->GetBasisKey());
706 }
707 
709 {
711  m_base[0]->GetPointsKey());
713  m_base[1]->GetPointsKey());
714 
716  bkey1);
717 }
718 
720  Array<OneD, NekDouble> &coords)
721 {
722  int i;
723 
724  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[1] <= 1.0 && Lcoords[1] >= -1.0 &&
725  Lcoords[1] <= 1.0,
726  "Local coordinates are not in region [-1,1]");
727 
728  m_geom->FillGeom();
729 
730  for (i = 0; i < m_geom->GetCoordim(); ++i)
731  {
732  coords[i] = m_geom->GetCoord(i, Lcoords);
733  }
734 }
735 
737  Array<OneD, NekDouble> &coords_1,
738  Array<OneD, NekDouble> &coords_2)
739 {
740  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
741 }
742 
743 /**
744  * Given the local cartesian coordinate \a Lcoord evaluate the
745  * value of physvals at this point by calling through to the
746  * StdExpansion method
747  */
749  const Array<OneD, const NekDouble> &Lcoord,
750  const Array<OneD, const NekDouble> &physvals)
751 {
752  // Evaluate point in local (eta) coordinates.
753  return StdTriExp::v_PhysEvaluate(Lcoord, physvals);
754 }
755 
757  const Array<OneD, const NekDouble> &physvals)
758 {
760 
761  ASSERTL0(m_geom, "m_geom not defined");
762  m_geom->GetLocCoords(coord, Lcoord);
763 
764  return StdTriExp::v_PhysEvaluate(Lcoord, physvals);
765 }
766 
768  const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp,
769  const Array<OneD, const NekDouble> &inarray,
771 {
772  int nquad0 = m_base[0]->GetNumPoints();
773  int nquad1 = m_base[1]->GetNumPoints();
774  int nt = 0;
775  // Extract in Cartesian direction because we have to deal with
776  // e.g. Gauss-Radau points.
777  switch (edge)
778  {
779  case 0:
780  Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
781  nt = nquad0;
782  break;
783  case 1:
784  Vmath::Vcopy(nquad1, &(inarray[0]) + (nquad0 - 1), nquad0,
785  &(outarray[0]), 1);
786  nt = nquad1;
787  break;
788  case 2:
789  Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
790  nt = nquad1;
791  break;
792  default:
793  ASSERTL0(false, "edge value (< 3) is out of range");
794  break;
795  }
796 
797  ASSERTL1(EdgeExp->GetBasis(0)->GetPointsType() ==
799  "Edge expansion should be GLL");
800 
801  // Interpolate if required
802  if (m_base[edge ? 1 : 0]->GetPointsKey() !=
803  EdgeExp->GetBasis(0)->GetPointsKey())
804  {
805  Array<OneD, NekDouble> outtmp(max(nquad0, nquad1));
806 
807  Vmath::Vcopy(nt, outarray, 1, outtmp, 1);
808 
809  LibUtilities::Interp1D(m_base[edge ? 1 : 0]->GetPointsKey(), outtmp,
810  EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
811  }
812 
813  if (orient == StdRegions::eNoOrientation)
814  {
815  orient = GetTraceOrient(edge);
816  }
817 
818  // Reverse data if necessary
819  if (orient == StdRegions::eBackwards)
820  {
821  Vmath::Reverse(EdgeExp->GetNumPoints(0), &outarray[0], 1, &outarray[0],
822  1);
823  }
824 }
825 
826 void TriExp::v_GetEdgeInterpVals(const int edge,
827  const Array<OneD, const NekDouble> &inarray,
828  Array<OneD, NekDouble> &outarray)
829 {
830  boost::ignore_unused(edge, inarray, outarray);
831  ASSERTL0(false, "Routine not implemented for triangular elements");
832 }
833 
834 void TriExp::v_GetTraceQFactors(const int edge,
835  Array<OneD, NekDouble> &outarray)
836 {
837  boost::ignore_unused(edge, outarray);
838  ASSERTL0(false, "Routine not implemented for triangular elements");
839 }
840 
841 void TriExp::v_GetTracePhysMap(const int edge, Array<OneD, int> &outarray)
842 {
843  int nquad0 = m_base[0]->GetNumPoints();
844  int nquad1 = m_base[1]->GetNumPoints();
845 
846  // Get points in Cartesian orientation
847  switch (edge)
848  {
849  case 0:
850  outarray = Array<OneD, int>(nquad0);
851  for (int i = 0; i < nquad0; ++i)
852  {
853  outarray[i] = i;
854  }
855  break;
856  case 1:
857  outarray = Array<OneD, int>(nquad1);
858  for (int i = 0; i < nquad1; ++i)
859  {
860  outarray[i] = (nquad0 - 1) + i * nquad0;
861  }
862  break;
863  case 2:
864  outarray = Array<OneD, int>(nquad1);
865  for (int i = 0; i < nquad1; ++i)
866  {
867  outarray[i] = i * nquad0;
868  }
869  break;
870  default:
871  ASSERTL0(false, "edge value (< 3) is out of range");
872  break;
873  }
874 }
875 
876 void TriExp::v_ComputeTraceNormal(const int edge)
877 {
878  int i;
879  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
880  GetGeom()->GetMetricInfo();
881 
883  for (i = 0; i < ptsKeys.size(); ++i)
884  {
885  // Need at least 2 points for computing normals
886  if (ptsKeys[i].GetNumPoints() == 1)
887  {
888  LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
889  ptsKeys[i] = pKey;
890  }
891  }
892 
893  const SpatialDomains::GeomType type = geomFactors->GetGtype();
894  const Array<TwoD, const NekDouble> &df =
895  geomFactors->GetDerivFactors(ptsKeys);
896  const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
897  int nqe = m_base[0]->GetNumPoints();
898  int dim = GetCoordim();
899 
902  for (i = 0; i < dim; ++i)
903  {
904  normal[i] = Array<OneD, NekDouble>(nqe);
905  }
906 
907  size_t nqb = nqe;
908  size_t nbnd = edge;
911 
912  // Regular geometry case
913  if ((type == SpatialDomains::eRegular) ||
915  {
916  NekDouble fac;
917  // Set up normals
918  switch (edge)
919  {
920  case 0:
921  for (i = 0; i < GetCoordim(); ++i)
922  {
923  Vmath::Fill(nqe, -df[2 * i + 1][0], normal[i], 1);
924  }
925  break;
926  case 1:
927  for (i = 0; i < GetCoordim(); ++i)
928  {
929  Vmath::Fill(nqe, df[2 * i + 1][0] + df[2 * i][0], normal[i],
930  1);
931  }
932  break;
933  case 2:
934  for (i = 0; i < GetCoordim(); ++i)
935  {
936  Vmath::Fill(nqe, -df[2 * i][0], normal[i], 1);
937  }
938  break;
939  default:
940  ASSERTL0(false, "Edge is out of range (edge < 3)");
941  }
942 
943  // normalise
944  fac = 0.0;
945  for (i = 0; i < GetCoordim(); ++i)
946  {
947  fac += normal[i][0] * normal[i][0];
948  }
949  fac = 1.0 / sqrt(fac);
950 
951  Vmath::Fill(nqb, fac, length, 1);
952 
953  for (i = 0; i < GetCoordim(); ++i)
954  {
955  Vmath::Smul(nqe, fac, normal[i], 1, normal[i], 1);
956  }
957  }
958  else // Set up deformed normals
959  {
960  int j;
961 
962  int nquad0 = ptsKeys[0].GetNumPoints();
963  int nquad1 = ptsKeys[1].GetNumPoints();
964 
965  LibUtilities::PointsKey from_key;
966 
967  Array<OneD, NekDouble> normals(GetCoordim() * max(nquad0, nquad1), 0.0);
968  Array<OneD, NekDouble> edgejac(GetCoordim() * max(nquad0, nquad1), 0.0);
969 
970  // Extract Jacobian along edges and recover local
971  // derivates (dx/dr) for polynomial interpolation by
972  // multiplying m_gmat by jacobian
973  switch (edge)
974  {
975  case 0:
976  for (j = 0; j < nquad0; ++j)
977  {
978  edgejac[j] = jac[j];
979  for (i = 0; i < GetCoordim(); ++i)
980  {
981  normals[i * nquad0 + j] =
982  -df[2 * i + 1][j] * edgejac[j];
983  }
984  }
985  from_key = ptsKeys[0];
986  break;
987  case 1:
988  for (j = 0; j < nquad1; ++j)
989  {
990  edgejac[j] = jac[nquad0 * j + nquad0 - 1];
991  for (i = 0; i < GetCoordim(); ++i)
992  {
993  normals[i * nquad1 + j] =
994  (df[2 * i][nquad0 * j + nquad0 - 1] +
995  df[2 * i + 1][nquad0 * j + nquad0 - 1]) *
996  edgejac[j];
997  }
998  }
999  from_key = ptsKeys[1];
1000  break;
1001  case 2:
1002  for (j = 0; j < nquad1; ++j)
1003  {
1004  edgejac[j] = jac[nquad0 * j];
1005  for (i = 0; i < GetCoordim(); ++i)
1006  {
1007  normals[i * nquad1 + j] =
1008  -df[2 * i][nquad0 * j] * edgejac[j];
1009  }
1010  }
1011  from_key = ptsKeys[1];
1012  break;
1013  default:
1014  ASSERTL0(false, "edge is out of range (edge < 3)");
1015  }
1016 
1017  int nq = from_key.GetNumPoints();
1018  Array<OneD, NekDouble> work(nqe, 0.0);
1019 
1020  // interpolate Jacobian and invert
1021  LibUtilities::Interp1D(from_key, jac, m_base[0]->GetPointsKey(), work);
1022  Vmath::Sdiv(nqe, 1.0, &work[0], 1, &work[0], 1);
1023 
1024  // interpolate
1025  for (i = 0; i < GetCoordim(); ++i)
1026  {
1027  LibUtilities::Interp1D(from_key, &normals[i * nq],
1028  m_base[0]->GetPointsKey(), &normal[i][0]);
1029  Vmath::Vmul(nqe, work, 1, normal[i], 1, normal[i], 1);
1030  }
1031 
1032  // normalise normal vectors
1033  Vmath::Zero(nqe, work, 1);
1034  for (i = 0; i < GetCoordim(); ++i)
1035  {
1036  Vmath::Vvtvp(nqe, normal[i], 1, normal[i], 1, work, 1, work, 1);
1037  }
1038 
1039  Vmath::Vsqrt(nqe, work, 1, work, 1);
1040  Vmath::Sdiv(nqe, 1.0, work, 1, work, 1);
1041 
1042  Vmath::Vcopy(nqb, work, 1, length, 1);
1043 
1044  for (i = 0; i < GetCoordim(); ++i)
1045  {
1046  Vmath::Vmul(nqe, normal[i], 1, work, 1, normal[i], 1);
1047  }
1048  }
1049 
1050  if (GetGeom()->GetEorient(edge) == StdRegions::eBackwards)
1051  {
1052  for (i = 0; i < GetCoordim(); ++i)
1053  {
1054  if (geomFactors->GetGtype() == SpatialDomains::eDeformed)
1055  {
1056  Vmath::Reverse(nqe, normal[i], 1, normal[i], 1);
1057  }
1058  }
1059  }
1060 }
1061 
1063 {
1064  return m_geom->GetCoordim();
1065 }
1066 
1068  const NekDouble *data, const std::vector<unsigned int> &nummodes,
1069  const int mode_offset, NekDouble *coeffs,
1070  std::vector<LibUtilities::BasisType> &fromType)
1071 {
1072  boost::ignore_unused(fromType);
1073 
1074  int data_order0 = nummodes[mode_offset];
1075  int fillorder0 = min(m_base[0]->GetNumModes(), data_order0);
1076  int data_order1 = nummodes[mode_offset + 1];
1077  int order1 = m_base[1]->GetNumModes();
1078  int fillorder1 = min(order1, data_order1);
1079 
1080  switch (m_base[0]->GetBasisType())
1081  {
1084  {
1085  int i;
1086  int cnt = 0;
1087  int cnt1 = 0;
1088 
1091  "Extraction routine not set up for this basis");
1092 
1093  Vmath::Zero(m_ncoeffs, coeffs, 1);
1094  for (i = 0; i < fillorder0; ++i)
1095  {
1096  Vmath::Vcopy(fillorder1 - i, &data[cnt], 1, &coeffs[cnt1], 1);
1097  cnt += data_order1 - i;
1098  cnt1 += order1 - i;
1099  }
1100  }
1101  break;
1102  default:
1103  ASSERTL0(false, "basis is either not set up or not hierarchicial");
1104  }
1105 }
1106 
1108 {
1109  return GetGeom2D()->GetEorient(edge);
1110 }
1111 
1113 {
1114  ASSERTL1(dir >= 0 && dir <= 1, "input dir is out of range");
1115  return m_base[dir];
1116 }
1117 
1118 int TriExp::v_GetNumPoints(const int dir) const
1119 {
1120  return GetNumPoints(dir);
1121 }
1122 
1124 {
1125  DNekMatSharedPtr returnval;
1126  switch (mkey.GetMatrixType())
1127  {
1135  returnval = Expansion2D::v_GenMatrix(mkey);
1136  break;
1137  default:
1138  returnval = StdTriExp::v_GenMatrix(mkey);
1139  break;
1140  }
1141 
1142  return returnval;
1143 }
1144 
1146 {
1147  LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1148  LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1151 
1152  return tmp->GetStdMatrix(mkey);
1153 }
1154 
1156 {
1157  return m_matrixManager[mkey];
1158 }
1159 
1161 {
1162  return m_staticCondMatrixManager[mkey];
1163 }
1164 
1166 {
1167  m_staticCondMatrixManager.DeleteObject(mkey);
1168 }
1169 
1171  Array<OneD, NekDouble> &outarray,
1172  const StdRegions::StdMatrixKey &mkey)
1173 {
1174  StdExpansion::MassMatrixOp_MatFree(inarray, outarray, mkey);
1175 }
1176 
1178  Array<OneD, NekDouble> &outarray,
1179  const StdRegions::StdMatrixKey &mkey)
1180 {
1181  TriExp::LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
1182 }
1183 
1184 void TriExp::v_LaplacianMatrixOp(const int k1, const int k2,
1185  const Array<OneD, const NekDouble> &inarray,
1186  Array<OneD, NekDouble> &outarray,
1187  const StdRegions::StdMatrixKey &mkey)
1188 {
1189  StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
1190 }
1191 
1193  const Array<OneD, const NekDouble> &inarray,
1194  Array<OneD, NekDouble> &outarray,
1195  const StdRegions::StdMatrixKey &mkey)
1196 {
1197  StdExpansion::WeakDerivMatrixOp_MatFree(i, inarray, outarray, mkey);
1198 }
1199 
1201  const Array<OneD, const NekDouble> &inarray,
1202  Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1203 {
1204  StdExpansion::WeakDirectionalDerivMatrixOp_MatFree(inarray, outarray, mkey);
1205 }
1206 
1208  const Array<OneD, const NekDouble> &inarray,
1209  Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1210 {
1211  StdExpansion::MassLevelCurvatureMatrixOp_MatFree(inarray, outarray, mkey);
1212 }
1213 
1215  Array<OneD, NekDouble> &outarray,
1216  const StdRegions::StdMatrixKey &mkey)
1217 {
1218  TriExp::HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
1219 }
1220 
1222  const Array<OneD, const NekDouble> &inarray,
1223  Array<OneD, NekDouble> &outarray, const StdRegions::StdMatrixKey &mkey)
1224 {
1225  DNekScalMatSharedPtr mat = GetLocMatrix(mkey);
1226 
1227  if (inarray.get() == outarray.get())
1228  {
1230  Vmath::Vcopy(m_ncoeffs, inarray.get(), 1, tmp.get(), 1);
1231 
1232  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, mat->Scale(),
1233  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
1234  tmp.get(), 1, 0.0, outarray.get(), 1);
1235  }
1236  else
1237  {
1238  Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, mat->Scale(),
1239  (mat->GetOwnedMatrix())->GetPtr().get(), m_ncoeffs,
1240  inarray.get(), 1, 0.0, outarray.get(), 1);
1241  }
1242 }
1243 
1245  const Array<OneD, const NekDouble> &inarray,
1247 {
1248  if (m_metrics.count(eMetricLaplacian00) == 0)
1249  {
1251  }
1252 
1253  int nquad0 = m_base[0]->GetNumPoints();
1254  int nquad1 = m_base[1]->GetNumPoints();
1255  int nqtot = nquad0 * nquad1;
1256  int nmodes0 = m_base[0]->GetNumModes();
1257  int nmodes1 = m_base[1]->GetNumModes();
1258  int wspsize =
1259  max(max(max(nqtot, m_ncoeffs), nquad1 * nmodes0), nquad0 * nmodes1);
1260 
1261  ASSERTL1(wsp.size() >= 3 * wspsize, "Workspace is of insufficient size.");
1262 
1263  const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1264  const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1265  const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1266  const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1267  const Array<OneD, const NekDouble> &metric00 =
1269  const Array<OneD, const NekDouble> &metric01 =
1271  const Array<OneD, const NekDouble> &metric11 =
1273 
1274  // Allocate temporary storage
1275  Array<OneD, NekDouble> wsp0(wsp);
1276  Array<OneD, NekDouble> wsp1(wsp + wspsize);
1277  Array<OneD, NekDouble> wsp2(wsp + 2 * wspsize);
1278 
1279  StdExpansion2D::PhysTensorDeriv(inarray, wsp1, wsp2);
1280 
1281  // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1282  // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1283  // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1284  // especially for this purpose
1285  Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp1[0], 1, &metric01[0], 1,
1286  &wsp2[0], 1, &wsp0[0], 1);
1287  Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp1[0], 1, &metric11[0], 1,
1288  &wsp2[0], 1, &wsp2[0], 1);
1289 
1290  // outarray = m = (D_xi1 * B)^T * k
1291  // wsp1 = n = (D_xi2 * B)^T * l
1292  IProductWRTBase_SumFacKernel(dbase0, base1, wsp0, outarray, wsp1);
1293  IProductWRTBase_SumFacKernel(base0, dbase1, wsp2, wsp1, wsp0);
1294 
1295  // outarray = outarray + wsp1
1296  // = L * u_hat
1297  Vmath::Vadd(m_ncoeffs, wsp1.get(), 1, outarray.get(), 1, outarray.get(), 1);
1298 }
1299 
1301 {
1302  if (m_metrics.count(eMetricQuadrature) == 0)
1303  {
1305  }
1306 
1307  unsigned int i, j;
1308  const SpatialDomains::GeomType type = m_metricinfo->GetGtype();
1309  const unsigned int nqtot = GetTotPoints();
1310  const unsigned int dim = 2;
1311  const MetricType m[3][3] = {
1315 
1316  Array<OneD, NekDouble> dEta_dXi[2] = {Array<OneD, NekDouble>(nqtot, 1.0),
1317  Array<OneD, NekDouble>(nqtot, 1.0)};
1318 
1319  for (i = 0; i < dim; ++i)
1320  {
1321  for (j = i; j < dim; ++j)
1322  {
1323  m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1324  }
1325  }
1326 
1327  const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
1328  const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1329  const unsigned int nquad0 = m_base[0]->GetNumPoints();
1330  const unsigned int nquad1 = m_base[1]->GetNumPoints();
1331  const Array<TwoD, const NekDouble> &df =
1332  m_metricinfo->GetDerivFactors(GetPointsKeys());
1333 
1334  for (i = 0; i < nquad1; i++)
1335  {
1336  Blas::Dscal(nquad0, 2.0 / (1 - z1[i]), &dEta_dXi[0][0] + i * nquad0, 1);
1337  Blas::Dscal(nquad0, 2.0 / (1 - z1[i]), &dEta_dXi[1][0] + i * nquad0, 1);
1338  }
1339  for (i = 0; i < nquad0; i++)
1340  {
1341  Blas::Dscal(nquad1, 0.5 * (1 + z0[i]), &dEta_dXi[1][0] + i, nquad0);
1342  }
1343 
1344  Array<OneD, NekDouble> tmp(nqtot);
1345  if ((type == SpatialDomains::eRegular ||
1347  {
1348  Vmath::Smul(nqtot, df[0][0], &dEta_dXi[0][0], 1, &tmp[0], 1);
1349  Vmath::Svtvp(nqtot, df[1][0], &dEta_dXi[1][0], 1, &tmp[0], 1, &tmp[0],
1350  1);
1351 
1352  Vmath::Vmul(nqtot, &tmp[0], 1, &tmp[0], 1,
1353  &m_metrics[eMetricLaplacian00][0], 1);
1354  Vmath::Smul(nqtot, df[1][0], &tmp[0], 1,
1355  &m_metrics[eMetricLaplacian01][0], 1);
1356 
1357  Vmath::Smul(nqtot, df[2][0], &dEta_dXi[0][0], 1, &tmp[0], 1);
1358  Vmath::Svtvp(nqtot, df[3][0], &dEta_dXi[1][0], 1, &tmp[0], 1, &tmp[0],
1359  1);
1360 
1361  Vmath::Vvtvp(nqtot, &tmp[0], 1, &tmp[0], 1,
1362  &m_metrics[eMetricLaplacian00][0], 1,
1363  &m_metrics[eMetricLaplacian00][0], 1);
1364  Vmath::Svtvp(nqtot, df[3][0], &tmp[0], 1,
1365  &m_metrics[eMetricLaplacian01][0], 1,
1366  &m_metrics[eMetricLaplacian01][0], 1);
1367 
1368  if (GetCoordim() == 3)
1369  {
1370  Vmath::Smul(nqtot, df[4][0], &dEta_dXi[0][0], 1, &tmp[0], 1);
1371  Vmath::Svtvp(nqtot, df[5][0], &dEta_dXi[1][0], 1, &tmp[0], 1,
1372  &tmp[0], 1);
1373 
1374  Vmath::Vvtvp(nqtot, &tmp[0], 1, &tmp[0], 1,
1375  &m_metrics[eMetricLaplacian00][0], 1,
1376  &m_metrics[eMetricLaplacian00][0], 1);
1377  Vmath::Svtvp(nqtot, df[5][0], &tmp[0], 1,
1378  &m_metrics[eMetricLaplacian01][0], 1,
1379  &m_metrics[eMetricLaplacian01][0], 1);
1380  }
1381 
1382  NekDouble g2 = df[1][0] * df[1][0] + df[3][0] * df[3][0];
1383  if (GetCoordim() == 3)
1384  {
1385  g2 += df[5][0] * df[5][0];
1386  }
1387  Vmath::Fill(nqtot, g2, &m_metrics[eMetricLaplacian11][0], 1);
1388  }
1389  else
1390  {
1391 
1392  Vmath::Vmul(nqtot, &df[0][0], 1, &dEta_dXi[0][0], 1, &tmp[0], 1);
1393  Vmath::Vvtvp(nqtot, &df[1][0], 1, &dEta_dXi[1][0], 1, &tmp[0], 1,
1394  &tmp[0], 1);
1395 
1396  Vmath::Vmul(nqtot, &tmp[0], 1, &tmp[0], 1,
1397  &m_metrics[eMetricLaplacian00][0], 1);
1398  Vmath::Vmul(nqtot, &df[1][0], 1, &tmp[0], 1,
1399  &m_metrics[eMetricLaplacian01][0], 1);
1400  Vmath::Vmul(nqtot, &df[1][0], 1, &df[1][0], 1,
1401  &m_metrics[eMetricLaplacian11][0], 1);
1402 
1403  Vmath::Vmul(nqtot, &df[2][0], 1, &dEta_dXi[0][0], 1, &tmp[0], 1);
1404  Vmath::Vvtvp(nqtot, &df[3][0], 1, &dEta_dXi[1][0], 1, &tmp[0], 1,
1405  &tmp[0], 1);
1406 
1407  Vmath::Vvtvp(nqtot, &tmp[0], 1, &tmp[0], 1,
1408  &m_metrics[eMetricLaplacian00][0], 1,
1409  &m_metrics[eMetricLaplacian00][0], 1);
1410  Vmath::Vvtvp(nqtot, &df[3][0], 1, &tmp[0], 1,
1411  &m_metrics[eMetricLaplacian01][0], 1,
1412  &m_metrics[eMetricLaplacian01][0], 1);
1413  Vmath::Vvtvp(nqtot, &df[3][0], 1, &df[3][0], 1,
1414  &m_metrics[eMetricLaplacian11][0], 1,
1415  &m_metrics[eMetricLaplacian11][0], 1);
1416 
1417  if (GetCoordim() == 3)
1418  {
1419  Vmath::Vmul(nqtot, &df[4][0], 1, &dEta_dXi[0][0], 1, &tmp[0], 1);
1420  Vmath::Vvtvp(nqtot, &df[5][0], 1, &dEta_dXi[1][0], 1, &tmp[0], 1,
1421  &tmp[0], 1);
1422 
1423  Vmath::Vvtvp(nqtot, &tmp[0], 1, &tmp[0], 1,
1424  &m_metrics[eMetricLaplacian00][0], 1,
1425  &m_metrics[eMetricLaplacian00][0], 1);
1426  Vmath::Vvtvp(nqtot, &df[5][0], 1, &tmp[0], 1,
1427  &m_metrics[eMetricLaplacian01][0], 1,
1428  &m_metrics[eMetricLaplacian01][0], 1);
1429  Vmath::Vvtvp(nqtot, &df[5][0], 1, &df[5][0], 1,
1430  &m_metrics[eMetricLaplacian11][0], 1,
1431  &m_metrics[eMetricLaplacian11][0], 1);
1432  }
1433  }
1434 
1435  for (unsigned int i = 0; i < dim; ++i)
1436  {
1437  for (unsigned int j = i; j < dim; ++j)
1438  {
1439  MultiplyByQuadratureMetric(m_metrics[m[i][j]], m_metrics[m[i][j]]);
1440  }
1441  }
1442 }
1443 
1444 /**
1445  * Function is used to compute exactly the advective numerical flux on
1446  * theinterface of two elements with different expansions, hence an
1447  * appropriate number of Gauss points has to be used. The number of
1448  * Gauss points has to be equal to the number used by the highest
1449  * polynomial degree of the two adjacent elements. Furthermore, this
1450  * function is used to compute the sensor value in each element.
1451  *
1452  * @param numMin Is the reduced polynomial order
1453  * @param inarray Input array of coefficients
1454  * @param dumpVar Output array of reduced coefficients.
1455  */
1457  const Array<OneD, const NekDouble> &inarray,
1458  Array<OneD, NekDouble> &outarray)
1459 {
1460  int n_coeffs = inarray.size();
1461  int nquad0 = m_base[0]->GetNumPoints();
1462  int nquad1 = m_base[1]->GetNumPoints();
1463  int nqtot = nquad0 * nquad1;
1464  int nmodes0 = m_base[0]->GetNumModes();
1465  int nmodes1 = m_base[1]->GetNumModes();
1466  int numMin2 = nmodes0, i;
1467 
1468  Array<OneD, NekDouble> coeff(n_coeffs, 0.0);
1469  Array<OneD, NekDouble> phys_tmp(nqtot, 0.0);
1470  Array<OneD, NekDouble> tmp, tmp2;
1471 
1472  const LibUtilities::PointsKey Pkey0 = m_base[0]->GetPointsKey();
1473  const LibUtilities::PointsKey Pkey1 = m_base[1]->GetPointsKey();
1474 
1475  LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
1476  LibUtilities::BasisKey b1(m_base[1]->GetBasisType(), nmodes1, Pkey1);
1477  LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
1478  LibUtilities::BasisKey bortho1(LibUtilities::eOrtho_B, nmodes1, Pkey1);
1479 
1480  // Check if it is also possible to use the same InterCoeff routine
1481  // which is also used for Quadrilateral and Hexagonal shaped
1482  // elements
1483 
1484  // For now, set up the used basis on the standard element to
1485  // calculate the phys values, set up the orthogonal basis to do a
1486  // forward transform, to obtain the coefficients in orthogonal
1487  // coefficient space
1488  StdRegions::StdTriExpSharedPtr m_OrthoTriExp;
1490 
1493  bortho0, bortho1);
1494 
1495  m_TriExp->BwdTrans(inarray, phys_tmp);
1496  m_OrthoTriExp->FwdTrans(phys_tmp, coeff);
1497 
1498  for (i = 0; i < n_coeffs; i++)
1499  {
1500  if (i == numMin)
1501  {
1502  coeff[i] = 0.0;
1503  numMin += numMin2 - 1;
1504  numMin2 -= 1.0;
1505  }
1506  }
1507 
1508  m_OrthoTriExp->BwdTrans(coeff, phys_tmp);
1509  m_TriExp->FwdTrans(phys_tmp, outarray);
1510 }
1511 
1513  const StdRegions::StdMatrixKey &mkey)
1514 {
1515  int nq = GetTotPoints();
1516 
1517  // Calculate sqrt of the Jacobian
1519  Array<OneD, NekDouble> sqrt_jac(nq);
1520  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1521  {
1522  Vmath::Vsqrt(nq, jac, 1, sqrt_jac, 1);
1523  }
1524  else
1525  {
1526  Vmath::Fill(nq, sqrt(jac[0]), sqrt_jac, 1);
1527  }
1528 
1529  // Multiply array by sqrt(Jac)
1530  Vmath::Vmul(nq, sqrt_jac, 1, array, 1, array, 1);
1531 
1532  // Apply std region filter
1533  StdTriExp::v_SVVLaplacianFilter(array, mkey);
1534 
1535  // Divide by sqrt(Jac)
1536  Vmath::Vdiv(nq, array, 1, sqrt_jac, 1, array, 1);
1537 }
1538 
1539 /** @brief: This method gets all of the factors which are
1540  required as part of the Gradient Jump Penalty
1541  stabilisation and involves the product of the normal and
1542  geometric factors along the element trace.
1543 */
1545  Array<OneD, Array<OneD, NekDouble>> &d0factors,
1546  Array<OneD, Array<OneD, NekDouble>> &d1factors,
1547  Array<OneD, Array<OneD, NekDouble>> &d2factors)
1548 {
1549  boost::ignore_unused(d2factors); // for 3D shapes
1550  int nquad0 = GetNumPoints(0);
1551  int nquad1 = GetNumPoints(1);
1552 
1553  const Array<TwoD, const NekDouble> &df =
1554  m_metricinfo->GetDerivFactors(GetPointsKeys());
1555 
1556  if (d0factors.size() != 3)
1557  {
1558  d0factors = Array<OneD, Array<OneD, NekDouble>>(3);
1559  d1factors = Array<OneD, Array<OneD, NekDouble>>(3);
1560  }
1561 
1562  if (d0factors[0].size() != nquad0)
1563  {
1564  d0factors[0] = Array<OneD, NekDouble>(nquad0);
1565  d1factors[0] = Array<OneD, NekDouble>(nquad0);
1566  }
1567 
1568  if (d0factors[1].size() != nquad1)
1569  {
1570  d0factors[1] = Array<OneD, NekDouble>(nquad1);
1571  d0factors[2] = Array<OneD, NekDouble>(nquad1);
1572  d1factors[1] = Array<OneD, NekDouble>(nquad1);
1573  d1factors[2] = Array<OneD, NekDouble>(nquad1);
1574  }
1575 
1576  // Outwards normals
1577  const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1578  GetTraceNormal(0);
1579  const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1580  GetTraceNormal(1);
1581  const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1582  GetTraceNormal(2);
1583 
1584  int ncoords = normal_0.size();
1585 
1586  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1587  {
1588 
1589  // d xi_2/dx n_x
1590  for (int i = 0; i < nquad0; ++i)
1591  {
1592  d1factors[0][i] = df[1][i] * normal_0[0][i];
1593  }
1594 
1595  // d xi_1/dx n_x
1596  for (int i = 0; i < nquad1; ++i)
1597  {
1598  d0factors[1][i] = df[0][(i + 1) * nquad0 - 1] * normal_1[0][i];
1599  d0factors[2][i] = df[0][i * nquad0] * normal_2[0][i];
1600  }
1601 
1602  for (int n = 1; n < ncoords; ++n)
1603  {
1604  // d xi_2/dy n_y
1605  // needs checking for 3D coords
1606  for (int i = 0; i < nquad0; ++i)
1607  {
1608  d1factors[0][i] += df[2 * n + 1][i] * normal_0[n][i];
1609  }
1610 
1611  // d xi_1/dy n_y
1612  // needs checking for 3D coords
1613  for (int i = 0; i < nquad1; ++i)
1614  {
1615  d0factors[1][i] +=
1616  df[2 * n][(i + 1) * nquad0 - 1] * normal_1[n][i];
1617  d0factors[2][i] += df[2 * n][i * nquad0] * normal_2[n][i];
1618  }
1619  }
1620 
1621  // d0 factors
1622  // d xi_1/dx n_x
1623  for (int i = 0; i < nquad0; ++i)
1624  {
1625  d0factors[0][i] = df[0][i] * normal_0[0][i];
1626  }
1627 
1628  // d xi_2/dx n_x
1629  for (int i = 0; i < nquad1; ++i)
1630  {
1631  d1factors[1][i] = df[1][(i + 1) * nquad0 - 1] * normal_1[0][i];
1632  d1factors[2][i] = df[1][i * nquad0] * normal_2[0][i];
1633  }
1634 
1635  for (int n = 1; n < ncoords; ++n)
1636  {
1637  // d xi_1/dy n_y
1638  // needs checking for 3D coords
1639  for (int i = 0; i < nquad0; ++i)
1640  {
1641  d0factors[0][i] += df[2 * n][i] * normal_0[n][i];
1642  }
1643 
1644  // d xi_2/dy n_y
1645  // needs checking for 3D coords
1646  for (int i = 0; i < nquad1; ++i)
1647  {
1648  d1factors[1][i] +=
1649  df[2 * n + 1][(i + 1) * nquad0 - 1] * normal_1[n][i];
1650  d1factors[2][i] += df[2 * n + 1][i * nquad0] * normal_2[n][i];
1651  }
1652  }
1653  }
1654  else
1655  {
1656  // d xi_2/dx n_x
1657  for (int i = 0; i < nquad0; ++i)
1658  {
1659  d1factors[0][i] = df[1][0] * normal_0[0][i];
1660  }
1661 
1662  // d xi_1/dx n_x
1663  for (int i = 0; i < nquad1; ++i)
1664  {
1665  d0factors[1][i] = df[0][0] * normal_1[0][i];
1666  d0factors[2][i] = df[0][0] * normal_2[0][i];
1667  }
1668 
1669  for (int n = 1; n < ncoords; ++n)
1670  {
1671  // d xi_2/dy n_y
1672  // needs checking for 3D coords
1673  for (int i = 0; i < nquad0; ++i)
1674  {
1675  d1factors[0][i] += df[2 * n + 1][0] * normal_0[n][i];
1676  }
1677 
1678  // d xi_1/dy n_y
1679  // needs checking for 3D coords
1680  for (int i = 0; i < nquad1; ++i)
1681  {
1682  d0factors[1][i] += df[2 * n][0] * normal_1[n][i];
1683  d0factors[2][i] += df[2 * n][0] * normal_2[n][i];
1684  }
1685  }
1686 
1687  // d1factors
1688  // d xi_1/dx n_x
1689  for (int i = 0; i < nquad0; ++i)
1690  {
1691  d0factors[0][i] = df[0][0] * normal_0[0][i];
1692  }
1693 
1694  // d xi_2/dx n_x
1695  for (int i = 0; i < nquad1; ++i)
1696  {
1697  d1factors[1][i] = df[1][0] * normal_1[0][i];
1698  d1factors[2][i] = df[1][0] * normal_2[0][i];
1699  }
1700 
1701  for (int n = 1; n < ncoords; ++n)
1702  {
1703  // d xi_1/dy n_y
1704  // needs checking for 3D coords
1705  for (int i = 0; i < nquad0; ++i)
1706  {
1707  d0factors[0][i] += df[2 * n][0] * normal_0[n][i];
1708  }
1709 
1710  // d xi_2/dy n_y
1711  // needs checking for 3D coords
1712  for (int i = 0; i < nquad1; ++i)
1713  {
1714  d1factors[1][i] += df[2 * n + 1][0] * normal_1[n][i];
1715  d1factors[2][i] += df[2 * n + 1][0] * normal_2[n][i];
1716  }
1717  }
1718  }
1719 }
1720 } // namespace LocalRegions
1721 } // 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 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
unsigned int GetNumPoints() const
Definition: Points.h:104
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
Definition: Expansion2D.h:180
std::map< int, NormalVector > m_traceNormals
Definition: Expansion.h:275
std::map< int, Array< OneD, NekDouble > > m_elmtBndNormDirElmtLen
the element length in each element boundary(Vertex, edge or face) normal direction calculated based o...
Definition: Expansion.h:285
SpatialDomains::GeometrySharedPtr GetGeom() const
Definition: Expansion.cpp:166
SpatialDomains::GeometrySharedPtr m_geom
Definition: Expansion.h:272
ExpansionSharedPtr GetLeftAdjacentElementExp() const
Definition: Expansion.h:433
void ComputeGmatcdotMF(const Array< TwoD, const NekDouble > &df, const Array< OneD, const NekDouble > &direction, Array< OneD, Array< OneD, NekDouble >> &dfdir)
Definition: Expansion.cpp:597
int GetLeftAdjacentElementTrace() const
Definition: Expansion.h:446
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:273
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: Expansion.cpp:524
StdRegions::Orientation GetTraceOrient(int trace)
Definition: Expansion.h:167
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:88
const NormalVector & GetTraceNormal(const int id)
Definition: Expansion.cpp:250
virtual void v_ComputeTraceNormal(const int edge)
Definition: TriExp.cpp:876
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: TriExp.h:249
virtual void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &direction, Array< OneD, NekDouble > &out)
Physical derivative along a direction vector.
Definition: TriExp.cpp:192
virtual void v_NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, const Array< OneD, const NekDouble > &Fz, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:666
virtual void v_GetTracePhysMap(const int edge, Array< OneD, int > &outarray)
Definition: TriExp.cpp:841
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords)
Definition: TriExp.cpp:719
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1177
virtual void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType)
Definition: TriExp.cpp:1067
virtual void v_WeakDerivMatrixOp(const int i, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1192
virtual void v_MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1170
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1155
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:273
virtual void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:1456
virtual const LibUtilities::BasisSharedPtr & v_GetBasis(int dir) const
Definition: TriExp.cpp:1112
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:406
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1160
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: TriExp.cpp:102
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: TriExp.h:251
virtual void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble >> &outarray)
Definition: TriExp.cpp:484
virtual int v_GetCoordim()
Definition: TriExp.cpp:1062
virtual int v_GetNumPoints(const int dir) const
Definition: TriExp.cpp:1118
virtual void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
Definition: TriExp.cpp:413
virtual void v_GetEdgeInterpVals(const int edge, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:826
virtual StdRegions::Orientation v_GetTraceOrient(int edge)
Definition: TriExp.cpp:1107
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const
Definition: TriExp.cpp:701
virtual void v_GetTraceQFactors(const int edge, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:834
void v_DropLocStaticCondMatrix(const MatrixKey &mkey)
Definition: TriExp.cpp:1165
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
Integrates the specified function over the domain.
Definition: TriExp.cpp:79
virtual void v_MassLevelCurvatureMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1207
virtual void v_ComputeLaplacianMetric()
Definition: TriExp.cpp:1300
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1123
virtual void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1512
virtual void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp)
Definition: TriExp.cpp:1244
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1145
virtual void v_IProductWRTBase_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:440
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: TriExp.cpp:257
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: TriExp.cpp:400
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals)
This function evaluates the expansion at a single (arbitrary) point of the domain.
Definition: TriExp.cpp:756
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals)
Definition: TriExp.cpp:748
virtual void v_NormalTraceDerivFactors(Array< OneD, Array< OneD, NekDouble >> &factors, Array< OneD, Array< OneD, NekDouble >> &d0factors, Array< OneD, Array< OneD, NekDouble >> &d1factors)
: This method gets all of the factors which are required as part of the Gradient Jump Penalty stabili...
Definition: TriExp.cpp:1544
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3)
Definition: TriExp.cpp:736
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1214
virtual void v_GetTracePhysVals(const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient)
Definition: TriExp.cpp:767
virtual void v_IProductWRTDirectionalDerivBase(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:588
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const
Definition: TriExp.cpp:708
virtual void v_IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Directinoal Derivative in the modal space in the dir direction of varcoeffs.
Definition: TriExp.cpp:600
virtual void v_GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1221
virtual void v_WeakDirectionalDerivMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TriExp.cpp:1200
virtual void v_IProductWRTDerivBase_MatOp(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:549
virtual void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: TriExp.cpp:453
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
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)
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:140
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:163
void HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
const LibUtilities::PointsKeyVector GetPointsKeys() const
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:760
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:621
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
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:213
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:226
void IProductWRTDirectionalDerivBase_SumFac(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true)
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:731
void PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray)
Definition: StdExpansion.h:850
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
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 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
std::shared_ptr< Basis > BasisSharedPtr
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
Definition: Interp.cpp:52
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:250
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:53
@ 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
std::shared_ptr< SegExp > SegExpSharedPtr
Definition: SegExp.h:255
std::shared_ptr< GeomFactors > GeomFactorsSharedPtr
Pointer to a GeomFactors object.
Definition: GeomFactors.h:62
GeomType
Indicates the type of element geometry.
@ eRegular
Geometry is straight-sided with constant geometric factors.
@ eMovingRegular
Currently unused.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
Definition: Geometry.h:65
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< StdTriExp > StdTriExpSharedPtr
Definition: StdTriExp.h:235
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.cpp:534
void Svtsvtp(int n, const T alpha, const T *x, int incx, const T beta, const T *y, int incy, T *z, int incz)
vvtvvtp (scalar times vector plus scalar times vector):
Definition: Vmath.cpp:751
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 Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
Definition: Vmath.cpp:622
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 Vdiv(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:284
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 Reverse(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1286
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.cpp:692
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