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