Nektar++
SegExp.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File SegExp.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: SegExp routines
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <boost/core/ignore_unused.hpp>
36 
39 #include <LocalRegions/SegExp.h>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45 namespace LocalRegions
46 {
47 
48 /**
49  * @class SegExp
50  * Defines a Segment local expansion.
51  */
52 
53 /// Constructor using BasisKey class for quadrature points and
54 /// order definition.
55 /**
56  * @param Ba Basis key of segment expansion.
57  * @param geom Description of geometry.
58  */
59 SegExp::SegExp(const LibUtilities::BasisKey &Ba,
61  : StdExpansion(Ba.GetNumModes(), 1, Ba),
62  StdExpansion1D(Ba.GetNumModes(), Ba), StdRegions::StdSegExp(Ba),
63  Expansion(geom), Expansion1D(geom),
64  m_matrixManager(
65  std::bind(&SegExp::CreateMatrix, this, std::placeholders::_1),
66  std::string("SegExpMatrix")),
67  m_staticCondMatrixManager(std::bind(&Expansion::CreateStaticCondMatrix,
68  this, std::placeholders::_1),
69  std::string("SegExpStaticCondMatrix"))
70 {
71 }
72 
73 /// Copy Constructor
74 /**
75  * @param S Existing segment to duplicate.
76  */
78  : StdExpansion(S), StdExpansion1D(S), StdRegions::StdSegExp(S),
79  Expansion(S), Expansion1D(S), m_matrixManager(S.m_matrixManager),
80  m_staticCondMatrixManager(S.m_staticCondMatrixManager)
81 {
82 }
83 
84 /**
85  *
86  */
88 {
89 }
90 
91 //----------------------------
92 // Integration Methods
93 //----------------------------
94 
95 /** \brief Integrate the physical point list \a inarray over region
96  and return the value
97 
98  Inputs:\n
99 
100  - \a inarray: definition of function to be returned at
101  quadrature point of expansion.
102 
103  Outputs:\n
104 
105  - returns \f$\int^1_{-1} u(\xi_1)d \xi_1 \f$ where \f$inarray[i]
106  = u(\xi_{1i}) \f$
107 */
108 
110 {
111  int nquad0 = m_base[0]->GetNumPoints();
113  NekDouble ival;
114  Array<OneD, NekDouble> tmp(nquad0);
115 
116  // multiply inarray with Jacobian
117  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
118  {
119  Vmath::Vmul(nquad0, jac, 1, inarray, 1, tmp, 1);
120  }
121  else
122  {
123  Vmath::Smul(nquad0, jac[0], inarray, 1, tmp, 1);
124  }
125 
126  // call StdSegExp version;
127  ival = StdSegExp::v_Integral(tmp);
128  // ival = StdSegExp::Integral(tmp);
129  return ival;
130 }
131 
132 //-----------------------------
133 // Differentiation Methods
134 //-----------------------------
135 
136 /** \brief Evaluate the derivative \f$ d/d{\xi_1} \f$ at the
137  physical quadrature points given by \a inarray and return in \a
138  outarray.
139 
140  This is a wrapper around StdExpansion1D::Tensor_Deriv
141 
142  Input:\n
143 
144  - \a n: number of derivatives to be evaluated where \f$ n \leq dim\f$
145 
146  - \a inarray: array of function evaluated at the quadrature points
147 
148  Output: \n
149 
150  - \a outarray: array of the derivatives \f$
151  du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dx,
152  du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dy,
153  du/d_{\xi_1}|_{\xi_{1i}} d\xi_1/dz,
154  \f$ depending on value of \a dim
155 */
157  Array<OneD, NekDouble> &out_d0,
158  Array<OneD, NekDouble> &out_d1,
159  Array<OneD, NekDouble> &out_d2)
160 {
161  int nquad0 = m_base[0]->GetNumPoints();
163  m_metricinfo->GetDerivFactors(GetPointsKeys());
164  Array<OneD, NekDouble> diff(nquad0);
165 
166  // StdExpansion1D::PhysTensorDeriv(inarray,diff);
167  PhysTensorDeriv(inarray, diff);
168  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
169  {
170  if (out_d0.size())
171  {
172  Vmath::Vmul(nquad0, &gmat[0][0], 1, &diff[0], 1, &out_d0[0], 1);
173  }
174 
175  if (out_d1.size())
176  {
177  Vmath::Vmul(nquad0, &gmat[1][0], 1, &diff[0], 1, &out_d1[0], 1);
178  }
179 
180  if (out_d2.size())
181  {
182  Vmath::Vmul(nquad0, &gmat[2][0], 1, &diff[0], 1, &out_d2[0], 1);
183  }
184  }
185  else
186  {
187  if (out_d0.size())
188  {
189  Vmath::Smul(nquad0, gmat[0][0], diff, 1, out_d0, 1);
190  }
191 
192  if (out_d1.size())
193  {
194  Vmath::Smul(nquad0, gmat[1][0], diff, 1, out_d1, 1);
195  }
196 
197  if (out_d2.size())
198  {
199  Vmath::Smul(nquad0, gmat[2][0], diff, 1, out_d2, 1);
200  }
201  }
202 }
203 
204 /**
205  *\brief Evaluate the derivative along a line:
206  * \f$ d/ds=\frac{spacedim}{||tangent||}d/d{\xi} \f$.
207  * The derivative is calculated performing
208  *the product \f$ du/d{s}=\nabla u \cdot tangent \f$.
209  *\param inarray function to derive
210  *\param out_ds result of the derivative operation
211  **/
213  Array<OneD, NekDouble> &out_ds)
214 {
215  int nquad0 = m_base[0]->GetNumPoints();
216  int coordim = m_geom->GetCoordim();
217  Array<OneD, NekDouble> diff(nquad0);
218  // this operation is needed if you put out_ds==inarray
219  Vmath::Zero(nquad0, out_ds, 1);
220  switch (coordim)
221  {
222  case 2:
223  // diff= dU/de
224  Array<OneD, NekDouble> diff(nquad0);
225 
226  PhysTensorDeriv(inarray, diff);
227 
228  // get dS/de= (Jac)^-1
230  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
231  {
232  // calculate the derivative as (dU/de)*(Jac)^-1
233  Vmath::Vdiv(nquad0, diff, 1, Jac, 1, out_ds, 1);
234  }
235  else
236  {
237  NekDouble invJac = 1 / Jac[0];
238  Vmath::Smul(nquad0, invJac, diff, 1, out_ds, 1);
239  }
240  }
241 }
242 
243 /**
244  *\brief Evaluate the derivative normal to a line:
245  * \f$ d/dn=\frac{spacedim}{||normal||}d/d{\xi} \f$.
246  * The derivative is calculated performing
247  *the product \f$ du/d{s}=\nabla u \cdot normal \f$.
248  *\param inarray function to derive
249  *\param out_dn result of the derivative operation
250  **/
252  Array<OneD, NekDouble> &out_dn)
253 {
254  int nquad0 = m_base[0]->GetNumPoints();
256  m_metricinfo->GetDerivFactors(GetPointsKeys());
257  int coordim = m_geom->GetCoordim();
258  Array<OneD, NekDouble> out_dn_tmp(nquad0, 0.0);
259  switch (coordim)
260  {
261  case 2:
262 
263  Array<OneD, NekDouble> inarray_d0(nquad0);
264  Array<OneD, NekDouble> inarray_d1(nquad0);
265 
266  v_PhysDeriv(inarray, inarray_d0, inarray_d1);
268  normals = Array<OneD, Array<OneD, NekDouble>>(coordim);
269  cout << "der_n" << endl;
270  for (int k = 0; k < coordim; ++k)
271  {
272  normals[k] = Array<OneD, NekDouble>(nquad0);
273  }
274  // @TODO: this routine no longer makes sense, since normals are not
275  // unique on
276  // an edge
277  // normals = GetMetricInfo()->GetNormal();
278  for (int i = 0; i < nquad0; i++)
279  {
280  cout << "nx= " << normals[0][i] << " ny=" << normals[1][i]
281  << endl;
282  }
284  "normal vectors do not exist: check if a"
285  "boundary region is defined as I ");
286  // \nabla u \cdot normal
287  Vmath::Vmul(nquad0, normals[0], 1, inarray_d0, 1, out_dn_tmp, 1);
288  Vmath::Vadd(nquad0, out_dn_tmp, 1, out_dn, 1, out_dn, 1);
289  Vmath::Zero(nquad0, out_dn_tmp, 1);
290  Vmath::Vmul(nquad0, normals[1], 1, inarray_d1, 1, out_dn_tmp, 1);
291  Vmath::Vadd(nquad0, out_dn_tmp, 1, out_dn, 1, out_dn, 1);
292 
293  for (int i = 0; i < nquad0; i++)
294  {
295  cout << "deps/dx =" << inarray_d0[i]
296  << " deps/dy=" << inarray_d1[i] << endl;
297  }
298  }
299 }
300 void SegExp::v_PhysDeriv(const int dir,
301  const Array<OneD, const NekDouble> &inarray,
302  Array<OneD, NekDouble> &outarray)
303 {
304  switch (dir)
305  {
306  case 0:
307  {
308  PhysDeriv(inarray, outarray, NullNekDouble1DArray,
310  }
311  break;
312  case 1:
313  {
314  PhysDeriv(inarray, NullNekDouble1DArray, outarray,
316  }
317  break;
318  case 2:
319  {
321  outarray);
322  }
323  break;
324  default:
325  {
326  ASSERTL1(false, "input dir is out of range");
327  }
328  break;
329  }
330 }
331 
332 //-----------------------------
333 // Transforms
334 //-----------------------------
335 
336 /** \brief Forward transform from physical quadrature space
337  stored in \a inarray and evaluate the expansion coefficients and
338  store in \a outarray
339 
340  Perform a forward transform using a Galerkin projection by
341  taking the inner product of the physical points and multiplying
342  by the inverse of the mass matrix using the Solve method of the
343  standard matrix container holding the local mass matrix, i.e.
344  \f$ {\bf \hat{u}} = {\bf M}^{-1} {\bf I} \f$ where \f$ {\bf I}[p] =
345  \int^1_{-1} \phi_p(\xi_1) u(\xi_1) d\xi_1 \f$
346 
347  Inputs:\n
348 
349  - \a inarray: array of physical quadrature points to be transformed
350 
351  Outputs:\n
352 
353  - \a outarray: updated array of expansion coefficients.
354 
355 */
356 // need to sort out family of matrices
358  Array<OneD, NekDouble> &outarray)
359 {
360  if (m_base[0]->Collocation())
361  {
362  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
363  }
364  else
365  {
366  v_IProductWRTBase(inarray, outarray);
367 
368  // get Mass matrix inverse
369  MatrixKey masskey(StdRegions::eInvMass, DetShapeType(), *this);
370  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
371 
372  // copy inarray in case inarray == outarray
373  NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
374  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
375 
376  out = (*matsys) * in;
377  }
378 }
379 
381  const Array<OneD, const NekDouble> &inarray,
382  Array<OneD, NekDouble> &outarray)
383 {
384  if (m_base[0]->Collocation())
385  {
386  Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
387  }
388  else
389  {
390  int nInteriorDofs = m_ncoeffs - 2;
391  int offset = 0;
392 
393  switch (m_base[0]->GetBasisType())
394  {
396  {
397  offset = 1;
398  }
399  break;
401  {
402  nInteriorDofs = m_ncoeffs;
403  offset = 0;
404  }
405  break;
408  {
409  ASSERTL1(
410  m_base[0]->GetPointsType() ==
412  m_base[0]->GetPointsType() ==
414  "Cannot use FwdTrans_BndConstrained with these points.");
415  offset = 2;
416  }
417  break;
418  default:
419  ASSERTL0(false, "This type of FwdTrans is not defined"
420  "for this expansion type");
421  }
422 
423  fill(outarray.get(), outarray.get() + m_ncoeffs, 0.0);
424 
426  {
427 
428  outarray[GetVertexMap(0)] = inarray[0];
429  outarray[GetVertexMap(1)] = inarray[m_base[0]->GetNumPoints() - 1];
430 
431  if (m_ncoeffs > 2)
432  {
433  // ideally, we would like to have tmp0 to be replaced
434  // by outarray (currently MassMatrixOp does not allow
435  // aliasing)
438 
440  DetShapeType(), *this);
441  MassMatrixOp(outarray, tmp0, stdmasskey);
442  v_IProductWRTBase(inarray, tmp1);
443 
444  Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
445 
446  // get Mass matrix inverse (only of interior DOF)
447  MatrixKey masskey(StdRegions::eMass, DetShapeType(), *this);
448  DNekScalMatSharedPtr matsys =
449  (m_staticCondMatrixManager[masskey])->GetBlock(1, 1);
450 
451  Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, matsys->Scale(),
452  &((matsys->GetOwnedMatrix())->GetPtr())[0],
453  nInteriorDofs, tmp1.get() + offset, 1, 0.0,
454  outarray.get() + offset, 1);
455  }
456  }
457  else
458  {
459  SegExp::v_FwdTrans(inarray, outarray);
460  }
461  }
462 }
463 
464 //-----------------------------
465 // Inner product functions
466 //-----------------------------
467 
468 /** \brief Inner product of \a inarray over region with respect to
469  the expansion basis (this)->_Base[0] and return in \a outarray
470 
471  Wrapper call to SegExp::IProduct_WRT_B
472 
473  Input:\n
474 
475  - \a inarray: array of function evaluated at the physical
476  collocation points
477 
478  Output:\n
479 
480  - \a outarray: array of inner product with respect to each
481  basis over region
482 */
484  Array<OneD, NekDouble> &outarray)
485 {
486  v_IProductWRTBase(m_base[0]->GetBdata(), inarray, outarray, 1);
487 }
488 
489 /**
490  \brief Inner product of \a inarray over region with respect to
491  expansion basis \a base and return in \a outarray
492 
493  Calculate \f$ I[p] = \int^{1}_{-1} \phi_p(\xi_1) u(\xi_1) d\xi_1
494  = \sum_{i=0}^{nq-1} \phi_p(\xi_{1i}) u(\xi_{1i}) w_i \f$ where
495  \f$ outarray[p] = I[p], inarray[i] = u(\xi_{1i}), base[p*nq+i] =
496  \phi_p(\xi_{1i}) \f$.
497 
498  Inputs: \n
499 
500  - \a base: an array definiing the local basis for the inner
501  product usually passed from Basis->get_bdata() or
502  Basis->get_Dbdata()
503  - \a inarray: physical point array of function to be integrated
504  \f$ u(\xi_1) \f$
505  - \a coll_check: Flag to identify when a Basis->collocation()
506  call should be performed to see if this is a GLL_Lagrange basis
507  with a collocation property. (should be set to 0 if taking the
508  inner product with respect to the derivative of basis)
509 
510  Output: \n
511 
512  - \a outarray: array of coefficients representing the inner
513  product of function with ever mode in the exapnsion
514 
515 **/
517  const Array<OneD, const NekDouble> &inarray,
518  Array<OneD, NekDouble> &outarray, int coll_check)
519 {
520  int nquad0 = m_base[0]->GetNumPoints();
522  Array<OneD, NekDouble> tmp(nquad0);
523 
524  // multiply inarray with Jacobian
525  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
526  {
527  Vmath::Vmul(nquad0, jac, 1, inarray, 1, tmp, 1);
528  }
529  else
530  {
531  Vmath::Smul(nquad0, jac[0], inarray, 1, tmp, 1);
532  }
533  StdSegExp::v_IProductWRTBase(base, tmp, outarray, coll_check);
534 }
535 
537  const Array<OneD, const NekDouble> &inarray,
538  Array<OneD, NekDouble> &outarray)
539 {
540  ASSERTL1(dir < 3, "input dir is out of range");
541  ASSERTL1((dir == 2) ? m_geom->GetCoordim() == 3 : true,
542  "input dir is out of range");
543 
544  int nquad = m_base[0]->GetNumPoints();
545  const Array<TwoD, const NekDouble> &gmat =
546  m_metricinfo->GetDerivFactors(GetPointsKeys());
547 
548  Array<OneD, NekDouble> tmp1(nquad);
549 
550  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
551  {
552  Vmath::Vmul(nquad, gmat[dir], 1, inarray, 1, tmp1, 1);
553  }
554  else
555  {
556  Vmath::Smul(nquad, gmat[dir][0], inarray, 1, tmp1, 1);
557  }
558 
559  v_IProductWRTBase(m_base[0]->GetDbdata(), tmp1, outarray, 1);
560 }
561 
564  Array<OneD, NekDouble> &outarray)
565 {
566  int nq = m_base[0]->GetNumPoints();
567  Array<OneD, NekDouble> Fn(nq);
568 
569  // @TODO: This routine no longer makes sense as a normal is not unique to an
570  // edge
572  GetLeftAdjacentElementExp()->GetTraceNormal(
574  Vmath::Vmul(nq, &Fx[0], 1, &normals[0][0], 1, &Fn[0], 1);
575  Vmath::Vvtvp(nq, &Fy[0], 1, &normals[1][0], 1, &Fn[0], 1, &Fn[0], 1);
576 
577  v_IProductWRTBase(Fn, outarray);
578 }
579 
581  const Array<OneD, const Array<OneD, NekDouble>> &Fvec,
582  Array<OneD, NekDouble> &outarray)
583 {
584  NormVectorIProductWRTBase(Fvec[0], Fvec[1], outarray);
585 }
586 
587 //-----------------------------
588 // Evaluation functions
589 //-----------------------------
590 
591 /**
592  * Given the local cartesian coordinate \a Lcoord evaluate the
593  * value of physvals at this point by calling through to the
594  * StdExpansion method
595  */
597  const Array<OneD, const NekDouble> &Lcoord,
598  const Array<OneD, const NekDouble> &physvals)
599 {
600  // Evaluate point in local (eta) coordinates.
601  return StdSegExp::v_PhysEvaluate(Lcoord, physvals);
602 }
603 
605  const Array<OneD, const NekDouble> &physvals)
606 {
608 
609  ASSERTL0(m_geom, "m_geom not defined");
610  m_geom->GetLocCoords(coord, Lcoord);
611 
612  return StdSegExp::v_PhysEvaluate(Lcoord, physvals);
613 }
614 
616  Array<OneD, NekDouble> &coords)
617 {
618  int i;
619 
620  ASSERTL1(Lcoords[0] >= -1.0 && Lcoords[0] <= 1.0,
621  "Local coordinates are not in region [-1,1]");
622 
623  m_geom->FillGeom();
624  for (i = 0; i < m_geom->GetCoordim(); ++i)
625  {
626  coords[i] = m_geom->GetCoord(i, Lcoords);
627  }
628 }
629 
631  Array<OneD, NekDouble> &coords_1,
632  Array<OneD, NekDouble> &coords_2)
633 {
634  Expansion::v_GetCoords(coords_0, coords_1, coords_2);
635 }
636 
637 // Get vertex value from the 1D Phys space.
638 void SegExp::v_GetVertexPhysVals(const int vertex,
639  const Array<OneD, const NekDouble> &inarray,
640  NekDouble &outarray)
641 {
642  int nquad = m_base[0]->GetNumPoints();
643 
645  {
646  switch (vertex)
647  {
648  case 0:
649  outarray = inarray[0];
650  break;
651  case 1:
652  outarray = inarray[nquad - 1];
653  break;
654  }
655  }
656  else
657  {
659  factors[StdRegions::eFactorGaussVertex] = vertex;
660 
662  *this, factors);
663 
664  DNekScalMatSharedPtr mat_gauss = m_matrixManager[key];
665 
666  outarray =
667  Blas::Ddot(nquad, mat_gauss->GetOwnedMatrix()->GetPtr().get(), 1,
668  &inarray[0], 1);
669  }
670 }
671 
672 // Add vertex value of the 1D Phys space to outarray.
673 void SegExp::v_AddVertexPhysVals(const int vertex, const NekDouble &inarray,
674  Array<OneD, NekDouble> &outarray)
675 {
676  size_t nquad = m_base[0]->GetNumPoints();
677 
679  {
680  switch (vertex)
681  {
682  case 0:
683  outarray[0] += inarray;
684  break;
685  case 1:
686  outarray[nquad - 1] += inarray;
687  break;
688  }
689  }
690  else
691  {
693  factors[StdRegions::eFactorGaussVertex] = vertex;
694 
696  *this, factors);
697 
698  DNekScalMatSharedPtr mat_gauss = m_matrixManager[key];
699 
700  Vmath::Svtvp(nquad, inarray,
701  mat_gauss->GetOwnedMatrix()->GetPtr().get(), 1,
702  &outarray[0], 1, &outarray[0], 1);
703  }
704 }
705 // Get vertex value from the 1D Phys space.
707  const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp,
708  const Array<OneD, const NekDouble> &inarray,
710 {
711  boost::ignore_unused(EdgeExp, orient);
712 
713  NekDouble result;
714  v_GetVertexPhysVals(edge, inarray, result);
715  outarray[0] = result;
716 }
717 
718 // Get vertex map from the 1D Phys space.
719 void SegExp::v_GetTracePhysMap(const int vertex, Array<OneD, int> &map)
720 {
721  int nquad = m_base[0]->GetNumPoints();
722 
723  ASSERTL1(vertex == 0 || vertex == 1, "Vertex value should be 0 or 1");
724 
725  map = Array<OneD, int>(1);
726 
727  map[0] = vertex == 0 ? 0 : nquad - 1;
728 }
729 
730 //-----------------------------
731 // Helper functions
732 //-----------------------------
733 
736  Array<OneD, NekDouble> &outarray)
737 {
738 
739  if (dir == StdRegions::eBackwards)
740  {
741  if (&inarray[0] != &outarray[0])
742  {
743  Array<OneD, NekDouble> intmp(inarray);
744  ReverseCoeffsAndSign(intmp, outarray);
745  }
746  else
747  {
748  ReverseCoeffsAndSign(inarray, outarray);
749  }
750  }
751 }
752 
754 {
756  m_base[0]->GetBasisKey());
757 }
758 
760 {
762  m_base[0]->GetPointsKey());
763 
765 }
766 
768 {
769  return m_geom->GetCoordim();
770 }
771 
773 {
774  NEKERROR(ErrorUtil::efatal, "Got to SegExp");
775  return NullNekDouble1DArray;
776 }
777 
779 {
780  return m_metricinfo->GetGtype();
781 }
782 
783 int SegExp::v_GetNumPoints(const int dir) const
784 {
785  return GetNumPoints(dir);
786 }
787 
788 int SegExp::v_GetNcoeffs(void) const
789 {
790  return m_ncoeffs;
791 }
792 
794 {
795  return GetBasis(dir);
796 }
797 
799 {
800  return 2;
801 }
802 
804 {
805  return 2;
806 }
807 
808 /// Unpack data from input file assuming it comes from
809 // the same expansion type
811  const NekDouble *data, const std::vector<unsigned int> &nummodes,
812  const int mode_offset, NekDouble *coeffs,
813  std::vector<LibUtilities::BasisType> &fromType)
814 {
815  boost::ignore_unused(fromType);
816 
817  switch (m_base[0]->GetBasisType())
818  {
820  {
821  int fillorder = min((int)nummodes[mode_offset], m_ncoeffs);
822 
823  Vmath::Zero(m_ncoeffs, coeffs, 1);
824  Vmath::Vcopy(fillorder, &data[0], 1, &coeffs[0], 1);
825  }
826  break;
828  {
829  // Assume that input is also Gll_Lagrange
830  // but no way to check;
831  LibUtilities::PointsKey f0(nummodes[mode_offset],
833  LibUtilities::PointsKey t0(m_base[0]->GetNumModes(),
835  LibUtilities::Interp1D(f0, data, t0, coeffs);
836  }
837  break;
839  {
840  // Assume that input is also Gauss_Lagrange
841  // but no way to check;
842  LibUtilities::PointsKey f0(nummodes[mode_offset],
844  LibUtilities::PointsKey t0(m_base[0]->GetNumModes(),
846  LibUtilities::Interp1D(f0, data, t0, coeffs);
847  }
848  break;
849  default:
850  ASSERTL0(false, "basis is either not set up or not hierarchicial");
851  }
852 }
853 
854 void SegExp::v_ComputeTraceNormal(const int vertex)
855 {
856  int i;
857  const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
858  GetGeom()->GetMetricInfo();
859  SpatialDomains::GeomType type = geomFactors->GetGtype();
860  const Array<TwoD, const NekDouble> &gmat =
861  geomFactors->GetDerivFactors(GetPointsKeys());
862  int nqe = 1;
863  int vCoordDim = GetCoordim();
864 
867  for (i = 0; i < vCoordDim; ++i)
868  {
869  normal[i] = Array<OneD, NekDouble>(nqe);
870  }
871 
872  size_t nqb = nqe;
873  size_t nbnd = vertex;
876 
877  // Regular geometry case
878  if ((type == SpatialDomains::eRegular) ||
880  {
881  NekDouble vert;
882  // Set up normals
883  switch (vertex)
884  {
885  case 0:
886  for (i = 0; i < vCoordDim; ++i)
887  {
888  Vmath::Fill(nqe, -gmat[i][0], normal[i], 1);
889  }
890  break;
891  case 1:
892  for (i = 0; i < vCoordDim; ++i)
893  {
894  Vmath::Fill(nqe, gmat[i][0], normal[i], 1);
895  }
896  break;
897  default:
898  ASSERTL0(false, "point is out of range (point < 2)");
899  }
900 
901  // normalise
902  vert = 0.0;
903  for (i = 0; i < vCoordDim; ++i)
904  {
905  vert += normal[i][0] * normal[i][0];
906  }
907  vert = 1.0 / sqrt(vert);
908 
909  Vmath::Fill(nqb, vert, length, 1);
910 
911  for (i = 0; i < vCoordDim; ++i)
912  {
913  Vmath::Smul(nqe, vert, normal[i], 1, normal[i], 1);
914  }
915  }
916 }
917 
918 //-----------------------------
919 // Operator creation functions
920 //-----------------------------
921 
923  Array<OneD, NekDouble> &outarray,
924  const StdRegions::StdMatrixKey &mkey)
925 {
926  boost::ignore_unused(mkey);
927 
928  int nquad = m_base[0]->GetNumPoints();
929  const Array<TwoD, const NekDouble> &gmat =
930  m_metricinfo->GetDerivFactors(GetPointsKeys());
931 
932  Array<OneD, NekDouble> physValues(nquad);
933  Array<OneD, NekDouble> dPhysValuesdx(nquad);
934 
935  BwdTrans(inarray, physValues);
936 
937  // Laplacian matrix operation
938  switch (m_geom->GetCoordim())
939  {
940  case 1:
941  {
942  PhysDeriv(physValues, dPhysValuesdx);
943 
944  // multiply with the proper geometric factors
945  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
946  {
947  Vmath::Vmul(nquad, &gmat[0][0], 1, dPhysValuesdx.get(), 1,
948  dPhysValuesdx.get(), 1);
949  }
950  else
951  {
952  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
953  }
954  }
955  break;
956  case 2:
957  {
958  Array<OneD, NekDouble> dPhysValuesdy(nquad);
959 
960  PhysDeriv(physValues, dPhysValuesdx, dPhysValuesdy);
961 
962  // multiply with the proper geometric factors
963  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
964  {
965  Vmath::Vmul(nquad, &gmat[0][0], 1, dPhysValuesdx.get(), 1,
966  dPhysValuesdx.get(), 1);
967  Vmath::Vvtvp(nquad, &gmat[1][0], 1, dPhysValuesdy.get(), 1,
968  dPhysValuesdx.get(), 1, dPhysValuesdx.get(), 1);
969  }
970  else
971  {
972  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
973  Blas::Daxpy(nquad, gmat[1][0], dPhysValuesdy.get(), 1,
974  dPhysValuesdx.get(), 1);
975  }
976  }
977  break;
978  case 3:
979  {
980  Array<OneD, NekDouble> dPhysValuesdy(nquad);
981  Array<OneD, NekDouble> dPhysValuesdz(nquad);
982 
983  PhysDeriv(physValues, dPhysValuesdx, dPhysValuesdy, dPhysValuesdz);
984 
985  // multiply with the proper geometric factors
986  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
987  {
988  Vmath::Vmul(nquad, &gmat[0][0], 1, dPhysValuesdx.get(), 1,
989  dPhysValuesdx.get(), 1);
990  Vmath::Vvtvp(nquad, &gmat[1][0], 1, dPhysValuesdy.get(), 1,
991  dPhysValuesdx.get(), 1, dPhysValuesdx.get(), 1);
992  Vmath::Vvtvp(nquad, &gmat[2][0], 1, dPhysValuesdz.get(), 1,
993  dPhysValuesdx.get(), 1, dPhysValuesdx.get(), 1);
994  }
995  else
996  {
997  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
998  Blas::Daxpy(nquad, gmat[1][0], dPhysValuesdy.get(), 1,
999  dPhysValuesdx.get(), 1);
1000  Blas::Daxpy(nquad, gmat[2][0], dPhysValuesdz.get(), 1,
1001  dPhysValuesdx.get(), 1);
1002  }
1003  }
1004  break;
1005  default:
1006  ASSERTL0(false, "Wrong number of dimensions");
1007  break;
1008  }
1009 
1010  v_IProductWRTBase(m_base[0]->GetDbdata(), dPhysValuesdx, outarray, 1);
1011 }
1012 
1014  Array<OneD, NekDouble> &outarray,
1015  const StdRegions::StdMatrixKey &mkey)
1016 {
1017  int nquad = m_base[0]->GetNumPoints();
1018  const Array<TwoD, const NekDouble> &gmat =
1019  m_metricinfo->GetDerivFactors(GetPointsKeys());
1020  const NekDouble lambda = mkey.GetConstFactor(StdRegions::eFactorLambda);
1021 
1022  Array<OneD, NekDouble> physValues(nquad);
1023  Array<OneD, NekDouble> dPhysValuesdx(nquad);
1025 
1026  BwdTrans(inarray, physValues);
1027 
1028  // mass matrix operation
1029  v_IProductWRTBase((m_base[0]->GetBdata()), physValues, wsp, 1);
1030 
1031  // Laplacian matrix operation
1032  switch (m_geom->GetCoordim())
1033  {
1034  case 1:
1035  {
1036  PhysDeriv(physValues, dPhysValuesdx);
1037 
1038  // multiply with the proper geometric factors
1039  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1040  {
1041  Vmath::Vmul(nquad, &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1042  dPhysValuesdx.get(), 1);
1043  }
1044  else
1045  {
1046  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1047  }
1048  }
1049  break;
1050  case 2:
1051  {
1052  Array<OneD, NekDouble> dPhysValuesdy(nquad);
1053 
1054  PhysDeriv(physValues, dPhysValuesdx, dPhysValuesdy);
1055 
1056  // multiply with the proper geometric factors
1057  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1058  {
1059  Vmath::Vmul(nquad, &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1060  dPhysValuesdx.get(), 1);
1061  Vmath::Vvtvp(nquad, &gmat[1][0], 1, dPhysValuesdy.get(), 1,
1062  dPhysValuesdx.get(), 1, dPhysValuesdx.get(), 1);
1063  }
1064  else
1065  {
1066  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1067  Blas::Daxpy(nquad, gmat[1][0], dPhysValuesdy.get(), 1,
1068  dPhysValuesdx.get(), 1);
1069  }
1070  }
1071  break;
1072  case 3:
1073  {
1074  Array<OneD, NekDouble> dPhysValuesdy(nquad);
1075  Array<OneD, NekDouble> dPhysValuesdz(nquad);
1076 
1077  PhysDeriv(physValues, dPhysValuesdx, dPhysValuesdy, dPhysValuesdz);
1078 
1079  // multiply with the proper geometric factors
1080  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1081  {
1082  Vmath::Vmul(nquad, &gmat[0][0], 1, dPhysValuesdx.get(), 1,
1083  dPhysValuesdx.get(), 1);
1084  Vmath::Vvtvp(nquad, &gmat[1][0], 1, dPhysValuesdy.get(), 1,
1085  dPhysValuesdx.get(), 1, dPhysValuesdx.get(), 1);
1086  Vmath::Vvtvp(nquad, &gmat[2][0], 1, dPhysValuesdz.get(), 1,
1087  dPhysValuesdx.get(), 1, dPhysValuesdx.get(), 1);
1088  }
1089  else
1090  {
1091  Blas::Dscal(nquad, gmat[0][0], dPhysValuesdx.get(), 1);
1092  Blas::Daxpy(nquad, gmat[1][0], dPhysValuesdy.get(), 1,
1093  dPhysValuesdx.get(), 1);
1094  Blas::Daxpy(nquad, gmat[2][0], dPhysValuesdz.get(), 1,
1095  dPhysValuesdx.get(), 1);
1096  }
1097  }
1098  break;
1099  default:
1100  ASSERTL0(false, "Wrong number of dimensions");
1101  break;
1102  }
1103 
1104  v_IProductWRTBase(m_base[0]->GetDbdata(), dPhysValuesdx, outarray, 1);
1105  Blas::Daxpy(m_ncoeffs, lambda, wsp.get(), 1, outarray.get(), 1);
1106 }
1107 
1108 //-----------------------------
1109 // Matrix creation functions
1110 //-----------------------------
1111 
1113 {
1114  return m_staticCondMatrixManager[mkey];
1115 }
1116 
1118 {
1119  m_staticCondMatrixManager.DeleteObject(mkey);
1120 }
1121 
1123 {
1124  return m_matrixManager[mkey];
1125 }
1126 
1128 {
1129  LibUtilities::BasisKey bkey = m_base[0]->GetBasisKey();
1132 
1133  return tmp->GetStdMatrix(mkey);
1134 }
1135 
1137 {
1138  DNekScalMatSharedPtr returnval;
1139  NekDouble fac;
1141 
1143  "Geometric information is not set up");
1144 
1145  switch (mkey.GetMatrixType())
1146  {
1147  case StdRegions::eMass:
1148  {
1149  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
1150  (mkey.GetNVarCoeff()))
1151  {
1152  fac = 1.0;
1153  goto UseLocRegionsMatrix;
1154  }
1155  else
1156  {
1157  fac = (m_metricinfo->GetJac(ptsKeys))[0];
1158  goto UseStdRegionsMatrix;
1159  }
1160  }
1161  break;
1162  case StdRegions::eInvMass:
1163  {
1164  if ((m_metricinfo->GetGtype() == SpatialDomains::eDeformed) ||
1165  (mkey.GetNVarCoeff()))
1166  {
1167  NekDouble one = 1.0;
1169  DetShapeType(), *this);
1170  DNekMatSharedPtr mat = GenMatrix(masskey);
1171  mat->Invert();
1172 
1173  returnval =
1175  }
1176  else
1177  {
1178  fac = 1.0 / (m_metricinfo->GetJac(ptsKeys))[0];
1179  goto UseStdRegionsMatrix;
1180  }
1181  }
1182  break;
1186  {
1187  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed ||
1188  mkey.GetNVarCoeff())
1189  {
1190  fac = 1.0;
1191  goto UseLocRegionsMatrix;
1192  }
1193  else
1194  {
1195  int dir = 0;
1196  switch (mkey.GetMatrixType())
1197  {
1199  dir = 0;
1200  break;
1202  ASSERTL1(m_geom->GetCoordim() >= 2,
1203  "Cannot call eWeakDeriv2 in a "
1204  "coordinate system which is not at "
1205  "least two-dimensional");
1206  dir = 1;
1207  break;
1209  ASSERTL1(m_geom->GetCoordim() == 3,
1210  "Cannot call eWeakDeriv2 in a "
1211  "coordinate system which is not "
1212  "three-dimensional");
1213  dir = 2;
1214  break;
1215  default:
1216  break;
1217  }
1218 
1220  mkey.GetShapeType(), *this);
1221 
1222  DNekMatSharedPtr WeakDerivStd = GetStdMatrix(deriv0key);
1223  fac = m_metricinfo->GetDerivFactors(ptsKeys)[dir][0] *
1224  m_metricinfo->GetJac(ptsKeys)[0];
1225 
1227  fac, WeakDerivStd);
1228  }
1229  }
1230  break;
1232  {
1233  if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1234  {
1235  fac = 1.0;
1236  goto UseLocRegionsMatrix;
1237  }
1238  else
1239  {
1240  int coordim = m_geom->GetCoordim();
1241  fac = 0.0;
1242  for (int i = 0; i < coordim; ++i)
1243  {
1244  fac += m_metricinfo->GetDerivFactors(ptsKeys)[i][0] *
1245  m_metricinfo->GetDerivFactors(ptsKeys)[i][0];
1246  }
1247  fac *= m_metricinfo->GetJac(ptsKeys)[0];
1248  goto UseStdRegionsMatrix;
1249  }
1250  }
1251  break;
1253  {
1255  MatrixKey masskey(StdRegions::eMass, mkey.GetShapeType(), *this);
1256  DNekScalMat &MassMat = *(this->m_matrixManager[masskey]);
1257  MatrixKey lapkey(StdRegions::eLaplacian, mkey.GetShapeType(), *this,
1258  mkey.GetConstFactors(), mkey.GetVarCoeffs());
1259  DNekScalMat &LapMat = *(this->m_matrixManager[lapkey]);
1260 
1261  int rows = LapMat.GetRows();
1262  int cols = LapMat.GetColumns();
1263 
1264  DNekMatSharedPtr helm =
1266 
1267  NekDouble one = 1.0;
1268  (*helm) = LapMat + factor * MassMat;
1269 
1270  returnval =
1272  }
1273  break;
1278  {
1279  NekDouble one = 1.0;
1280 
1281  DNekMatSharedPtr mat = GenMatrix(mkey);
1282  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
1283  }
1284  break;
1286  {
1287  NekDouble one = 1.0;
1288 
1289  // StdRegions::StdMatrixKey
1290  // hkey(StdRegions::eHybridDGHelmholtz,
1291  // DetShapeType(),*this,
1292  // mkey.GetConstant(0),
1293  // mkey.GetConstant(1));
1295  *this, mkey.GetConstFactors(), mkey.GetVarCoeffs());
1296  DNekMatSharedPtr mat = GenMatrix(hkey);
1297 
1298  mat->Invert();
1299  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(one, mat);
1300  }
1301  break;
1303  {
1304  DNekMatSharedPtr m_Ix;
1305  Array<OneD, NekDouble> coords(1, 0.0);
1307  int vertex = (int)factors[StdRegions::eFactorGaussVertex];
1308 
1309  coords[0] = (vertex == 0) ? -1.0 : 1.0;
1310 
1311  m_Ix = m_base[0]->GetI(coords);
1312  returnval =
1314  }
1315  break;
1316 
1317  UseLocRegionsMatrix:
1318  {
1319  DNekMatSharedPtr mat = GenMatrix(mkey);
1320  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac, mat);
1321  }
1322  break;
1323  UseStdRegionsMatrix:
1324  {
1325  DNekMatSharedPtr mat = GetStdMatrix(mkey);
1326  returnval = MemoryManager<DNekScalMat>::AllocateSharedPtr(fac, mat);
1327  }
1328  break;
1329  default:
1330  NEKERROR(ErrorUtil::efatal, "Matrix creation not defined");
1331  break;
1332  }
1333 
1334  return returnval;
1335 }
1336 
1338 {
1339  DNekMatSharedPtr returnval;
1340 
1341  switch (mkey.GetMatrixType())
1342  {
1349  returnval = Expansion1D::v_GenMatrix(mkey);
1350  break;
1351  default:
1352  returnval = StdSegExp::v_GenMatrix(mkey);
1353  break;
1354  }
1355 
1356  return returnval;
1357 }
1358 
1359 //-----------------------------
1360 // Private methods
1361 //-----------------------------
1362 
1363 /// Reverse the coefficients in a boundary interior expansion
1364 /// this routine is of use when we need the segment
1365 /// coefficients corresponding to a expansion in the reverse
1366 /// coordinate direction
1368  Array<OneD, NekDouble> &outarray)
1369 {
1370 
1371  int m;
1372  NekDouble sgn = 1;
1373 
1374  ASSERTL1(&inarray[0] != &outarray[0],
1375  "inarray and outarray can not be the same");
1376  switch (GetBasisType(0))
1377  {
1379  // Swap vertices
1380  outarray[0] = inarray[1];
1381  outarray[1] = inarray[0];
1382  // negate odd modes
1383  for (m = 2; m < m_ncoeffs; ++m)
1384  {
1385  outarray[m] = sgn * inarray[m];
1386  sgn = -sgn;
1387  }
1388  break;
1391  for (m = 0; m < m_ncoeffs; ++m)
1392  {
1393  outarray[m_ncoeffs - 1 - m] = inarray[m];
1394  }
1395  break;
1396  default:
1397  ASSERTL0(false, "This basis is not allowed in this method");
1398  break;
1399  }
1400 }
1401 
1402 /* \brief Mass inversion product from \a inarray to \a outarray
1403  *
1404  * Multiply by the inverse of the mass matrix
1405  * \f$ {\bf \hat{u}} = {\bf M}^{-1} {\bf I} \f$
1406  *
1407  **/
1409  Array<OneD, NekDouble> &outarray)
1410 {
1411  // get Mass matrix inverse
1412  MatrixKey masskey(StdRegions::eInvMass, DetShapeType(), *this);
1413  DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
1414 
1415  NekVector<NekDouble> in(m_ncoeffs, inarray, eCopy);
1416  NekVector<NekDouble> out(m_ncoeffs, outarray, eWrapper);
1417 
1418  out = (*matsys) * in;
1419 }
1420 
1421 } // namespace LocalRegions
1422 } // 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
Describes the specification for a Basis.
Definition: Basis.h:50
Defines a specification for a set of points.
Definition: Points.h:59
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
Definition: Expansion1D.cpp:54
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
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
virtual void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: SegExp.cpp:922
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Forward transform from physical quadrature space stored in inarray and evaluate the expansion coeffic...
Definition: SegExp.cpp:357
virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
Definition: SegExp.cpp:759
virtual void v_GetTracePhysVals(const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, StdRegions::Orientation orient) override
Definition: SegExp.cpp:706
virtual void v_PhysDeriv_s(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_ds) override
Evaluate the derivative along a line: . The derivative is calculated performing the product .
Definition: SegExp.cpp:212
virtual int v_GetNumPoints(const int dir) const
Definition: SegExp.cpp:783
virtual void v_GetVertexPhysVals(const int vertex, const Array< OneD, const NekDouble > &inarray, NekDouble &outarray) override
Definition: SegExp.cpp:638
virtual NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) override
Definition: SegExp.cpp:596
void ReverseCoeffsAndSign(const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Reverse the coefficients in a boundary interior expansion this routine is of use when we need the seg...
Definition: SegExp.cpp:1367
virtual void v_NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray) override
Definition: SegExp.cpp:562
DNekScalMatSharedPtr CreateMatrix(const MatrixKey &mkey)
Definition: SegExp.cpp:1136
virtual int v_GetCoordim() override
Definition: SegExp.cpp:767
virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
Definition: SegExp.cpp:753
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: SegExp.cpp:380
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Inner product of inarray over region with respect to the expansion basis (this)->_Base[0] and return ...
Definition: SegExp.cpp:483
virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: SegExp.cpp:1337
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: SegExp.h:240
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrate the physical point list inarray over region and return the value.
Definition: SegExp.cpp:109
virtual void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: SegExp.cpp:1013
virtual void v_PhysDeriv_n(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_dn) override
Evaluate the derivative normal to a line: . The derivative is calculated performing the product .
Definition: SegExp.cpp:251
virtual void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType) override
Unpack data from input file assuming it comes from.
Definition: SegExp.cpp:810
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray) override
Evaluate the derivative at the physical quadrature points given by inarray and return in outarray.
Definition: SegExp.cpp:156
virtual void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
Definition: SegExp.cpp:630
virtual DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: SegExp.cpp:1112
void MultiplyByElmtInvMass(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: SegExp.cpp:1408
virtual int v_NumBndryCoeffs() const override
Definition: SegExp.cpp:798
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: SegExp.h:238
virtual NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coord, const Array< OneD, const NekDouble > &physvals) override
Definition: SegExp.cpp:604
virtual void v_GetTracePhysMap(const int vertex, Array< OneD, int > &map) override
Definition: SegExp.cpp:719
virtual int v_NumDGBndryCoeffs() const override
Definition: SegExp.cpp:803
virtual DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: SegExp.cpp:1127
virtual void v_SetCoeffsToOrientation(StdRegions::Orientation dir, Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: SegExp.cpp:734
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: SegExp.cpp:1117
virtual DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
Definition: SegExp.cpp:1122
virtual void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
Definition: SegExp.cpp:615
virtual void v_AddVertexPhysVals(const int vertex, const NekDouble &inarray, Array< OneD, NekDouble > &outarray)
Definition: SegExp.cpp:673
virtual int v_GetNcoeffs(void) const
Definition: SegExp.cpp:788
virtual const LibUtilities::BasisSharedPtr & v_GetBasis(int dir) const
Definition: SegExp.cpp:793
virtual SpatialDomains::GeomType v_MetricInfoType()
Definition: SegExp.cpp:778
virtual void v_ComputeTraceNormal(const int vertex) override
Definition: SegExp.cpp:854
virtual const Array< OneD, const NekDouble > & v_GetPhysNormals(void) override
Definition: SegExp.cpp:772
virtual void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: SegExp.cpp:536
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)
Evaluate the derivative at the physical quadrature points given by inarray and return in outarray.
const LibUtilities::BasisSharedPtr & GetBasis(int dir) const
This function gets the shared point to basis in the dir direction.
Definition: StdExpansion.h:118
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
const LibUtilities::PointsKeyVector GetPointsKeys() const
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:760
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:687
void NormVectorIProductWRTBase(const Array< OneD, const NekDouble > &Fx, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:621
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
DNekMatSharedPtr GenMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:845
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:213
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:226
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
const ConstFactorMap & GetConstFactors() const
Definition: StdMatrixKey.h:140
LibUtilities::ShapeType GetShapeType() const
Definition: StdMatrixKey.h:90
const VarCoeffMap & GetVarCoeffs() const
Definition: StdMatrixKey.h:167
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:85
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:126
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 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
@ ePolyEvenlySpaced
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:75
@ eGaussGaussLegendre
1D Gauss-Gauss-Legendre quadrature points
Definition: PointsType.h:48
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:51
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:59
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:58
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:50
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.
@ eNoGeomType
No type defined.
@ eMovingRegular
Currently unused.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
Definition: Geometry.h:63
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:282
std::shared_ptr< StdSegExp > StdSegExpSharedPtr
Definition: StdSegExp.h:46
The above copyright notice and this permission notice shall be included.
Definition: CoupledSolver.h:1
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
static Array< OneD, Array< OneD, NekDouble > > NullNekDoubleArrayOfArray
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 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 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 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