Nektar++
TetExp.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: TetExp.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:
32//
33///////////////////////////////////////////////////////////////////////////////
34
37#include <LocalRegions/TetExp.h>
39
40using namespace std;
41
43{
44/**
45 * @class TetExp
46 * Defines a Tetrahedral local expansion.
47 */
48
49/**
50 * \brief Constructor using BasisKey class for quadrature points and
51 * order definition
52 *
53 * @param Ba Basis key for first coordinate.
54 * @param Bb Basis key for second coordinate.
55 * @param Bc Basis key for third coordinate.
56 */
58 const LibUtilities::BasisKey &Bb,
59 const LibUtilities::BasisKey &Bc,
61 : StdExpansion(LibUtilities::StdTetData::getNumberOfCoefficients(
62 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
63 3, Ba, Bb, Bc),
64 StdExpansion3D(LibUtilities::StdTetData::getNumberOfCoefficients(
65 Ba.GetNumModes(), Bb.GetNumModes(), Bc.GetNumModes()),
66 Ba, Bb, Bc),
67 StdTetExp(Ba, Bb, Bc), Expansion(geom), Expansion3D(geom),
68 m_matrixManager(
69 std::bind(&Expansion3D::CreateMatrix, this, std::placeholders::_1),
70 std::string("TetExpMatrix")),
71 m_staticCondMatrixManager(std::bind(&Expansion::CreateStaticCondMatrix,
72 this, std::placeholders::_1),
73 std::string("TetExpStaticCondMatrix"))
74{
75}
76
77/**
78 * \brief Copy Constructor
79 */
81 : StdRegions::StdExpansion(T), StdRegions::StdExpansion3D(T),
82 StdRegions::StdTetExp(T), Expansion(T), Expansion3D(T),
83 m_matrixManager(T.m_matrixManager),
84 m_staticCondMatrixManager(T.m_staticCondMatrixManager)
85{
86}
87
88//-----------------------------
89// Integration Methods
90//-----------------------------
91/**
92 * \brief Integrate the physical point list \a inarray over region
93 *
94 * @param inarray Definition of function to be returned at
95 * quadrature point of expansion.
96 * @returns \f$\int^1_{-1}\int^1_{-1} \int^1_{-1}
97 * u(\eta_1, \eta_2, \eta_3) J[i,j,k] d \eta_1 d \eta_2 d \eta_3 \f$
98 * where \f$inarray[i,j,k] = u(\eta_{1i},\eta_{2j},\eta_{3k})
99 * \f$ and \f$ J[i,j,k] \f$ is the Jacobian evaluated at the quadrature
100 * point.
101 */
103{
104 int nquad0 = m_base[0]->GetNumPoints();
105 int nquad1 = m_base[1]->GetNumPoints();
106 int nquad2 = m_base[2]->GetNumPoints();
108 NekDouble retrunVal;
109 Array<OneD, NekDouble> tmp(nquad0 * nquad1 * nquad2);
110
111 // multiply inarray with Jacobian
112 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
113 {
114 Vmath::Vmul(nquad0 * nquad1 * nquad2, &jac[0], 1,
115 (NekDouble *)&inarray[0], 1, &tmp[0], 1);
116 }
117 else
118 {
119 Vmath::Smul(nquad0 * nquad1 * nquad2, (NekDouble)jac[0],
120 (NekDouble *)&inarray[0], 1, &tmp[0], 1);
121 }
122
123 // call StdTetExp version;
124 retrunVal = StdTetExp::v_Integral(tmp);
125
126 return retrunVal;
127}
128
129//-----------------------------
130// Differentiation Methods
131//-----------------------------
132/**
133 * \brief Differentiate \a inarray in the three coordinate directions.
134 *
135 * @param inarray Input array of values at quadrature points to
136 * be differentiated.
137 * @param out_d0 Derivative in first coordinate direction.
138 * @param out_d1 Derivative in second coordinate direction.
139 * @param out_d2 Derivative in third coordinate direction.
140 */
145{
146 int TotPts = m_base[0]->GetNumPoints() * m_base[1]->GetNumPoints() *
147 m_base[2]->GetNumPoints();
148
150 m_metricinfo->GetDerivFactors(GetPointsKeys());
152 Array<OneD, NekDouble> Diff1 = Diff0 + TotPts;
153 Array<OneD, NekDouble> Diff2 = Diff1 + TotPts;
154
155 StdTetExp::v_PhysDeriv(inarray, Diff0, Diff1, Diff2);
156
157 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
158 {
159 if (out_d0.size())
160 {
161 Vmath::Vmul(TotPts, &df[0][0], 1, &Diff0[0], 1, &out_d0[0], 1);
162 Vmath::Vvtvp(TotPts, &df[1][0], 1, &Diff1[0], 1, &out_d0[0], 1,
163 &out_d0[0], 1);
164 Vmath::Vvtvp(TotPts, &df[2][0], 1, &Diff2[0], 1, &out_d0[0], 1,
165 &out_d0[0], 1);
166 }
167
168 if (out_d1.size())
169 {
170 Vmath::Vmul(TotPts, &df[3][0], 1, &Diff0[0], 1, &out_d1[0], 1);
171 Vmath::Vvtvp(TotPts, &df[4][0], 1, &Diff1[0], 1, &out_d1[0], 1,
172 &out_d1[0], 1);
173 Vmath::Vvtvp(TotPts, &df[5][0], 1, &Diff2[0], 1, &out_d1[0], 1,
174 &out_d1[0], 1);
175 }
176
177 if (out_d2.size())
178 {
179 Vmath::Vmul(TotPts, &df[6][0], 1, &Diff0[0], 1, &out_d2[0], 1);
180 Vmath::Vvtvp(TotPts, &df[7][0], 1, &Diff1[0], 1, &out_d2[0], 1,
181 &out_d2[0], 1);
182 Vmath::Vvtvp(TotPts, &df[8][0], 1, &Diff2[0], 1, &out_d2[0], 1,
183 &out_d2[0], 1);
184 }
185 }
186 else // regular geometry
187 {
188 if (out_d0.size())
189 {
190 Vmath::Smul(TotPts, df[0][0], &Diff0[0], 1, &out_d0[0], 1);
191 Blas::Daxpy(TotPts, df[1][0], &Diff1[0], 1, &out_d0[0], 1);
192 Blas::Daxpy(TotPts, df[2][0], &Diff2[0], 1, &out_d0[0], 1);
193 }
194
195 if (out_d1.size())
196 {
197 Vmath::Smul(TotPts, df[3][0], &Diff0[0], 1, &out_d1[0], 1);
198 Blas::Daxpy(TotPts, df[4][0], &Diff1[0], 1, &out_d1[0], 1);
199 Blas::Daxpy(TotPts, df[5][0], &Diff2[0], 1, &out_d1[0], 1);
200 }
201
202 if (out_d2.size())
203 {
204 Vmath::Smul(TotPts, df[6][0], &Diff0[0], 1, &out_d2[0], 1);
205 Blas::Daxpy(TotPts, df[7][0], &Diff1[0], 1, &out_d2[0], 1);
206 Blas::Daxpy(TotPts, df[8][0], &Diff2[0], 1, &out_d2[0], 1);
207 }
208 }
209}
210
211//-----------------------------
212// Transforms
213//-----------------------------
214/**
215 * \brief Forward transform from physical quadrature space stored in
216 * \a inarray and evaluate the expansion coefficients and store
217 * in \a (this)->_coeffs
218 *
219 * @param inarray Array of physical quadrature points to be
220 * transformed.
221 * @param outarray Array of coefficients to update.
222 */
224 Array<OneD, NekDouble> &outarray)
225{
226 if ((m_base[0]->Collocation()) && (m_base[1]->Collocation()) &&
227 (m_base[2]->Collocation()))
228 {
229 Vmath::Vcopy(GetNcoeffs(), &inarray[0], 1, &outarray[0], 1);
230 }
231 else
232 {
233 IProductWRTBase(inarray, outarray);
234
235 // get Mass matrix inverse
237 DNekScalMatSharedPtr matsys = m_matrixManager[masskey];
238
239 // copy inarray in case inarray == outarray
240 DNekVec in(m_ncoeffs, outarray);
241 DNekVec out(m_ncoeffs, outarray, eWrapper);
242
243 out = (*matsys) * in;
244 }
245}
246
247//-----------------------------
248// Inner product functions
249//-----------------------------
250/**
251 * \brief Calculate the inner product of inarray with respect to the
252 * basis B=m_base0*m_base1*m_base2 and put into outarray:
253 *
254 * \f$ \begin{array}{rcl} I_{pqr} = (\phi_{pqr}, u)_{\delta}
255 * & = & \sum_{i=0}^{nq_0} \sum_{j=0}^{nq_1} \sum_{k=0}^{nq_2}
256 * \psi_{p}^{a} (\eta_{1i}) \psi_{pq}^{b} (\eta_{2j}) \psi_{pqr}^{c}
257 * (\eta_{3k}) w_i w_j w_k u(\eta_{1,i} \eta_{2,j} \eta_{3,k})
258 * J_{i,j,k}\\ & = & \sum_{i=0}^{nq_0} \psi_p^a(\eta_{1,i})
259 * \sum_{j=0}^{nq_1} \psi_{pq}^b(\eta_{2,j}) \sum_{k=0}^{nq_2}
260 * \psi_{pqr}^c u(\eta_{1i},\eta_{2j},\eta_{3k}) J_{i,j,k}
261 * \end{array} \f$ \n
262 * where
263 * \f$ \phi_{pqr} (\xi_1 , \xi_2 , \xi_3)
264 * = \psi_p^a (\eta_1) \psi_{pq}^b (\eta_2) \psi_{pqr}^c (\eta_3) \f$
265 * which can be implemented as \n
266 * \f$f_{pqr} (\xi_{3k})
267 * = \sum_{k=0}^{nq_3} \psi_{pqr}^c u(\eta_{1i},\eta_{2j},\eta_{3k})
268 * J_{i,j,k} = {\bf B_3 U} \f$ \n
269 * \f$ g_{pq} (\xi_{3k})
270 * = \sum_{j=0}^{nq_1} \psi_{pq}^b (\xi_{2j}) f_{pqr} (\xi_{3k})
271 * = {\bf B_2 F} \f$ \n
272 * \f$ (\phi_{pqr}, u)_{\delta}
273 * = \sum_{k=0}^{nq_0} \psi_{p}^a (\xi_{3k}) g_{pq} (\xi_{3k})
274 * = {\bf B_1 G} \f$
275 */
277 Array<OneD, NekDouble> &outarray)
278{
279 v_IProductWRTBase_SumFac(inarray, outarray);
280}
281
283 const Array<OneD, const NekDouble> &inarray,
284 Array<OneD, NekDouble> &outarray, bool multiplybyweights)
285{
286 const int nquad0 = m_base[0]->GetNumPoints();
287 const int nquad1 = m_base[1]->GetNumPoints();
288 const int nquad2 = m_base[2]->GetNumPoints();
289 const int order0 = m_base[0]->GetNumModes();
290 const int order1 = m_base[1]->GetNumModes();
291 Array<OneD, NekDouble> wsp(nquad1 * nquad2 * order0 +
292 nquad2 * order0 * (order1 + 1) / 2);
293
294 if (multiplybyweights)
295 {
296 Array<OneD, NekDouble> tmp(nquad0 * nquad1 * nquad2);
297
298 MultiplyByQuadratureMetric(inarray, tmp);
300 m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
301 tmp, outarray, wsp, true, true, true);
302 }
303 else
304 {
306 m_base[0]->GetBdata(), m_base[1]->GetBdata(), m_base[2]->GetBdata(),
307 inarray, outarray, wsp, true, true, true);
308 }
309}
310
311/**
312 * @brief Calculates the inner product \f$ I_{pqr} = (u,
313 * \partial_{x_i} \phi_{pqr}) \f$.
314 *
315 * The derivative of the basis functions is performed using the chain
316 * rule in order to incorporate the geometric factors. Assuming that
317 * the basis functions are a tensor product
318 * \f$\phi_{pqr}(\eta_1,\eta_2,\eta_3) =
319 * \phi_1(\eta_1)\phi_2(\eta_2)\phi_3(\eta_3)\f$, this yields the
320 * result
321 *
322 * \f[
323 * I_{pqr} = \sum_{j=1}^3 \left(u, \frac{\partial u}{\partial \eta_j}
324 * \frac{\partial \eta_j}{\partial x_i}\right)
325 * \f]
326 *
327 * In the prismatic element, we must also incorporate a second set of
328 * geometric factors which incorporate the collapsed co-ordinate
329 * system, so that
330 *
331 * \f[ \frac{\partial\eta_j}{\partial x_i} = \sum_{k=1}^3
332 * \frac{\partial\eta_j}{\partial\xi_k}\frac{\partial\xi_k}{\partial
333 * x_i} \f]
334 *
335 * These derivatives can be found on p152 of Sherwin & Karniadakis.
336 *
337 * @param dir Direction in which to take the derivative.
338 * @param inarray The function \f$ u \f$.
339 * @param outarray Value of the inner product.
340 */
342 const Array<OneD, const NekDouble> &inarray,
343 Array<OneD, NekDouble> &outarray)
344{
345 const int nquad0 = m_base[0]->GetNumPoints();
346 const int nquad1 = m_base[1]->GetNumPoints();
347 const int nquad2 = m_base[2]->GetNumPoints();
348 const int order0 = m_base[0]->GetNumModes();
349 const int order1 = m_base[1]->GetNumModes();
350 const int nqtot = nquad0 * nquad1 * nquad2;
351
352 Array<OneD, NekDouble> tmp1(nqtot);
353 Array<OneD, NekDouble> tmp2(nqtot);
354 Array<OneD, NekDouble> tmp3(nqtot);
355 Array<OneD, NekDouble> tmp4(nqtot);
357 Array<OneD, NekDouble> wsp(nquad1 * nquad2 * order0 +
358 nquad2 * order0 * (order1 + 1) / 2);
359
360 MultiplyByQuadratureMetric(inarray, tmp1);
361
363 tmp2D[0] = tmp2;
364 tmp2D[1] = tmp3;
365 tmp2D[2] = tmp4;
366
367 TetExp::v_AlignVectorToCollapsedDir(dir, tmp1, tmp2D);
368
369 IProductWRTBase_SumFacKernel(m_base[0]->GetDbdata(), m_base[1]->GetBdata(),
370 m_base[2]->GetBdata(), tmp2, outarray, wsp,
371 false, true, true);
372
373 IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetDbdata(),
374 m_base[2]->GetBdata(), tmp3, tmp6, wsp, true,
375 false, true);
376
377 Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
378
379 IProductWRTBase_SumFacKernel(m_base[0]->GetBdata(), m_base[1]->GetBdata(),
380 m_base[2]->GetDbdata(), tmp4, tmp6, wsp, true,
381 true, false);
382
383 Vmath::Vadd(m_ncoeffs, tmp6, 1, outarray, 1, outarray, 1);
384}
385
387 const int dir, const Array<OneD, const NekDouble> &inarray,
389{
390 int i, j;
391
392 const int nquad0 = m_base[0]->GetNumPoints();
393 const int nquad1 = m_base[1]->GetNumPoints();
394 const int nquad2 = m_base[2]->GetNumPoints();
395 const int nqtot = nquad0 * nquad1 * nquad2;
396
397 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
398 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
399 const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
400
401 Array<OneD, NekDouble> tmp2(nqtot);
402 Array<OneD, NekDouble> tmp3(nqtot);
403
405 m_metricinfo->GetDerivFactors(GetPointsKeys());
406
407 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
408 {
409 Vmath::Vmul(nqtot, &df[3 * dir][0], 1, inarray.data(), 1, tmp2.data(),
410 1);
411 Vmath::Vmul(nqtot, &df[3 * dir + 1][0], 1, inarray.data(), 1,
412 tmp3.data(), 1);
413 Vmath::Vmul(nqtot, &df[3 * dir + 2][0], 1, inarray.data(), 1,
414 outarray[2].data(), 1);
415 }
416 else
417 {
418 Vmath::Smul(nqtot, df[3 * dir][0], inarray.data(), 1, tmp2.data(), 1);
419 Vmath::Smul(nqtot, df[3 * dir + 1][0], inarray.data(), 1, tmp3.data(),
420 1);
421 Vmath::Smul(nqtot, df[3 * dir + 2][0], inarray.data(), 1,
422 outarray[2].data(), 1);
423 }
424
425 NekDouble g0, g1, g1a, g2, g3;
426 int k, cnt;
427
428 for (cnt = 0, k = 0; k < nquad2; ++k)
429 {
430 for (j = 0; j < nquad1; ++j)
431 {
432 g2 = 2.0 / (1.0 - z2[k]);
433 g1 = g2 / (1.0 - z1[j]);
434 g0 = 2.0 * g1;
435 g3 = (1.0 + z1[j]) * g2 * 0.5;
436
437 for (i = 0; i < nquad0; ++i, ++cnt)
438 {
439 g1a = g1 * (1 + z0[i]);
440
441 outarray[0][cnt] =
442 g0 * tmp2[cnt] + g1a * (tmp3[cnt] + outarray[2][cnt]);
443
444 outarray[1][cnt] = g2 * tmp3[cnt] + g3 * outarray[2][cnt];
445 }
446 }
447 }
448}
449
450//-----------------------------
451// Evaluation functions
452//-----------------------------
453
454/**
455 * Given the local cartesian coordinate \a Lcoord evaluate the
456 * value of physvals at this point by calling through to the
457 * StdExpansion method
458 */
460 const Array<OneD, const NekDouble> &Lcoord,
461 const Array<OneD, const NekDouble> &physvals)
462{
463 // Evaluate point in local (eta) coordinates.
464 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
465}
466
467/**
468 * @param coord Physical space coordinate
469 * @returns Evaluation of expansion at given coordinate.
470 */
472 const Array<OneD, const NekDouble> &physvals)
473{
474 ASSERTL0(m_geom, "m_geom not defined");
475
477
478 // Get the local (eta) coordinates of the point
479 m_geom->GetLocCoords(coord, Lcoord);
480
481 // Evaluate point in local (eta) coordinates.
482 return StdExpansion3D::v_PhysEvaluate(Lcoord, physvals);
483}
484
486 const Array<OneD, NekDouble> &coord,
487 const Array<OneD, const NekDouble> &inarray,
488 std::array<NekDouble, 3> &firstOrderDerivs)
489{
490 Array<OneD, NekDouble> Lcoord(3);
491 ASSERTL0(m_geom, "m_geom not defined");
492 m_geom->GetLocCoords(coord, Lcoord);
493 return StdTetExp::v_PhysEvalFirstDeriv(Lcoord, inarray, firstOrderDerivs);
494}
495
496/**
497 * \brief Get the coordinates "coords" at the local coordinates "Lcoords"
498 */
501{
502 int i;
503
504 ASSERTL1(Lcoords[0] <= -1.0 && Lcoords[0] >= 1.0 && Lcoords[1] <= -1.0 &&
505 Lcoords[1] >= 1.0 && Lcoords[2] <= -1.0 && Lcoords[2] >= 1.0,
506 "Local coordinates are not in region [-1,1]");
507
508 // m_geom->FillGeom(); // TODO: implement FillGeom()
509
510 for (i = 0; i < m_geom->GetCoordim(); ++i)
511 {
512 coords[i] = m_geom->GetCoord(i, Lcoords);
513 }
514}
515
517 Array<OneD, NekDouble> &coords_1,
518 Array<OneD, NekDouble> &coords_2)
519{
520 Expansion::v_GetCoords(coords_0, coords_1, coords_2);
521}
522
523//-----------------------------
524// Helper functions
525//-----------------------------
526
527/**
528 * \brief Return Shape of region, using ShapeType enum list.
529 */
531{
533}
534
536{
538 m_base[0]->GetBasisKey(), m_base[1]->GetBasisKey(),
539 m_base[2]->GetBasisKey());
540}
541
543{
545 m_base[0]->GetPointsKey());
547 m_base[1]->GetPointsKey());
549 m_base[2]->GetPointsKey());
550
552 bkey2);
553}
554
556 const NekDouble *data, const std::vector<unsigned int> &nummodes,
557 const int mode_offset, NekDouble *coeffs,
558 [[maybe_unused]] std::vector<LibUtilities::BasisType> &fromType)
559{
560 int data_order0 = nummodes[mode_offset];
561 int fillorder0 = min(m_base[0]->GetNumModes(), data_order0);
562 int data_order1 = nummodes[mode_offset + 1];
563 int order1 = m_base[1]->GetNumModes();
564 int fillorder1 = min(order1, data_order1);
565 int data_order2 = nummodes[mode_offset + 2];
566 int order2 = m_base[2]->GetNumModes();
567 int fillorder2 = min(order2, data_order2);
568
569 switch (m_base[0]->GetBasisType())
570 {
572 {
573 int i, j;
574 int cnt = 0;
575 int cnt1 = 0;
576
578 "Extraction routine not set up for this basis");
580 "Extraction routine not set up for this basis");
581
582 Vmath::Zero(m_ncoeffs, coeffs, 1);
583 for (j = 0; j < fillorder0; ++j)
584 {
585 for (i = 0; i < fillorder1 - j; ++i)
586 {
587 Vmath::Vcopy(fillorder2 - j - i, &data[cnt], 1,
588 &coeffs[cnt1], 1);
589 cnt += data_order2 - j - i;
590 cnt1 += order2 - j - i;
591 }
592
593 // count out data for j iteration
594 for (i = fillorder1 - j; i < data_order1 - j; ++i)
595 {
596 cnt += data_order2 - j - i;
597 }
598
599 for (i = fillorder1 - j; i < order1 - j; ++i)
600 {
601 cnt1 += order2 - j - i;
602 }
603 }
604 }
605 break;
606 default:
607 ASSERTL0(false, "basis is either not set up or not "
608 "hierarchicial");
609 }
610}
611
612/**
613 * \brief Returns the physical values at the quadrature points of a face
614 */
615void TetExp::v_GetTracePhysMap(const int face, Array<OneD, int> &outarray)
616{
617 int nquad0 = m_base[0]->GetNumPoints();
618 int nquad1 = m_base[1]->GetNumPoints();
619 int nquad2 = m_base[2]->GetNumPoints();
620
621 int nq0 = 0;
622 int nq1 = 0;
623
624 // get forward aligned faces.
625 switch (face)
626 {
627 case 0:
628 {
629 nq0 = nquad0;
630 nq1 = nquad1;
631 if (outarray.size() != nq0 * nq1)
632 {
633 outarray = Array<OneD, int>(nq0 * nq1);
634 }
635
636 for (int i = 0; i < nquad0 * nquad1; ++i)
637 {
638 outarray[i] = i;
639 }
640
641 break;
642 }
643 case 1:
644 {
645 nq0 = nquad0;
646 nq1 = nquad2;
647 if (outarray.size() != nq0 * nq1)
648 {
649 outarray = Array<OneD, int>(nq0 * nq1);
650 }
651
652 // Direction A and B positive
653 for (int k = 0; k < nquad2; k++)
654 {
655 for (int i = 0; i < nquad0; ++i)
656 {
657 outarray[k * nquad0 + i] = (nquad0 * nquad1 * k) + i;
658 }
659 }
660 break;
661 }
662 case 2:
663 {
664 nq0 = nquad1;
665 nq1 = nquad2;
666 if (outarray.size() != nq0 * nq1)
667 {
668 outarray = Array<OneD, int>(nq0 * nq1);
669 }
670
671 // Directions A and B positive
672 for (int j = 0; j < nquad1 * nquad2; ++j)
673 {
674 outarray[j] = nquad0 - 1 + j * nquad0;
675 }
676 break;
677 }
678 case 3:
679 {
680 nq0 = nquad1;
681 nq1 = nquad2;
682 if (outarray.size() != nq0 * nq1)
683 {
684 outarray = Array<OneD, int>(nq0 * nq1);
685 }
686
687 // Directions A and B positive
688 for (int j = 0; j < nquad1 * nquad2; ++j)
689 {
690 outarray[j] = j * nquad0;
691 }
692 }
693 break;
694 default:
695 ASSERTL0(false, "face value (> 3) is out of range");
696 break;
697 }
698}
699
700/**
701 * \brief Compute the normal of a triangular face
702 */
704{
705 int i;
706 const SpatialDomains::GeomFactorsSharedPtr &geomFactors =
707 GetGeom()->GetMetricInfo();
708
710 for (int i = 0; i < ptsKeys.size(); ++i)
711 {
712 // Need at least 2 points for computing normals
713 if (ptsKeys[i].GetNumPoints() == 1)
714 {
715 LibUtilities::PointsKey pKey(2, ptsKeys[i].GetPointsType());
716 ptsKeys[i] = pKey;
717 }
718 }
719
720 SpatialDomains::GeomType type = geomFactors->GetGtype();
722 geomFactors->GetDerivFactors(ptsKeys);
723 const Array<OneD, const NekDouble> &jac = geomFactors->GetJac(ptsKeys);
724
725 LibUtilities::BasisKey tobasis0 = GetTraceBasisKey(face, 0);
726 LibUtilities::BasisKey tobasis1 = GetTraceBasisKey(face, 1);
727
728 // number of face quadrature points
729 int nq_face = tobasis0.GetNumPoints() * tobasis1.GetNumPoints();
730
731 int vCoordDim = GetCoordim();
732
735 for (i = 0; i < vCoordDim; ++i)
736 {
737 normal[i] = Array<OneD, NekDouble>(nq_face);
738 }
739
740 size_t nqb = nq_face;
741 size_t nbnd = face;
744
745 // Regular geometry case
746 if (type == SpatialDomains::eRegular ||
748 {
749 NekDouble fac;
750
751 // Set up normals
752 switch (face)
753 {
754 case 0:
755 {
756 for (i = 0; i < vCoordDim; ++i)
757 {
758 normal[i][0] = -df[3 * i + 2][0];
759 }
760
761 break;
762 }
763 case 1:
764 {
765 for (i = 0; i < vCoordDim; ++i)
766 {
767 normal[i][0] = -df[3 * i + 1][0];
768 }
769
770 break;
771 }
772 case 2:
773 {
774 for (i = 0; i < vCoordDim; ++i)
775 {
776 normal[i][0] =
777 df[3 * i][0] + df[3 * i + 1][0] + df[3 * i + 2][0];
778 }
779
780 break;
781 }
782 case 3:
783 {
784 for (i = 0; i < vCoordDim; ++i)
785 {
786 normal[i][0] = -df[3 * i][0];
787 }
788 break;
789 }
790 default:
791 ASSERTL0(false, "face is out of range (edge < 3)");
792 }
793
794 // normalise
795 fac = 0.0;
796 for (i = 0; i < vCoordDim; ++i)
797 {
798 fac += normal[i][0] * normal[i][0];
799 }
800 fac = 1.0 / sqrt(fac);
801 Vmath::Fill(nqb, fac, length, 1);
802
803 for (i = 0; i < vCoordDim; ++i)
804 {
805 Vmath::Fill(nq_face, fac * normal[i][0], normal[i], 1);
806 }
807 }
808 else
809 {
810 // Set up deformed normals
811 int j, k;
812
813 int nq0 = ptsKeys[0].GetNumPoints();
814 int nq1 = ptsKeys[1].GetNumPoints();
815 int nq2 = ptsKeys[2].GetNumPoints();
816 int nqtot;
817 int nq01 = nq0 * nq1;
818
819 // number of elemental quad points
820 if (face == 0)
821 {
822 nqtot = nq01;
823 }
824 else if (face == 1)
825 {
826 nqtot = nq0 * nq2;
827 }
828 else
829 {
830 nqtot = nq1 * nq2;
831 }
832
835
836 Array<OneD, NekDouble> faceJac(nqtot);
837 Array<OneD, NekDouble> normals(vCoordDim * nqtot, 0.0);
838
839 // Extract Jacobian along face and recover local derivates
840 // (dx/dr) for polynomial interpolation by multiplying m_gmat by
841 // jacobian
842 switch (face)
843 {
844 case 0:
845 {
846 for (j = 0; j < nq01; ++j)
847 {
848 normals[j] = -df[2][j] * jac[j];
849 normals[nqtot + j] = -df[5][j] * jac[j];
850 normals[2 * nqtot + j] = -df[8][j] * jac[j];
851 faceJac[j] = jac[j];
852 }
853
854 points0 = ptsKeys[0];
855 points1 = ptsKeys[1];
856 break;
857 }
858
859 case 1:
860 {
861 for (j = 0; j < nq0; ++j)
862 {
863 for (k = 0; k < nq2; ++k)
864 {
865 int tmp = j + nq01 * k;
866 normals[j + k * nq0] = -df[1][tmp] * jac[tmp];
867 normals[nqtot + j + k * nq0] = -df[4][tmp] * jac[tmp];
868 normals[2 * nqtot + j + k * nq0] =
869 -df[7][tmp] * jac[tmp];
870 faceJac[j + k * nq0] = jac[tmp];
871 }
872 }
873
874 points0 = ptsKeys[0];
875 points1 = ptsKeys[2];
876 break;
877 }
878
879 case 2:
880 {
881 for (j = 0; j < nq1; ++j)
882 {
883 for (k = 0; k < nq2; ++k)
884 {
885 int tmp = nq0 - 1 + nq0 * j + nq01 * k;
886 normals[j + k * nq1] =
887 (df[0][tmp] + df[1][tmp] + df[2][tmp]) * jac[tmp];
888 normals[nqtot + j + k * nq1] =
889 (df[3][tmp] + df[4][tmp] + df[5][tmp]) * jac[tmp];
890 normals[2 * nqtot + j + k * nq1] =
891 (df[6][tmp] + df[7][tmp] + df[8][tmp]) * jac[tmp];
892 faceJac[j + k * nq1] = jac[tmp];
893 }
894 }
895
896 points0 = ptsKeys[1];
897 points1 = ptsKeys[2];
898 break;
899 }
900
901 case 3:
902 {
903 for (j = 0; j < nq1; ++j)
904 {
905 for (k = 0; k < nq2; ++k)
906 {
907 int tmp = j * nq0 + nq01 * k;
908 normals[j + k * nq1] = -df[0][tmp] * jac[tmp];
909 normals[nqtot + j + k * nq1] = -df[3][tmp] * jac[tmp];
910 normals[2 * nqtot + j + k * nq1] =
911 -df[6][tmp] * jac[tmp];
912 faceJac[j + k * nq1] = jac[tmp];
913 }
914 }
915
916 points0 = ptsKeys[1];
917 points1 = ptsKeys[2];
918 break;
919 }
920
921 default:
922 ASSERTL0(false, "face is out of range (face < 3)");
923 }
924
925 Array<OneD, NekDouble> work(nq_face, 0.0);
926 // Interpolate Jacobian and invert
927 LibUtilities::Interp2D(points0, points1, faceJac,
928 tobasis0.GetPointsKey(), tobasis1.GetPointsKey(),
929 work);
930 Vmath::Sdiv(nq_face, 1.0, &work[0], 1, &work[0], 1);
931
932 // Interpolate normal and multiply by inverse Jacobian.
933 for (i = 0; i < vCoordDim; ++i)
934 {
935 LibUtilities::Interp2D(points0, points1, &normals[i * nqtot],
936 tobasis0.GetPointsKey(),
937 tobasis1.GetPointsKey(), &normal[i][0]);
938 Vmath::Vmul(nq_face, work, 1, normal[i], 1, normal[i], 1);
939 }
940
941 // Normalise to obtain unit normals.
942 Vmath::Zero(nq_face, work, 1);
943 for (i = 0; i < GetCoordim(); ++i)
944 {
945 Vmath::Vvtvp(nq_face, normal[i], 1, normal[i], 1, work, 1, work, 1);
946 }
947
948 Vmath::Vsqrt(nq_face, work, 1, work, 1);
949 Vmath::Sdiv(nq_face, 1.0, work, 1, work, 1);
950
951 Vmath::Vcopy(nqb, work, 1, length, 1);
952
953 for (i = 0; i < GetCoordim(); ++i)
954 {
955 Vmath::Vmul(nq_face, normal[i], 1, work, 1, normal[i], 1);
956 }
957 }
958}
959
960//-----------------------------
961// Operator creation functions
962//-----------------------------
964 Array<OneD, NekDouble> &outarray,
965 const StdRegions::StdMatrixKey &mkey)
966{
967 TetExp::v_HelmholtzMatrixOp_MatFree(inarray, outarray, mkey);
968}
969
971 Array<OneD, NekDouble> &outarray,
972 const StdRegions::StdMatrixKey &mkey)
973{
974 TetExp::v_LaplacianMatrixOp_MatFree(inarray, outarray, mkey);
975}
976
977void TetExp::v_LaplacianMatrixOp(const int k1, const int k2,
978 const Array<OneD, const NekDouble> &inarray,
979 Array<OneD, NekDouble> &outarray,
980 const StdRegions::StdMatrixKey &mkey)
981{
982 StdExpansion::LaplacianMatrixOp_MatFree(k1, k2, inarray, outarray, mkey);
983}
984
986 const StdRegions::StdMatrixKey &mkey)
987{
988 int nq = GetTotPoints();
989
990 // Calculate sqrt of the Jacobian
992 Array<OneD, NekDouble> sqrt_jac(nq);
993 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
994 {
995 Vmath::Vsqrt(nq, jac, 1, sqrt_jac, 1);
996 }
997 else
998 {
999 Vmath::Fill(nq, sqrt(jac[0]), sqrt_jac, 1);
1000 }
1001
1002 // Multiply array by sqrt(Jac)
1003 Vmath::Vmul(nq, sqrt_jac, 1, array, 1, array, 1);
1004
1005 // Apply std region filter
1006 StdTetExp::v_SVVLaplacianFilter(array, mkey);
1007
1008 // Divide by sqrt(Jac)
1009 Vmath::Vdiv(nq, array, 1, sqrt_jac, 1, array, 1);
1010}
1011
1012//-----------------------------
1013// Matrix creation functions
1014//-----------------------------
1016{
1017 DNekMatSharedPtr returnval;
1018
1019 switch (mkey.GetMatrixType())
1020 {
1028 returnval = Expansion3D::v_GenMatrix(mkey);
1029 break;
1030 default:
1031 returnval = StdTetExp::v_GenMatrix(mkey);
1032 }
1033
1034 return returnval;
1035}
1036
1038{
1039 LibUtilities::BasisKey bkey0 = m_base[0]->GetBasisKey();
1040 LibUtilities::BasisKey bkey1 = m_base[1]->GetBasisKey();
1041 LibUtilities::BasisKey bkey2 = m_base[2]->GetBasisKey();
1044
1045 return tmp->GetStdMatrix(mkey);
1046}
1047
1049{
1050 return m_matrixManager[mkey];
1051}
1052
1054{
1055 m_matrixManager.DeleteObject(mkey);
1056}
1057
1059{
1060 return m_staticCondMatrixManager[mkey];
1061}
1062
1064{
1065 m_staticCondMatrixManager.DeleteObject(mkey);
1066}
1067
1069 Array<OneD, NekDouble> &outarray,
1070 const StdRegions::StdMatrixKey &mkey)
1071{
1073
1074 if (inarray.data() == outarray.data())
1075 {
1077 Vmath::Vcopy(m_ncoeffs, inarray.data(), 1, tmp.data(), 1);
1078
1079 Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, mat->Scale(),
1080 (mat->GetOwnedMatrix())->GetPtr().data(), m_ncoeffs,
1081 tmp.data(), 1, 0.0, outarray.data(), 1);
1082 }
1083 else
1084 {
1085 Blas::Dgemv('N', m_ncoeffs, m_ncoeffs, mat->Scale(),
1086 (mat->GetOwnedMatrix())->GetPtr().data(), m_ncoeffs,
1087 inarray.data(), 1, 0.0, outarray.data(), 1);
1088 }
1089}
1090
1092 const Array<OneD, const NekDouble> &inarray,
1094{
1095 // This implementation is only valid when there are no
1096 // coefficients associated to the Laplacian operator
1097 if (m_metrics.count(eMetricLaplacian00) == 0)
1098 {
1100 }
1101
1102 int nquad0 = m_base[0]->GetNumPoints();
1103 int nquad1 = m_base[1]->GetNumPoints();
1104 int nquad2 = m_base[2]->GetNumPoints();
1105 int nqtot = nquad0 * nquad1 * nquad2;
1106
1107 ASSERTL1(wsp.size() >= 6 * nqtot, "Insufficient workspace size.");
1108 ASSERTL1(m_ncoeffs <= nqtot, "Workspace not set up for ncoeffs > nqtot");
1109
1110 const Array<OneD, const NekDouble> &base0 = m_base[0]->GetBdata();
1111 const Array<OneD, const NekDouble> &base1 = m_base[1]->GetBdata();
1112 const Array<OneD, const NekDouble> &base2 = m_base[2]->GetBdata();
1113 const Array<OneD, const NekDouble> &dbase0 = m_base[0]->GetDbdata();
1114 const Array<OneD, const NekDouble> &dbase1 = m_base[1]->GetDbdata();
1115 const Array<OneD, const NekDouble> &dbase2 = m_base[2]->GetDbdata();
1116 const Array<OneD, const NekDouble> &metric00 =
1117 m_metrics[eMetricLaplacian00];
1118 const Array<OneD, const NekDouble> &metric01 =
1119 m_metrics[eMetricLaplacian01];
1120 const Array<OneD, const NekDouble> &metric02 =
1121 m_metrics[eMetricLaplacian02];
1122 const Array<OneD, const NekDouble> &metric11 =
1123 m_metrics[eMetricLaplacian11];
1124 const Array<OneD, const NekDouble> &metric12 =
1125 m_metrics[eMetricLaplacian12];
1126 const Array<OneD, const NekDouble> &metric22 =
1127 m_metrics[eMetricLaplacian22];
1128
1129 // Allocate temporary storage
1130 Array<OneD, NekDouble> wsp0(2 * nqtot, wsp);
1131 Array<OneD, NekDouble> wsp1(nqtot, wsp + 1 * nqtot);
1132 Array<OneD, NekDouble> wsp2(nqtot, wsp + 2 * nqtot);
1133 Array<OneD, NekDouble> wsp3(nqtot, wsp + 3 * nqtot);
1134 Array<OneD, NekDouble> wsp4(nqtot, wsp + 4 * nqtot);
1135 Array<OneD, NekDouble> wsp5(nqtot, wsp + 5 * nqtot);
1136
1137 // LAPLACIAN MATRIX OPERATION
1138 // wsp1 = du_dxi1 = D_xi1 * inarray = D_xi1 * u
1139 // wsp2 = du_dxi2 = D_xi2 * inarray = D_xi2 * u
1140 // wsp2 = du_dxi3 = D_xi3 * inarray = D_xi3 * u
1141 StdExpansion3D::PhysTensorDeriv(inarray, wsp0, wsp1, wsp2);
1142
1143 // wsp0 = k = g0 * wsp1 + g1 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1144 // wsp2 = l = g1 * wsp1 + g2 * wsp2 = g0 * du_dxi1 + g1 * du_dxi2
1145 // where g0, g1 and g2 are the metric terms set up in the GeomFactors class
1146 // especially for this purpose
1147 Vmath::Vvtvvtp(nqtot, &metric00[0], 1, &wsp0[0], 1, &metric01[0], 1,
1148 &wsp1[0], 1, &wsp3[0], 1);
1149 Vmath::Vvtvp(nqtot, &metric02[0], 1, &wsp2[0], 1, &wsp3[0], 1, &wsp3[0], 1);
1150 Vmath::Vvtvvtp(nqtot, &metric01[0], 1, &wsp0[0], 1, &metric11[0], 1,
1151 &wsp1[0], 1, &wsp4[0], 1);
1152 Vmath::Vvtvp(nqtot, &metric12[0], 1, &wsp2[0], 1, &wsp4[0], 1, &wsp4[0], 1);
1153 Vmath::Vvtvvtp(nqtot, &metric02[0], 1, &wsp0[0], 1, &metric12[0], 1,
1154 &wsp1[0], 1, &wsp5[0], 1);
1155 Vmath::Vvtvp(nqtot, &metric22[0], 1, &wsp2[0], 1, &wsp5[0], 1, &wsp5[0], 1);
1156
1157 // outarray = m = (D_xi1 * B)^T * k
1158 // wsp1 = n = (D_xi2 * B)^T * l
1159 IProductWRTBase_SumFacKernel(dbase0, base1, base2, wsp3, outarray, wsp0,
1160 false, true, true);
1161 IProductWRTBase_SumFacKernel(base0, dbase1, base2, wsp4, wsp2, wsp0, true,
1162 false, true);
1163 Vmath::Vadd(m_ncoeffs, wsp2.data(), 1, outarray.data(), 1, outarray.data(),
1164 1);
1165 IProductWRTBase_SumFacKernel(base0, base1, dbase2, wsp5, wsp2, wsp0, true,
1166 true, false);
1167 Vmath::Vadd(m_ncoeffs, wsp2.data(), 1, outarray.data(), 1, outarray.data(),
1168 1);
1169}
1170
1171void TetExp::v_ComputeLaplacianMetric()
1172{
1173 if (m_metrics.count(eMetricQuadrature) == 0)
1174 {
1175 ComputeQuadratureMetric();
1176 }
1177
1178 int i, j;
1179 const unsigned int nqtot = GetTotPoints();
1180 const unsigned int dim = 3;
1181 const MetricType m[3][3] = {
1182 {eMetricLaplacian00, eMetricLaplacian01, eMetricLaplacian02},
1183 {eMetricLaplacian01, eMetricLaplacian11, eMetricLaplacian12},
1184 {eMetricLaplacian02, eMetricLaplacian12, eMetricLaplacian22}};
1185
1186 for (unsigned int i = 0; i < dim; ++i)
1187 {
1188 for (unsigned int j = i; j < dim; ++j)
1189 {
1190 m_metrics[m[i][j]] = Array<OneD, NekDouble>(nqtot);
1191 }
1192 }
1193
1194 // Define shorthand synonyms for m_metrics storage
1195 Array<OneD, NekDouble> g0(m_metrics[m[0][0]]);
1196 Array<OneD, NekDouble> g1(m_metrics[m[1][1]]);
1197 Array<OneD, NekDouble> g2(m_metrics[m[2][2]]);
1198 Array<OneD, NekDouble> g3(m_metrics[m[0][1]]);
1199 Array<OneD, NekDouble> g4(m_metrics[m[0][2]]);
1200 Array<OneD, NekDouble> g5(m_metrics[m[1][2]]);
1201
1202 // Allocate temporary storage
1203 Array<OneD, NekDouble> alloc(7 * nqtot, 0.0);
1204 Array<OneD, NekDouble> h0(alloc); // h0
1205 Array<OneD, NekDouble> h1(alloc + 1 * nqtot); // h1
1206 Array<OneD, NekDouble> h2(alloc + 2 * nqtot); // h2
1207 Array<OneD, NekDouble> h3(alloc + 3 * nqtot); // h3
1208 Array<OneD, NekDouble> wsp4(alloc + 4 * nqtot); // wsp4
1209 Array<OneD, NekDouble> wsp5(alloc + 5 * nqtot); // wsp5
1210 Array<OneD, NekDouble> wsp6(alloc + 6 * nqtot); // wsp6
1211 // Reuse some of the storage as workspace
1212 Array<OneD, NekDouble> wsp7(alloc); // wsp7
1213 Array<OneD, NekDouble> wsp8(alloc + 1 * nqtot); // wsp8
1214 Array<OneD, NekDouble> wsp9(alloc + 2 * nqtot); // wsp9
1215
1216 const Array<TwoD, const NekDouble> &df =
1217 m_metricinfo->GetDerivFactors(GetPointsKeys());
1218 const Array<OneD, const NekDouble> &z0 = m_base[0]->GetZ();
1219 const Array<OneD, const NekDouble> &z1 = m_base[1]->GetZ();
1220 const Array<OneD, const NekDouble> &z2 = m_base[2]->GetZ();
1221 const unsigned int nquad0 = m_base[0]->GetNumPoints();
1222 const unsigned int nquad1 = m_base[1]->GetNumPoints();
1223 const unsigned int nquad2 = m_base[2]->GetNumPoints();
1224
1225 for (j = 0; j < nquad2; ++j)
1226 {
1227 for (i = 0; i < nquad1; ++i)
1228 {
1229 Vmath::Fill(nquad0, 4.0 / (1.0 - z1[i]) / (1.0 - z2[j]),
1230 &h0[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1231 Vmath::Fill(nquad0, 2.0 / (1.0 - z1[i]) / (1.0 - z2[j]),
1232 &h1[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1233 Vmath::Fill(nquad0, 2.0 / (1.0 - z2[j]),
1234 &h2[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1235 Vmath::Fill(nquad0, (1.0 + z1[i]) / (1.0 - z2[j]),
1236 &h3[0] + i * nquad0 + j * nquad0 * nquad1, 1);
1237 }
1238 }
1239 for (i = 0; i < nquad0; i++)
1240 {
1241 Blas::Dscal(nquad1 * nquad2, 1 + z0[i], &h1[0] + i, nquad0);
1242 }
1243
1244 // Step 3. Construct combined metric terms for physical space to
1245 // collapsed coordinate system.
1246 // Order of construction optimised to minimise temporary storage
1247 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1248 {
1249 // wsp4
1250 Vmath::Vadd(nqtot, &df[1][0], 1, &df[2][0], 1, &wsp4[0], 1);
1251 Vmath::Vvtvvtp(nqtot, &df[0][0], 1, &h0[0], 1, &wsp4[0], 1, &h1[0], 1,
1252 &wsp4[0], 1);
1253 // wsp5
1254 Vmath::Vadd(nqtot, &df[4][0], 1, &df[5][0], 1, &wsp5[0], 1);
1255 Vmath::Vvtvvtp(nqtot, &df[3][0], 1, &h0[0], 1, &wsp5[0], 1, &h1[0], 1,
1256 &wsp5[0], 1);
1257 // wsp6
1258 Vmath::Vadd(nqtot, &df[7][0], 1, &df[8][0], 1, &wsp6[0], 1);
1259 Vmath::Vvtvvtp(nqtot, &df[6][0], 1, &h0[0], 1, &wsp6[0], 1, &h1[0], 1,
1260 &wsp6[0], 1);
1261
1262 // g0
1263 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1264 1, &g0[0], 1);
1265 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1266
1267 // g4
1268 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp4[0], 1, &df[5][0], 1, &wsp5[0],
1269 1, &g4[0], 1);
1270 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1271
1272 // overwrite h0, h1, h2
1273 // wsp7 (h2f1 + h3f2)
1274 Vmath::Vvtvvtp(nqtot, &df[1][0], 1, &h2[0], 1, &df[2][0], 1, &h3[0], 1,
1275 &wsp7[0], 1);
1276 // wsp8 (h2f4 + h3f5)
1277 Vmath::Vvtvvtp(nqtot, &df[4][0], 1, &h2[0], 1, &df[5][0], 1, &h3[0], 1,
1278 &wsp8[0], 1);
1279 // wsp9 (h2f7 + h3f8)
1280 Vmath::Vvtvvtp(nqtot, &df[7][0], 1, &h2[0], 1, &df[8][0], 1, &h3[0], 1,
1281 &wsp9[0], 1);
1282
1283 // g3
1284 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp7[0], 1, &wsp5[0], 1, &wsp8[0],
1285 1, &g3[0], 1);
1286 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp9[0], 1, &g3[0], 1, &g3[0], 1);
1287
1288 // overwrite wsp4, wsp5, wsp6
1289 // g1
1290 Vmath::Vvtvvtp(nqtot, &wsp7[0], 1, &wsp7[0], 1, &wsp8[0], 1, &wsp8[0],
1291 1, &g1[0], 1);
1292 Vmath::Vvtvp(nqtot, &wsp9[0], 1, &wsp9[0], 1, &g1[0], 1, &g1[0], 1);
1293
1294 // g5
1295 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &wsp7[0], 1, &df[5][0], 1, &wsp8[0],
1296 1, &g5[0], 1);
1297 Vmath::Vvtvp(nqtot, &df[8][0], 1, &wsp9[0], 1, &g5[0], 1, &g5[0], 1);
1298
1299 // g2
1300 Vmath::Vvtvvtp(nqtot, &df[2][0], 1, &df[2][0], 1, &df[5][0], 1,
1301 &df[5][0], 1, &g2[0], 1);
1302 Vmath::Vvtvp(nqtot, &df[8][0], 1, &df[8][0], 1, &g2[0], 1, &g2[0], 1);
1303 }
1304 else
1305 {
1306 // wsp4
1307 Vmath::Svtsvtp(nqtot, df[0][0], &h0[0], 1, df[1][0] + df[2][0], &h1[0],
1308 1, &wsp4[0], 1);
1309 // wsp5
1310 Vmath::Svtsvtp(nqtot, df[3][0], &h0[0], 1, df[4][0] + df[5][0], &h1[0],
1311 1, &wsp5[0], 1);
1312 // wsp6
1313 Vmath::Svtsvtp(nqtot, df[6][0], &h0[0], 1, df[7][0] + df[8][0], &h1[0],
1314 1, &wsp6[0], 1);
1315
1316 // g0
1317 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp4[0], 1, &wsp5[0], 1, &wsp5[0],
1318 1, &g0[0], 1);
1319 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp6[0], 1, &g0[0], 1, &g0[0], 1);
1320
1321 // g4
1322 Vmath::Svtsvtp(nqtot, df[2][0], &wsp4[0], 1, df[5][0], &wsp5[0], 1,
1323 &g4[0], 1);
1324 Vmath::Svtvp(nqtot, df[8][0], &wsp6[0], 1, &g4[0], 1, &g4[0], 1);
1325
1326 // overwrite h0, h1, h2
1327 // wsp7 (h2f1 + h3f2)
1328 Vmath::Svtsvtp(nqtot, df[1][0], &h2[0], 1, df[2][0], &h3[0], 1,
1329 &wsp7[0], 1);
1330 // wsp8 (h2f4 + h3f5)
1331 Vmath::Svtsvtp(nqtot, df[4][0], &h2[0], 1, df[5][0], &h3[0], 1,
1332 &wsp8[0], 1);
1333 // wsp9 (h2f7 + h3f8)
1334 Vmath::Svtsvtp(nqtot, df[7][0], &h2[0], 1, df[8][0], &h3[0], 1,
1335 &wsp9[0], 1);
1336
1337 // g3
1338 Vmath::Vvtvvtp(nqtot, &wsp4[0], 1, &wsp7[0], 1, &wsp5[0], 1, &wsp8[0],
1339 1, &g3[0], 1);
1340 Vmath::Vvtvp(nqtot, &wsp6[0], 1, &wsp9[0], 1, &g3[0], 1, &g3[0], 1);
1341
1342 // overwrite wsp4, wsp5, wsp6
1343 // g1
1344 Vmath::Vvtvvtp(nqtot, &wsp7[0], 1, &wsp7[0], 1, &wsp8[0], 1, &wsp8[0],
1345 1, &g1[0], 1);
1346 Vmath::Vvtvp(nqtot, &wsp9[0], 1, &wsp9[0], 1, &g1[0], 1, &g1[0], 1);
1347
1348 // g5
1349 Vmath::Svtsvtp(nqtot, df[2][0], &wsp7[0], 1, df[5][0], &wsp8[0], 1,
1350 &g5[0], 1);
1351 Vmath::Svtvp(nqtot, df[8][0], &wsp9[0], 1, &g5[0], 1, &g5[0], 1);
1352
1353 // g2
1354 Vmath::Fill(nqtot,
1355 df[2][0] * df[2][0] + df[5][0] * df[5][0] +
1356 df[8][0] * df[8][0],
1357 &g2[0], 1);
1358 }
1359
1360 for (unsigned int i = 0; i < dim; ++i)
1361 {
1362 for (unsigned int j = i; j < dim; ++j)
1363 {
1364 MultiplyByQuadratureMetric(m_metrics[m[i][j]], m_metrics[m[i][j]]);
1365 }
1366 }
1367}
1368
1369/** @brief: This method gets all of the factors which are
1370 required as part of the Gradient Jump Penalty
1371 stabilisation and involves the product of the normal and
1372 geometric factors along the element trace.
1373*/
1374void TetExp::v_NormalTraceDerivFactors(
1375 Array<OneD, Array<OneD, NekDouble>> &d0factors,
1376 Array<OneD, Array<OneD, NekDouble>> &d1factors,
1377 Array<OneD, Array<OneD, NekDouble>> &d2factors)
1378{
1379 int nquad0 = GetNumPoints(0);
1380 int nquad1 = GetNumPoints(1);
1381 int nquad2 = GetNumPoints(2);
1382
1383 const Array<TwoD, const NekDouble> &df =
1384 m_metricinfo->GetDerivFactors(GetPointsKeys());
1385
1386 if (d0factors.size() != 4)
1387 {
1388 d0factors = Array<OneD, Array<OneD, NekDouble>>(4);
1389 d1factors = Array<OneD, Array<OneD, NekDouble>>(4);
1390 d2factors = Array<OneD, Array<OneD, NekDouble>>(4);
1391 }
1392
1393 if (d0factors[0].size() != nquad0 * nquad1)
1394 {
1395 d0factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1396 d1factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1397 d2factors[0] = Array<OneD, NekDouble>(nquad0 * nquad1);
1398 }
1399
1400 if (d0factors[1].size() != nquad0 * nquad2)
1401 {
1402 d0factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1403 d1factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1404 d2factors[1] = Array<OneD, NekDouble>(nquad0 * nquad2);
1405 }
1406
1407 if (d0factors[2].size() != nquad1 * nquad2)
1408 {
1409 d0factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1410 d0factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1411 d1factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1412 d1factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1413 d2factors[2] = Array<OneD, NekDouble>(nquad1 * nquad2);
1414 d2factors[3] = Array<OneD, NekDouble>(nquad1 * nquad2);
1415 }
1416
1417 // Outwards normals
1418 const Array<OneD, const Array<OneD, NekDouble>> &normal_0 =
1419 GetTraceNormal(0);
1420 const Array<OneD, const Array<OneD, NekDouble>> &normal_1 =
1421 GetTraceNormal(1);
1422 const Array<OneD, const Array<OneD, NekDouble>> &normal_2 =
1423 GetTraceNormal(2);
1424 const Array<OneD, const Array<OneD, NekDouble>> &normal_3 =
1425 GetTraceNormal(3);
1426
1427 int ncoords = normal_0.size();
1428
1429 // first gather together standard cartesian inner products
1430 if (m_metricinfo->GetGtype() == SpatialDomains::eDeformed)
1431 {
1432 // face 0
1433 for (int i = 0; i < nquad0 * nquad1; ++i)
1434 {
1435 d0factors[0][i] = df[0][i] * normal_0[0][i];
1436 d1factors[0][i] = df[1][i] * normal_0[0][i];
1437 d2factors[0][i] = df[2][i] * normal_0[0][i];
1438 }
1439
1440 for (int n = 1; n < ncoords; ++n)
1441 {
1442 for (int i = 0; i < nquad0 * nquad1; ++i)
1443 {
1444 d0factors[0][i] += df[3 * n][i] * normal_0[n][i];
1445 d1factors[0][i] += df[3 * n + 1][i] * normal_0[n][i];
1446 d2factors[0][i] += df[3 * n + 2][i] * normal_0[n][i];
1447 }
1448 }
1449
1450 // face 1
1451 for (int j = 0; j < nquad2; ++j)
1452 {
1453 for (int i = 0; i < nquad0; ++i)
1454 {
1455 d0factors[1][j * nquad0 + i] = df[0][j * nquad0 * nquad1 + i] *
1456 normal_1[0][j * nquad0 + i];
1457 d1factors[1][j * nquad0 + i] = df[1][j * nquad0 * nquad1 + i] *
1458 normal_1[0][j * nquad0 + i];
1459 d2factors[1][j * nquad0 + i] = df[2][j * nquad0 * nquad1 + i] *
1460 normal_1[0][j * nquad0 + i];
1461 }
1462 }
1463
1464 for (int n = 1; n < ncoords; ++n)
1465 {
1466 for (int j = 0; j < nquad2; ++j)
1467 {
1468 for (int i = 0; i < nquad0; ++i)
1469 {
1470 d0factors[1][j * nquad0 + i] +=
1471 df[3 * n][j * nquad0 * nquad1 + i] *
1472 normal_1[0][j * nquad0 + i];
1473 d1factors[1][j * nquad0 + i] +=
1474 df[3 * n + 1][j * nquad0 * nquad1 + i] *
1475 normal_1[0][j * nquad0 + i];
1476 d2factors[1][j * nquad0 + i] +=
1477 df[3 * n + 2][j * nquad0 * nquad1 + i] *
1478 normal_1[0][j * nquad0 + i];
1479 }
1480 }
1481 }
1482
1483 // faces 2 and 3
1484 for (int j = 0; j < nquad2; ++j)
1485 {
1486 for (int i = 0; i < nquad1; ++i)
1487 {
1488 d0factors[2][j * nquad1 + i] =
1489 df[0][(j * nquad1 + i + 1) * nquad0 - 1] *
1490 normal_2[0][j * nquad1 + i];
1491 d1factors[2][j * nquad1 + i] =
1492 df[1][(j * nquad1 + i + 1) * nquad0 - 1] *
1493 normal_2[0][j * nquad1 + i];
1494 d2factors[2][j * nquad1 + i] =
1495 df[2][(j * nquad1 + i + 1) * nquad0 - 1] *
1496 normal_2[0][j * nquad1 + i];
1497
1498 d0factors[3][j * nquad1 + i] =
1499 df[0][(j * nquad1 + i) * nquad0] *
1500 normal_3[0][j * nquad1 + i];
1501 d1factors[3][j * nquad1 + i] =
1502 df[1][(j * nquad1 + i) * nquad0] *
1503 normal_3[0][j * nquad1 + i];
1504 d2factors[3][j * nquad1 + i] =
1505 df[2][(j * nquad1 + i) * nquad0] *
1506 normal_3[0][j * nquad1 + i];
1507 }
1508 }
1509
1510 for (int n = 1; n < ncoords; ++n)
1511 {
1512 for (int j = 0; j < nquad2; ++j)
1513 {
1514 for (int i = 0; i < nquad1; ++i)
1515 {
1516 d0factors[2][j * nquad1 + i] +=
1517 df[3 * n][(j * nquad1 + i + 1) * nquad0 - 1] *
1518 normal_2[n][j * nquad1 + i];
1519 d1factors[2][j * nquad1 + i] +=
1520 df[3 * n + 1][(j * nquad1 + i + 1) * nquad0 - 1] *
1521 normal_2[n][j * nquad1 + i];
1522 d2factors[2][j * nquad1 + i] +=
1523 df[3 * n + 2][(j * nquad1 + i + 1) * nquad0 - 1] *
1524 normal_2[n][j * nquad1 + i];
1525
1526 d0factors[3][j * nquad1 + i] +=
1527 df[3 * n][(j * nquad1 + i) * nquad0] *
1528 normal_3[n][j * nquad1 + i];
1529 d1factors[3][j * nquad1 + i] +=
1530 df[3 * n + 1][(j * nquad1 + i) * nquad0] *
1531 normal_3[n][j * nquad1 + i];
1532 d2factors[3][j * nquad1 + i] +=
1533 df[3 * n + 2][(j * nquad1 + i) * nquad0] *
1534 normal_3[n][j * nquad1 + i];
1535 }
1536 }
1537 }
1538 }
1539 else
1540 {
1541 // Face 0
1542 for (int i = 0; i < nquad0 * nquad1; ++i)
1543 {
1544 d0factors[0][i] = df[0][0] * normal_0[0][i];
1545 d1factors[0][i] = df[1][0] * normal_0[0][i];
1546 d2factors[0][i] = df[2][0] * normal_0[0][i];
1547 }
1548
1549 for (int n = 1; n < ncoords; ++n)
1550 {
1551 for (int i = 0; i < nquad0 * nquad1; ++i)
1552 {
1553 d0factors[0][i] += df[3 * n][0] * normal_0[n][i];
1554 d1factors[0][i] += df[3 * n + 1][0] * normal_0[n][i];
1555 d2factors[0][i] += df[3 * n + 2][0] * normal_0[n][i];
1556 }
1557 }
1558
1559 // face 1
1560 for (int i = 0; i < nquad0 * nquad2; ++i)
1561 {
1562 d0factors[1][i] = df[0][0] * normal_1[0][i];
1563 d1factors[1][i] = df[1][0] * normal_1[0][i];
1564 d2factors[1][i] = df[2][0] * normal_1[0][i];
1565 }
1566
1567 for (int n = 1; n < ncoords; ++n)
1568 {
1569 for (int i = 0; i < nquad0 * nquad2; ++i)
1570 {
1571 d0factors[1][i] += df[3 * n][0] * normal_1[n][i];
1572 d1factors[1][i] += df[3 * n + 1][0] * normal_1[n][i];
1573 d2factors[1][i] += df[3 * n + 2][0] * normal_1[n][i];
1574 }
1575 }
1576
1577 // faces 2 and 3
1578 for (int i = 0; i < nquad1 * nquad2; ++i)
1579 {
1580 d0factors[2][i] = df[0][0] * normal_2[0][i];
1581 d0factors[3][i] = df[0][0] * normal_3[0][i];
1582
1583 d1factors[2][i] = df[1][0] * normal_2[0][i];
1584 d1factors[3][i] = df[1][0] * normal_3[0][i];
1585
1586 d2factors[2][i] = df[2][0] * normal_2[0][i];
1587 d2factors[3][i] = df[2][0] * normal_3[0][i];
1588 }
1589
1590 for (int n = 1; n < ncoords; ++n)
1591 {
1592 for (int i = 0; i < nquad1 * nquad2; ++i)
1593 {
1594 d0factors[2][i] += df[3 * n][0] * normal_2[n][i];
1595 d0factors[3][i] += df[3 * n][0] * normal_3[n][i];
1596
1597 d1factors[2][i] += df[3 * n + 1][0] * normal_2[n][i];
1598 d1factors[3][i] += df[3 * n + 1][0] * normal_3[n][i];
1599
1600 d2factors[2][i] += df[3 * n + 2][0] * normal_2[n][i];
1601 d2factors[3][i] += df[3 * n + 2][0] * normal_3[n][i];
1602 }
1603 }
1604 }
1605}
1606} // namespace Nektar::LocalRegions
#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
Describes the specification for a Basis.
Definition: Basis.h:45
int GetNumPoints() const
Return points order at which basis is defined.
Definition: Basis.h:120
PointsKey GetPointsKey() const
Return distribution of points.
Definition: Basis.h:137
Defines a specification for a set of points.
Definition: Points.h:50
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
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
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
Definition: Expansion.cpp:534
SpatialDomains::GeomFactorsSharedPtr m_metricinfo
Definition: Expansion.h:274
DNekScalMatSharedPtr GetLocMatrix(const LocalRegions::MatrixKey &mkey)
Definition: Expansion.cpp:84
void v_ComputeTraceNormal(const int face) override
Compute the normal of a triangular face.
Definition: TetExp.cpp:703
NekDouble v_PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &physvals) override
Definition: TetExp.cpp:471
NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray) override
Integrate the physical point list inarray over region.
Definition: TetExp.cpp:102
DNekScalMatSharedPtr v_GetLocMatrix(const MatrixKey &mkey) override
Definition: TetExp.cpp:1048
NekDouble v_StdPhysEvaluate(const Array< OneD, const NekDouble > &Lcoord, const Array< OneD, const NekDouble > &physvals) override
Definition: TetExp.cpp:459
void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2) override
Differentiate inarray in the three coordinate directions.
Definition: TetExp.cpp:141
void v_DropLocMatrix(const MatrixKey &mkey) override
Definition: TetExp.cpp:1053
LibUtilities::NekManager< MatrixKey, DNekScalBlkMat, MatrixKey::opLess > m_staticCondMatrixManager
Definition: TetExp.h:202
LibUtilities::NekManager< MatrixKey, DNekScalMat, MatrixKey::opLess > m_matrixManager
Definition: TetExp.h:200
DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: TetExp.cpp:1015
void v_LaplacianMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: TetExp.cpp:970
void v_IProductWRTBase_SumFac(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool multiplybyweights=true) override
Definition: TetExp.cpp:282
NekDouble v_PhysEvalFirstDeriv(const Array< OneD, NekDouble > &coord, const Array< OneD, const NekDouble > &inarray, std::array< NekDouble, 3 > &firstOrderDerivs) override
Definition: TetExp.cpp:485
DNekScalBlkMatSharedPtr v_GetLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: TetExp.cpp:1058
void v_GetCoord(const Array< OneD, const NekDouble > &Lcoords, Array< OneD, NekDouble > &coords) override
Get the coordinates "coords" at the local coordinates "Lcoords".
Definition: TetExp.cpp:499
void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculate the inner product of inarray with respect to the basis B=m_base0*m_base1*m_base2 and put in...
Definition: TetExp.cpp:276
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: TetExp.cpp:223
StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const override
Definition: TetExp.cpp:535
LibUtilities::ShapeType v_DetShapeType() const override
Return Shape of region, using ShapeType enum list.
Definition: TetExp.cpp:530
DNekMatSharedPtr v_CreateStdMatrix(const StdRegions::StdMatrixKey &mkey) override
Definition: TetExp.cpp:1037
void v_IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray) override
Calculates the inner product .
Definition: TetExp.cpp:341
StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const override
Definition: TetExp.cpp:542
TetExp(const LibUtilities::BasisKey &Ba, const LibUtilities::BasisKey &Bb, const LibUtilities::BasisKey &Bc, const SpatialDomains::TetGeomSharedPtr &geom)
Constructor using BasisKey class for quadrature points and order definition.
Definition: TetExp.cpp:57
void v_GetCoords(Array< OneD, NekDouble > &coords_1, Array< OneD, NekDouble > &coords_2, Array< OneD, NekDouble > &coords_3) override
Definition: TetExp.cpp:516
void GeneralMatrixOp_MatOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey)
Definition: TetExp.cpp:1068
void v_DropLocStaticCondMatrix(const MatrixKey &mkey) override
Definition: TetExp.cpp:1063
void v_SVVLaplacianFilter(Array< OneD, NekDouble > &array, const StdRegions::StdMatrixKey &mkey) override
Definition: TetExp.cpp:985
void v_AlignVectorToCollapsedDir(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray) override
Definition: TetExp.cpp:386
void v_GetTracePhysMap(const int face, Array< OneD, int > &outarray) override
Returns the physical values at the quadrature points of a face.
Definition: TetExp.cpp:615
void v_HelmholtzMatrixOp(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
Definition: TetExp.cpp:963
void v_LaplacianMatrixOp_MatFree_Kernel(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp) override
Definition: TetExp.cpp:1091
void v_ExtractDataToCoeffs(const NekDouble *data, const std::vector< unsigned int > &nummodes, const int mode_offset, NekDouble *coeffs, std::vector< LibUtilities::BasisType > &fromType) override
Definition: TetExp.cpp:555
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
void IProductWRTBase_SumFacKernel(const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, Array< OneD, NekDouble > &wsp, bool doCheckCollDir0, bool doCheckCollDir1, bool doCheckCollDir2)
void v_HelmholtzMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
void v_LaplacianMatrixOp_MatFree(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::StdMatrixKey &mkey) override
int GetNcoeffs(void) const
This function returns the total number of coefficients used in the expansion.
Definition: StdExpansion.h:124
int GetTotPoints() const
This function returns the total number of quadrature points used in the element.
Definition: StdExpansion.h:134
LibUtilities::BasisType GetBasisType(const int dir) const
This function returns the type of basis used in the dir direction.
Definition: StdExpansion.h:156
const LibUtilities::PointsKeyVector GetPointsKeys() const
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
this function calculates the inner product of a given function f with the different modes of the expa...
Definition: StdExpansion.h:537
LibUtilities::PointsType GetPointsType(const int dir) const
This function returns the type of quadrature points used in the dir direction.
Definition: StdExpansion.h:205
int GetNumPoints(const int dir) const
This function returns the number of quadrature points in the dir direction.
Definition: StdExpansion.h:218
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Definition: StdExpansion.h:732
const LibUtilities::BasisKey GetTraceBasisKey(const int i, int k=-1, bool UseGLL=false) const
This function returns the basis key belonging to the i-th trace.
Definition: StdExpansion.h:301
Array< OneD, LibUtilities::BasisSharedPtr > m_base
MatrixType GetMatrixType() const
Definition: StdMatrixKey.h:83
LibUtilities::ShapeType DetShapeType() const
Definition: StdTetExp.h:56
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 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 Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
Definition: Interp.cpp:101
std::vector< PointsKey > PointsKeyVector
Definition: Points.h:231
@ eModified_B
Principle Modified Functions .
Definition: BasisType.h:49
@ eModified_C
Principle Modified Functions .
Definition: BasisType.h:50
@ 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.
@ eMovingRegular
Currently unused.
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< TetGeom > TetGeomSharedPtr
Definition: TetGeom.h:82
std::shared_ptr< StdTetExp > StdTetExpSharedPtr
Definition: StdTetExp.h:227
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
Definition: NekTypeDefs.hpp:79
std::shared_ptr< DNekMat > DNekMatSharedPtr
Definition: NekTypeDefs.hpp:75
double NekDouble
void Vsqrt(int n, const T *x, const int incx, T *y, const int incy)
sqrt y = sqrt(x)
Definition: Vmath.hpp:340
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 Sdiv(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha/x.
Definition: Vmath.hpp:154
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
STL namespace.
scalarT< T > sqrt(scalarT< T > in)
Definition: scalar.hpp:285