Nektar++
StdSegExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: StdSegExp.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: Routines within Standard Segment Expansions
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
38using namespace std;
39
40namespace Nektar::StdRegions
41{
42/** \brief Constructor using BasisKey class for quadrature points and
43 * order definition
44 *
45 * \param Ba BasisKey class definition containing order and quadrature
46 * points.
47 */
49 : StdExpansion(Ba.GetNumModes(), 1, Ba),
50 StdExpansion1D(Ba.GetNumModes(), Ba)
51{
52}
53
54/** \brief Return Shape of region, using ShapeType enum list.
55 * i.e. Segment
56 */
58{
60}
61
63{
64
65 bool returnval = false;
66
68 {
69 returnval = true;
70 }
71
73 {
74 returnval = true;
75 }
76
77 return returnval;
78}
79
80//---------------------------------------------------------------------
81// Integration Methods
82//---------------------------------------------------------------------
83
84/** \brief Integrate the physical point list \a inarray over region
85 * and return the value
86 *
87 * \param inarray definition of function to be integrated evauluated at
88 * quadrature point of expansion.
89 * \return returns \f$\int^1_{-1} u(\xi_1)d \xi_1 \f$ where \f$inarray[i]
90 * = u(\xi_{1i}) \f$
91 */
93{
94 NekDouble Int = 0.0;
95 int nquad0 = m_base[0]->GetNumPoints();
96 Array<OneD, NekDouble> tmp(nquad0);
99
100 // multiply by integration constants
101 Vmath::Vmul(nquad0, inarray, 1, w0, 1, tmp, 1);
102
103 Int = Vmath::Vsum(nquad0, tmp, 1);
104
105 return Int;
106}
107
108//---------------------------------------------------------------------
109// Differentiation Methods
110//---------------------------------------------------------------------
111
112/** \brief Evaluate the derivative \f$ d/d{\xi_1} \f$ at the physical
113 * quadrature points given by \a inarray and return in \a outarray.
114 *
115 * This is a wrapper around StdExpansion1D::Tensor_Deriv
116 * \param inarray array of a function evaluated at the quadrature points
117 * \param outarray the resulting array of the derivative \f$
118 * du/d_{\xi_1}|_{\xi_{1i}} \f$ will be stored in the array \a outarra
119 */
120
123 [[maybe_unused]] Array<OneD, NekDouble> &out_d1,
124 [[maybe_unused]] Array<OneD, NekDouble> &out_d2)
125{
126 PhysTensorDeriv(inarray, out_d0);
127}
128
129void StdSegExp::v_PhysDeriv([[maybe_unused]] const int dir,
130 const Array<OneD, const NekDouble> &inarray,
131 Array<OneD, NekDouble> &outarray)
132{
133 ASSERTL1(dir == 0, "input dir is out of range");
134 PhysTensorDeriv(inarray, outarray);
135 // PhysDeriv(inarray, outarray);
136}
137
140 [[maybe_unused]] Array<OneD, NekDouble> &out_d1,
141 [[maybe_unused]] Array<OneD, NekDouble> &out_d2)
142{
143 PhysTensorDeriv(inarray, out_d0);
144}
145
146void StdSegExp::v_StdPhysDeriv([[maybe_unused]] const int dir,
147 const Array<OneD, const NekDouble> &inarray,
148 Array<OneD, NekDouble> &outarray)
149{
150 ASSERTL1(dir == 0, "input dir is out of range");
151 PhysTensorDeriv(inarray, outarray);
152}
153
154//---------------------------------------------------------------------
155// Transforms
156//---------------------------------------------------------------------
157
158/** \brief Backward transform from coefficient space given
159 * in \a inarray and evaluate at the physical quadrature
160 * points \a outarray
161 *
162 * Operation can be evaluated as \f$ u(\xi_{1i}) =
163 * \sum_{p=0}^{order-1} \hat{u}_p \phi_p(\xi_{1i}) \f$ or equivalently
164 * \f$ {\bf u} = {\bf B}^T {\bf \hat{u}} \f$ where
165 * \f${\bf B}[i][j] = \phi_i(\xi_{1j}), \mbox{\_coeffs}[p] = {\bf
166 * \hat{u}}[p] \f$
167 *
168 * The function takes the coefficient array \a inarray as
169 * input for the transformation
170 *
171 * \param inarray: the coeffficients of the expansion
172 *
173 * \param outarray: the resulting array of the values of the function at
174 * the physical quadrature points will be stored in the array \a outarray
175 */
176
178 Array<OneD, NekDouble> &outarray)
179{
180 int nquad = m_base[0]->GetNumPoints();
181
182 if (m_base[0]->Collocation())
183 {
184 Vmath::Vcopy(nquad, inarray, 1, outarray, 1);
185 }
186 else
187 {
188 Blas::Dgemv('N', nquad, m_base[0]->GetNumModes(), 1.0,
189 (m_base[0]->GetBdata()).data(), nquad, &inarray[0], 1, 0.0,
190 &outarray[0], 1);
191 }
192}
193
195 const Array<OneD, const NekDouble> &inarray,
196 Array<OneD, NekDouble> &outarray)
197{
198 int n_coeffs = inarray.size();
199
200 Array<OneD, NekDouble> coeff(n_coeffs);
201 Array<OneD, NekDouble> coeff_tmp(n_coeffs, 0.0);
204
205 int nmodes0 = m_base[0]->GetNumModes();
206
207 Vmath::Vcopy(n_coeffs, inarray, 1, coeff_tmp, 1);
208
209 const LibUtilities::PointsKey Pkey0(nmodes0,
211
212 LibUtilities::BasisKey b0(m_base[0]->GetBasisType(), nmodes0, Pkey0);
213
214 LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A, nmodes0, Pkey0);
215
216 LibUtilities::InterpCoeff1D(b0, coeff_tmp, bortho0, coeff);
217
218 Vmath::Zero(n_coeffs, coeff_tmp, 1);
219
220 Vmath::Vcopy(numMin, tmp = coeff, 1, tmp2 = coeff_tmp, 1);
221
222 LibUtilities::InterpCoeff1D(bortho0, coeff_tmp, b0, outarray);
223}
224
225/**
226 * \brief Forward transform from physical quadrature space stored in \a
227 * inarray and evaluate the expansion coefficients and store in \a
228 * outarray
229 *
230 * Perform a forward transform using a Galerkin projection by taking the
231 * inner product of the physical points and multiplying by the inverse
232 * of the mass matrix using the Solve method of the standard matrix
233 * container holding the local mass matrix, i.e. \f$ {\bf \hat{u}} =
234 * {\bf M}^{-1} {\bf I} \f$ where \f$ {\bf I}[p] = \int^1_{-1}
235 * \phi_p(\xi_1) u(\xi_1) d\xi_1 \f$
236 *
237 * This function stores the expansion coefficients calculated by the
238 * transformation in the coefficient space array \a outarray
239 *
240 * \param inarray: array of physical quadrature points to be transformed
241 * \param outarray: the coeffficients of the expansion
242 */
244 Array<OneD, NekDouble> &outarray)
245{
246 if (m_base[0]->Collocation())
247 {
248 Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
249 }
250 else
251 {
252 v_IProductWRTBase(inarray, outarray);
253
254 // get Mass matrix inverse
255 StdMatrixKey masskey(eInvMass, v_DetShapeType(), *this);
256 DNekMatSharedPtr matsys = GetStdMatrix(masskey);
257
258 NekVector<NekDouble> in(m_ncoeffs, outarray, eCopy);
260
261 out = (*matsys) * in;
262 }
263}
264
266 const Array<OneD, const NekDouble> &inarray,
267 Array<OneD, NekDouble> &outarray)
268{
269 if (m_base[0]->Collocation())
270 {
271 Vmath::Vcopy(m_ncoeffs, inarray, 1, outarray, 1);
272 }
273 else
274 {
275 int nInteriorDofs = m_ncoeffs - 2;
276 int offset = 0;
277
278 switch (m_base[0]->GetBasisType())
279 {
281 {
282 offset = 1;
283 }
284 break;
286 {
287 nInteriorDofs = m_ncoeffs;
288 offset = 0;
289 }
290 break;
293 {
294 offset = 2;
295 }
296 break;
297 default:
298 ASSERTL0(false, "This type of FwdTrans is not defined for this "
299 "expansion type");
300 }
301
302 fill(outarray.data(), outarray.data() + m_ncoeffs, 0.0);
303
305 {
306 outarray[GetVertexMap(0)] = inarray[0];
307 outarray[GetVertexMap(1)] = inarray[m_base[0]->GetNumPoints() - 1];
308
309 if (m_ncoeffs > 2)
310 {
311 // ideally, we would like to have tmp0 to be replaced by
312 // outarray (currently MassMatrixOp does not allow aliasing)
315
316 StdMatrixKey masskey(eMass, v_DetShapeType(), *this);
317 MassMatrixOp(outarray, tmp0, masskey);
318 v_IProductWRTBase(inarray, tmp1);
319
320 Vmath::Vsub(m_ncoeffs, tmp1, 1, tmp0, 1, tmp1, 1);
321
322 // get Mass matrix inverse (only of interior DOF)
323 DNekMatSharedPtr matsys =
324 (m_stdStaticCondMatrixManager[masskey])->GetBlock(1, 1);
325
326 Blas::Dgemv('N', nInteriorDofs, nInteriorDofs, 1.0,
327 &(matsys->GetPtr())[0], nInteriorDofs,
328 tmp1.data() + offset, 1, 0.0,
329 outarray.data() + offset, 1);
330 }
331 }
332 else
333 {
334 StdSegExp::v_FwdTrans(inarray, outarray);
335 }
336 }
337}
338
340 Array<OneD, NekDouble> &outarray)
341{
342 v_BwdTrans(inarray, outarray);
343}
344
345//---------------------------------------------------------------------
346// Inner product functions
347//---------------------------------------------------------------------
348
349/** \brief Inner product of \a inarray over region with respect to
350 * expansion basis \a base and return in \a outarray
351 *
352 * Calculate \f$ I[p] = \int^{1}_{-1} \phi_p(\xi_1) u(\xi_1) d\xi_1
353 * = \sum_{i=0}^{nq-1} \phi_p(\xi_{1i}) u(\xi_{1i}) w_i \f$ where
354 * \f$ outarray[p] = I[p], inarray[i] = u(\xi_{1i}), base[p*nq+i] =
355 * \phi_p(\xi_{1i}) \f$.
356 *
357 * \param base an array defining the local basis for the inner product
358 * usually passed from Basis->GetBdata() or Basis->GetDbdata()
359 * \param inarray: physical point array of function to be integrated
360 * \f$ u(\xi_1) \f$
361 * \param coll_check flag to identify when a Basis->Collocation() call
362 * should be performed to see if this is a GLL_Lagrange basis with a
363 * collocation property. (should be set to 0 if taking the inner
364 * product with respect to the derivative of basis)
365 * \param outarray the values of the inner product with respect to
366 * each basis over region will be stored in the array \a outarray as
367 * output of the function
368 */
369
371 const Array<OneD, const NekDouble> &inarray,
372 Array<OneD, NekDouble> &outarray,
373 [[maybe_unused]] int coll_check)
374{
375 int nquad = m_base[0]->GetNumPoints();
376 Array<OneD, NekDouble> tmp(nquad);
378
379 Vmath::Vmul(nquad, inarray, 1, w, 1, tmp, 1);
380
381 /* Comment below was a bug for collocated basis
382 if(coll_check&&m_base[0]->Collocation())
383 {
384 Vmath::Vcopy(nquad, tmp, 1, outarray, 1);
385 }
386 else
387 {
388 Blas::Dgemv('T',nquad,m_ncoeffs,1.0,base.data(),nquad,
389 &tmp[0],1,0.0,outarray.data(),1);
390 }*/
391
392 // Correct implementation
393 Blas::Dgemv('T', nquad, m_ncoeffs, 1.0, base.data(), nquad, &tmp[0], 1, 0.0,
394 outarray.data(), 1);
395}
396
397/** \brief Inner product of \a inarray over region with respect to the
398 * expansion basis (this)->m_base[0] and return in \a outarray
399 *
400 * Wrapper call to StdSegExp::IProductWRTBase
401 * \param inarray array of function values evaluated at the physical
402 * collocation points
403 * \param outarray the values of the inner product with respect to
404 * each basis over region will be stored in the array \a outarray as
405 * output of the function
406 */
408 Array<OneD, NekDouble> &outarray)
409{
410 v_IProductWRTBase(m_base[0]->GetBdata(), inarray, outarray, 1);
411}
412
414 const int dir, const Array<OneD, const NekDouble> &inarray,
415 Array<OneD, NekDouble> &outarray)
416{
417 StdSegExp::IProductWRTDerivBase_SumFac(dir, inarray, outarray);
418}
419
421 [[maybe_unused]] const int dir, const Array<OneD, const NekDouble> &inarray,
422 Array<OneD, NekDouble> &outarray)
423{
424 ASSERTL1(dir == 0, "input dir is out of range");
425 StdSegExp::v_IProductWRTBase(m_base[0]->GetDbdata(), inarray, outarray, 1);
426}
427
429 const Array<OneD, const NekDouble> &inarray,
430 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
431{
432 int nquad = m_base[0]->GetNumPoints();
433 Array<OneD, NekDouble> tmp(nquad);
435 Array<OneD, const NekDouble> base = m_base[0]->GetBdata();
436
437 if (multiplybyweights)
438 {
439 Vmath::Vmul(nquad, inarray, 1, w, 1, tmp, 1);
440
441 Blas::Dgemv('T', nquad, m_ncoeffs, 1.0, base.data(), nquad, &tmp[0], 1,
442 0.0, outarray.data(), 1);
443 }
444 else
445 {
446 Blas::Dgemv('T', nquad, m_ncoeffs, 1.0, base.data(), nquad, &inarray[0],
447 1, 0.0, outarray.data(), 1);
448 }
449}
450
451//----------------------------
452// Evaluation
453//----------------------------
456{
457 eta[0] = xi[0];
458}
459
462{
463 xi[0] = eta[0];
464}
465
466void StdSegExp::v_FillMode(const int mode, Array<OneD, NekDouble> &outarray)
467{
468 int nquad = m_base[0]->GetNumPoints();
469 const NekDouble *base = m_base[0]->GetBdata().data();
470
471 ASSERTL2(mode <= m_ncoeffs,
472 "calling argument mode is larger than total expansion order");
473
474 Vmath::Vcopy(nquad, (NekDouble *)base + mode * nquad, 1, &outarray[0], 1);
475}
476
478 const Array<OneD, const NekDouble> &coords, int mode)
479{
480 return StdExpansion::BaryEvaluateBasis<0>(coords[0], mode);
481}
482
484 const Array<OneD, NekDouble> &coord,
485 const Array<OneD, const NekDouble> &inarray,
486 std::array<NekDouble, 3> &firstOrderDerivs)
487{
488 return StdExpansion1D::BaryTensorDeriv(coord, inarray, firstOrderDerivs);
489}
490
492 const Array<OneD, NekDouble> &coord,
493 const Array<OneD, const NekDouble> &inarray,
494 std::array<NekDouble, 3> &firstOrderDerivs,
495 std::array<NekDouble, 6> &secondOrderDerivs)
496{
497 return StdExpansion1D::BaryTensorDeriv(coord, inarray, firstOrderDerivs,
498 secondOrderDerivs);
499}
501 Array<OneD, NekDouble> &outarray,
502 [[maybe_unused]] const StdMatrixKey &mkey)
503{
504 int nquad = m_base[0]->GetNumPoints();
505
506 Array<OneD, NekDouble> physValues(nquad);
507 Array<OneD, NekDouble> dPhysValuesdx(nquad);
508
509 v_BwdTrans(inarray, physValues);
510
511 // Laplacian matrix operation
512 v_PhysDeriv(physValues, dPhysValuesdx);
513 v_IProductWRTBase(m_base[0]->GetDbdata(), dPhysValuesdx, outarray, 1);
514}
515
516void StdSegExp::v_LaplacianMatrixOp(const int k1, const int k2,
517 const Array<OneD, const NekDouble> &inarray,
518 Array<OneD, NekDouble> &outarray,
519 const StdMatrixKey &mkey)
520{
521 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
522}
523
525 Array<OneD, NekDouble> &outarray,
526 const StdMatrixKey &mkey)
527{
528 int nquad = m_base[0]->GetNumPoints();
529
530 Array<OneD, NekDouble> physValues(nquad);
531 Array<OneD, NekDouble> dPhysValuesdx(nquad);
533
534 v_BwdTrans(inarray, physValues);
535
536 // mass matrix operation
537 v_IProductWRTBase((m_base[0]->GetBdata()), physValues, wsp, 1);
538
539 // Laplacian matrix operation
540 v_PhysDeriv(physValues, dPhysValuesdx);
541 v_IProductWRTBase(m_base[0]->GetDbdata(), dPhysValuesdx, outarray, 1);
543 outarray.data(), 1);
544}
545
547 const StdMatrixKey &mkey)
548{
549 // Generate an orthogonal expansion
550 int nq = m_base[0]->GetNumPoints();
551 int nmodes = m_base[0]->GetNumModes();
552 // Declare orthogonal basis.
554
556 StdSegExp OrthoExp(B);
557
558 // SVV parameters loaded from the .xml case file
559 NekDouble SvvDiffCoeff = mkey.GetConstFactor(eFactorSVVDiffCoeff);
560 int cutoff = (int)(mkey.GetConstFactor(eFactorSVVCutoffRatio)) * nmodes;
561
562 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
563
564 // project onto modal space.
565 OrthoExp.FwdTrans(array, orthocoeffs);
566
567 //
568 for (int j = 0; j < nmodes; ++j)
569 {
570 if (j >= cutoff) // to filter out only the "high-modes"
571 {
572 orthocoeffs[j] *=
573 (SvvDiffCoeff *
574 exp(-(j - nmodes) * (j - nmodes) /
575 ((NekDouble)((j - cutoff + 1) * (j - cutoff + 1)))));
576 }
577 else
578 {
579 orthocoeffs[j] *= 0.0;
580 }
581 }
582
583 // backward transform to physical space
584 OrthoExp.BwdTrans(orthocoeffs, array);
585}
586
588 const NekDouble alpha,
589 const NekDouble exponent,
590 const NekDouble cutoff)
591{
592 // Generate an orthogonal expansion
593 int nq = m_base[0]->GetNumPoints();
594 int nmodes = m_base[0]->GetNumModes();
595 int P = nmodes - 1;
596 // Declare orthogonal basis.
598
600 StdSegExp OrthoExp(B);
601
602 // Cutoff
603 int Pcut = cutoff * P;
604
605 // Project onto orthogonal space.
606 Array<OneD, NekDouble> orthocoeffs(OrthoExp.GetNcoeffs());
607 OrthoExp.FwdTrans(array, orthocoeffs);
608
609 //
610 NekDouble fac;
611 for (int j = 0; j < nmodes; ++j)
612 {
613 // to filter out only the "high-modes"
614 if (j > Pcut)
615 {
616 fac = (NekDouble)(j - Pcut) / ((NekDouble)(P - Pcut));
617 fac = pow(fac, exponent);
618 orthocoeffs[j] *= exp(-alpha * fac);
619 }
620 }
621
622 // backward transform to physical space
623 OrthoExp.BwdTrans(orthocoeffs, array);
624}
625
626// up to here
628 const Array<OneD, const NekDouble> &inarray,
629 Array<OneD, NekDouble> &outarray)
630{
631 int nquad0 = m_base[0]->GetNumPoints();
632
633 const Array<OneD, const NekDouble> &w0 = m_base[0]->GetW();
634
635 Vmath::Vmul(nquad0, inarray.data(), 1, w0.data(), 1, outarray.data(), 1);
636}
637
639 [[maybe_unused]] Array<OneD, NekDouble> &coords_1,
640 [[maybe_unused]] Array<OneD, NekDouble> &coords_2)
641{
642 Blas::Dcopy(GetNumPoints(0), (m_base[0]->GetZ()).data(), 1, &coords_0[0],
643 1);
644}
645
646//---------------------------------------------------------------------
647// Helper functions
648//---------------------------------------------------------------------
649
651{
652 return 2;
653}
654
656{
657 return 2;
658}
659
660int StdSegExp::v_GetTraceNcoeffs([[maybe_unused]] const int i) const
661{
662 return 1;
663}
664
665int StdSegExp::v_GetTraceIntNcoeffs([[maybe_unused]] const int i) const
666{
667 return 0;
668}
669
670int StdSegExp::v_GetTraceNumPoints([[maybe_unused]] const int i) const
671{
672 return 1;
673}
674
676{
677 return 2;
678}
679
681{
682 return 2;
683}
684
686 const std::vector<unsigned int> &nummodes, int &modes_offset)
687{
688 int nmodes = nummodes[modes_offset];
689 modes_offset += 1;
690
691 return nmodes;
692}
693
694//---------------------------------------------------------------------
695// Wrapper functions
696//---------------------------------------------------------------------
697
699{
701 MatrixType mattype;
702
703 switch (mattype = mkey.GetMatrixType())
704 {
706 {
707 int nq = m_base[0]->GetNumPoints();
708
709 // take definition from key
711 {
712 nq = (int)mkey.GetConstFactor(eFactorConst);
713 }
714
716 Array<OneD, NekDouble> coords(1);
719
720 for (int i = 0; i < neq; ++i)
721 {
722 coords[0] = -1.0 + 2 * i / (NekDouble)(neq - 1);
723 I = m_base[0]->GetI(coords);
724 Vmath::Vcopy(nq, I->GetRawPtr(), 1, Mat->GetRawPtr() + i, neq);
725 }
726 }
727 break;
728 case eFwdTrans:
729 {
730 Mat =
732 StdMatrixKey iprodkey(eIProductWRTBase, v_DetShapeType(), *this);
733 DNekMat &Iprod = *GetStdMatrix(iprodkey);
734 StdMatrixKey imasskey(eInvMass, v_DetShapeType(), *this);
735 DNekMat &Imass = *GetStdMatrix(imasskey);
736
737 (*Mat) = Imass * Iprod;
738 }
739 break;
740 default:
741 {
743
744 if (mattype == eMass)
745 {
746 // For Fourier basis set the imaginary component
747 // of mean mode to have a unit diagonal component
748 // in mass matrix
750 {
751 (*Mat)(1, 1) = 1.0;
752 }
753 }
754 }
755 break;
756 }
757
758 return Mat;
759}
760
762{
763 return v_GenMatrix(mkey);
764}
765
766//---------------------------------------------------------------------
767// Mappings
768//---------------------------------------------------------------------
769
771{
772 if (outarray.size() != NumBndryCoeffs())
773 {
775 }
776 const LibUtilities::BasisType Btype = GetBasisType(0);
777 int nummodes = m_base[0]->GetNumModes();
778
779 outarray[0] = 0;
780
781 switch (Btype)
782 {
787 outarray[1] = nummodes - 1;
788 break;
791 outarray[1] = 1;
792 break;
793 default:
794 ASSERTL0(0, "Mapping array is not defined for this expansion");
795 break;
796 }
797}
798
800{
801 int i;
802 if (outarray.size() != GetNcoeffs() - NumBndryCoeffs())
803 {
805 }
806 const LibUtilities::BasisType Btype = GetBasisType(0);
807
808 switch (Btype)
809 {
814 for (i = 0; i < GetNcoeffs() - 2; i++)
815 {
816 outarray[i] = i + 1;
817 }
818 break;
821 for (i = 0; i < GetNcoeffs() - 2; i++)
822 {
823 outarray[i] = i + 2;
824 }
825 break;
826 default:
827 ASSERTL0(0, "Mapping array is not defined for this expansion");
828 break;
829 }
830}
831
832int StdSegExp::v_GetVertexMap(int localVertexId,
833 [[maybe_unused]] bool useCoeffPacking)
834{
835 ASSERTL0((localVertexId == 0) || (localVertexId == 1),
836 "local vertex id"
837 "must be between 0 or 1");
838
839 int localDOF = localVertexId;
840
842 (localVertexId == 1))
843 {
844 localDOF = m_base[0]->GetNumModes() - 1;
845 }
846 return localDOF;
847}
848
850 Array<OneD, int> &conn, [[maybe_unused]] bool standard)
851{
852 int np = m_base[0]->GetNumPoints();
853
854 conn = Array<OneD, int>(2 * (np - 1));
855 int cnt = 0;
856 for (int i = 0; i < np - 1; ++i)
857 {
858 conn[cnt++] = i;
859 conn[cnt++] = i + 1;
860 }
861}
862
863/** \brief Get the map of the coefficient location to teh
864 * local trace coefficients
865 */
866
867void StdSegExp::v_GetTraceCoeffMap(const unsigned int traceid,
869{
870 int order0 = m_base[0]->GetNumModes();
871
872 ASSERTL0(traceid < 2, "eid must be between 0 and 1");
873
874 if (maparray.size() != 1)
875 {
876 maparray = Array<OneD, unsigned int>(1);
877 }
878
879 const LibUtilities::BasisType bType = GetBasisType(0);
880
881 if (bType == LibUtilities::eModified_A)
882 {
883 maparray[0] = (traceid == 0) ? 0 : 1;
884 }
885 else if (bType == LibUtilities::eGLL_Lagrange ||
887 {
888 maparray[0] = (traceid == 0) ? 0 : order0 - 1;
889 }
890 else
891 {
892 ASSERTL0(false, "Unknown Basis");
893 }
894}
895
898 Array<OneD, int> &signarray,
899 [[maybe_unused]] Orientation orient,
900 [[maybe_unused]] int P,
901 [[maybe_unused]] int Q)
902{
903 v_GetTraceCoeffMap(tid, maparray);
904
905 if (signarray.size() != 1)
906 {
907 signarray = Array<OneD, int>(1, 1);
908 }
909 else
910 {
911 signarray[0] = 1;
912 }
913}
914
916 [[maybe_unused]] const unsigned int eid,
917 Array<OneD, unsigned int> &maparray, Array<OneD, int> &signarray,
918 [[maybe_unused]] Orientation orient, [[maybe_unused]] int P,
919 [[maybe_unused]] int Q)
920{
921 // parameters for higher dimnesion traces
922 if (maparray.size() != 1)
923 {
924 maparray = Array<OneD, unsigned int>(1);
925 }
926
927 maparray[0] = 0;
928
929 if (signarray.size() != 1)
930 {
931 signarray = Array<OneD, int>(1, 1);
932 }
933 else
934 {
935 signarray[0] = 1;
936 }
937}
938
939} // namespace Nektar::StdRegions
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
#define ASSERTL2(condition, msg)
Assert Level 2 – Debugging which is used FULLDEBUG compilation mode. This level assert is designed to...
Definition: ErrorUtil.hpp:265
Describes the specification for a Basis.
Definition: Basis.h:45
Defines a specification for a set of points.
Definition: Points.h:50
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.
NekDouble BaryTensorDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs)
The base class for all shapes.
Definition: StdExpansion.h:65
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:156
DNekMatSharedPtr GetStdMatrix(const StdMatrixKey &mkey)
Definition: StdExpansion.h:612
void MassMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
Definition: StdExpansion.h:761
void IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
int GetVertexMap(const int localVertexId, bool useCoeffPacking=false)
Definition: StdExpansion.h:688
DNekMatSharedPtr CreateGeneralMatrix(const StdMatrixKey &mkey)
this function generates the mass matrix
void LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey)
void 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:433
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:205
void FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function performs the Forward transformation from physical space to coefficient space.
LibUtilities::NekManager< StdMatrixKey, DNekBlkMat, StdMatrixKey::opLess > m_stdStaticCondMatrixManager
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:218
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:83
NekDouble GetConstFactor(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:124
bool ConstFactorExists(const ConstFactorType &factor) const
Definition: StdMatrixKey.h:133
Class representing a segment element in reference space All interface of this class sits in StdExpans...
Definition: StdSegExp.h:45
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)->m_base[0] and return...
Definition: StdSegExp.cpp:407
int v_NumDGBndryCoeffs() const override
Definition: StdSegExp.cpp:680
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: StdSegExp.cpp:483
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdSegExp.cpp:500
void v_ReduceOrderCoeffs(int numMin, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdSegExp.cpp:194
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdMatrixKey &mkey) override
Definition: StdSegExp.cpp:546
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Definition: StdSegExp.cpp:428
void v_GetInteriorMap(Array< OneD, unsigned int > &outarray) override
Definition: StdSegExp.cpp:799
NekDouble v_PhysEvaluateBasis(const Array< OneD, const NekDouble > &coords, int mode) final
Definition: StdSegExp.cpp:477
LibUtilities::ShapeType v_DetShapeType() const override
Return Shape of region, using ShapeType enum list. i.e. Segment.
Definition: StdSegExp.cpp:57
StdSegExp(const LibUtilities::BasisKey &Ba)
Constructor using BasisKey class for quadrature points and order definition.
Definition: StdSegExp.cpp:48
NekDouble v_PhysEvalFirstSecondDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs, std::array< NekDouble, 6 > &secondOrderDerivs) override
Definition: StdSegExp.cpp:491
int v_NumBndryCoeffs() const override
Definition: StdSegExp.cpp:675
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdMatrixKey &mkey) override
Definition: StdSegExp.cpp:524
void v_BwdTrans_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdSegExp.cpp:339
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdSegExp.cpp:413
void v_LocCollapsedToLocCoord(const Array< OneD, const NekDouble > &eta, Array< OneD, NekDouble > &xi) override
Definition: StdSegExp.cpp:460
void v_GetTraceCoeffMap(const unsigned int traceid, Array< OneD, unsigned int > &maparray) override
Get the map of the coefficient location to teh local trace coefficients.
Definition: StdSegExp.cpp:867
int v_GetNverts() const final
Definition: StdSegExp.cpp:650
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: StdSegExp.cpp:243
int v_GetTraceNumPoints(const int i) const final
Definition: StdSegExp.cpp:670
void v_FillMode(const int mode, Array< OneD, NekDouble > &outarray) override
Definition: StdSegExp.cpp:466
int v_GetTraceNcoeffs(const int i) const final
Definition: StdSegExp.cpp:660
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: StdSegExp.cpp:121
void v_GetSimplexEquiSpacedConnectivity(Array< OneD, int > &conn, bool standard=true) override
Definition: StdSegExp.cpp:849
void v_GetCoords(Array< OneD, NekDouble > &coords_0, Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2) override
Definition: StdSegExp.cpp:638
int v_GetTraceIntNcoeffs(const int i) const final
Definition: StdSegExp.cpp:665
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrate the physical point list inarray over region and return the value.
Definition: StdSegExp.cpp:92
int v_GetVertexMap(int localVertexId, bool useCoeffPacking=false) override
Definition: StdSegExp.cpp:832
void v_ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff) override
Definition: StdSegExp.cpp:587
DNekMatSharedPtr v_CreateStdMatrix(const StdMatrixKey &mkey) override
Definition: StdSegExp.cpp:761
void v_IProductWRTDerivBase_SumFac(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdSegExp.cpp:420
void v_MultiplyByStdQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdSegExp.cpp:627
void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Backward transform from coefficient space given in inarray and evaluate at the physical quadrature po...
Definition: StdSegExp.cpp:177
int v_GetNtraces() const final
Definition: StdSegExp.cpp:655
void v_GetBoundaryMap(Array< OneD, unsigned int > &outarray) override
Definition: StdSegExp.cpp:770
void v_GetTraceToElementMap(const int tid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation edgeOrient, int P, int Q) override
Definition: StdSegExp.cpp:896
void v_GetElmtTraceToTraceMap(const unsigned int eid, Array< OneD, unsigned int > &maparray, Array< OneD, int > &signarray, Orientation edgeOrient, int P, int Q) override
Definition: StdSegExp.cpp:915
void v_StdPhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1=NullNekDouble1DArray, Array< OneD, NekDouble > &out_d2=NullNekDouble1DArray) override
Definition: StdSegExp.cpp:138
void v_LocCoordToLocCollapsed(const Array< OneD, const NekDouble > &xi, Array< OneD, NekDouble > &eta) override
Definition: StdSegExp.cpp:454
bool v_IsBoundaryInteriorExpansion() const override
Definition: StdSegExp.cpp:62
int v_CalcNumberOfCoefficients(const std::vector< unsigned int > &nummodes, int &modes_offset) override
Definition: StdSegExp.cpp:685
DNekMatSharedPtr v_GenMatrix(const StdMatrixKey &mkey) override
Definition: StdSegExp.cpp:698
void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Definition: StdSegExp.cpp:265
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 = alpha A x plus beta y where A[m x n].
Definition: Blas.hpp:211
static void Dcopy(const int &n, const double *x, const int &incx, double *y, const int &incy)
BLAS level 1: Copy x to y.
Definition: Blas.hpp:128
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:135
void InterpCoeff1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
Definition: InterpCoeff.cpp:42
@ eGaussLobattoLegendre
1D Gauss-Lobatto-Legendre quadrature points
Definition: PointsType.h:51
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
@ P
Monomial polynomials .
Definition: BasisType.h:62
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
Definition: BasisType.h:57
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eGLL_Lagrange
Lagrange for SEM basis .
Definition: BasisType.h:56
@ eChebyshev
Chebyshev Polynomials.
Definition: BasisType.h:61
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
@ eFourier
Fourier Expansion .
Definition: BasisType.h:55
std::vector< double > w(NPUPPER)
std::vector< double > z(NPUPPER)
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.hpp:72
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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.hpp:220
STL namespace.