Nektar++
IProductWRTDerivBase.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: IProductWRTDerivBase.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: IProductWRTDerivBase operator implementations
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <Collections/Collection.h>
36#include <Collections/IProduct.h>
37#include <Collections/MatrixFreeBase.h>
38#include <Collections/Operator.h>
39#include <MatrixFreeOps/Operator.hpp>
40
41using namespace std;
42
43namespace Nektar::Collections
44{
45
46using LibUtilities::eHexahedron;
47using LibUtilities::ePrism;
48using LibUtilities::ePyramid;
49using LibUtilities::eQuadrilateral;
50using LibUtilities::eSegment;
51using LibUtilities::eTetrahedron;
52using LibUtilities::eTriangle;
53
54/**
55 * @brief Inner product deriv base help class to calculate the size of the
56 * collection that is given as an input and as an output to the
57 * IProductWRTDerivBase Operator. The size evaluation takes into account the
58 * conversion from the physical space to the coefficient space.
59 */
60class IProductWRTDerivBase_Helper : virtual public Operator
61{
62protected:
63 IProductWRTDerivBase_Helper()
64 {
65 // expect input to be number of elements by the number of quad points
66 m_inputSize = m_numElmt * m_stdExp->GetTotPoints();
67 // expect input to be number of elements by the number of coefficients
68 m_outputSize = m_numElmt * m_stdExp->GetNcoeffs();
69 }
70};
71
72/**
73 * @brief Inner product WRT deriv base operator using standard matrix approach
74 */
75class IProductWRTDerivBase_StdMat final : virtual public Operator,
76 IProductWRTDerivBase_Helper
77{
78public:
79 OPERATOR_CREATE(IProductWRTDerivBase_StdMat)
80
81 ~IProductWRTDerivBase_StdMat() final = default;
82
83 void operator()(const Array<OneD, const NekDouble> &entry0,
84 Array<OneD, NekDouble> &entry1,
85 Array<OneD, NekDouble> &entry2,
86 Array<OneD, NekDouble> &entry3,
87 Array<OneD, NekDouble> &wsp) final
88 {
89 int nPhys = m_stdExp->GetTotPoints();
90 int ntot = m_numElmt * nPhys;
91 int nmodes = m_stdExp->GetNcoeffs();
92 Array<OneD, Array<OneD, const NekDouble>> in(3);
93 Array<OneD, NekDouble> output;
94 Array<OneD, Array<OneD, NekDouble>> tmp(3);
95
96 in[0] = entry0;
97 in[1] = entry1;
98 in[2] = entry2;
99
100 output = (m_coordim == 3) ? entry3 : (m_coordim == 2) ? entry2 : entry1;
101
102 for (int i = 0; i < m_dim; ++i)
103 {
104 tmp[i] = wsp + i * ntot;
105 }
106
107 // calculate Iproduct WRT Std Deriv
108
109 // First component
110 if (m_isDeformed)
111 {
112 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
113 for (int i = 0; i < m_dim; ++i)
114 {
115 Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
116 for (int j = 1; j < m_coordim; ++j)
117 {
118 Vmath::Vvtvp(ntot, m_derivFac[i + j * m_dim], 1, in[j], 1,
119 tmp[i], 1, tmp[i], 1);
120 }
121 }
122
123 Vmath::Vmul(ntot, m_jac, 1, tmp[0], 1, tmp[0], 1);
124 }
125 else
126 {
127 Array<OneD, NekDouble> t;
128 for (int e = 0; e < m_numElmt; ++e)
129 {
130 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
131 for (int i = 0; i < m_dim; ++i)
132 {
133 Vmath::Smul(m_nqe, m_derivFac[i][e], in[0] + e * m_nqe, 1,
134 t = tmp[i] + e * m_nqe, 1);
135 for (int j = 1; j < m_coordim; ++j)
136 {
137 Vmath::Svtvp(m_nqe, m_derivFac[i + j * m_dim][e],
138 in[j] + e * m_nqe, 1, tmp[i] + e * m_nqe,
139 1, t = tmp[i] + e * m_nqe, 1);
140 }
141 }
142
143 Vmath::Smul(m_nqe, m_jac[e], tmp[0] + e * m_nqe, 1,
144 t = tmp[0] + e * m_nqe, 1);
145 }
146 }
147
148 Blas::Dgemm('N', 'N', m_iProdWRTStdDBase[0]->GetRows(), m_numElmt,
149 m_iProdWRTStdDBase[0]->GetColumns(), 1.0,
150 m_iProdWRTStdDBase[0]->GetRawPtr(),
151 m_iProdWRTStdDBase[0]->GetRows(), tmp[0].get(), nPhys, 0.0,
152 output.get(), nmodes);
153
154 // Other components
155 for (int i = 1; i < m_dim; ++i)
156 {
157 if (m_isDeformed)
158 {
159 Vmath::Vmul(ntot, m_jac, 1, tmp[i], 1, tmp[i], 1);
160 }
161 else
162 {
163 Array<OneD, NekDouble> t;
164 for (int e = 0; e < m_numElmt; ++e)
165 {
166 Vmath::Smul(m_nqe, m_jac[e], tmp[i] + e * m_nqe, 1,
167 t = tmp[i] + e * m_nqe, 1);
168 }
169 }
170 Blas::Dgemm('N', 'N', m_iProdWRTStdDBase[i]->GetRows(), m_numElmt,
171 m_iProdWRTStdDBase[i]->GetColumns(), 1.0,
172 m_iProdWRTStdDBase[i]->GetRawPtr(),
173 m_iProdWRTStdDBase[i]->GetRows(), tmp[i].get(), nPhys,
174 1.0, output.get(), nmodes);
175 }
176 }
177
178 void operator()([[maybe_unused]] int dir,
179 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
180 [[maybe_unused]] Array<OneD, NekDouble> &output,
181 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
182 {
183 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
184 }
185
186protected:
187 Array<OneD, DNekMatSharedPtr> m_iProdWRTStdDBase;
188 Array<TwoD, const NekDouble> m_derivFac;
189 Array<OneD, const NekDouble> m_jac;
190 int m_dim;
191 int m_coordim;
192
193private:
194 IProductWRTDerivBase_StdMat(
195 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
196 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
197 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper()
198 {
199 m_dim = pCollExp[0]->GetShapeDimension();
200 m_coordim = pCollExp[0]->GetCoordim();
201
202 m_nqe = m_stdExp->GetTotPoints();
203 int nmodes = m_stdExp->GetNcoeffs();
204
205 // set up a IProductWRTDerivBase StdMat.
206 m_iProdWRTStdDBase = Array<OneD, DNekMatSharedPtr>(m_dim);
207 for (int i = 0; i < m_dim; ++i)
208 {
209 Array<OneD, NekDouble> tmp(m_nqe), tmp1(nmodes);
210 m_iProdWRTStdDBase[i] =
211 MemoryManager<DNekMat>::AllocateSharedPtr(nmodes, m_nqe);
212 for (int j = 0; j < m_nqe; ++j)
213 {
214 Vmath::Zero(m_nqe, tmp, 1);
215 tmp[j] = 1.0;
216 m_stdExp->IProductWRTDerivBase(i, tmp, tmp1);
217 Vmath::Vcopy(nmodes, &tmp1[0], 1,
218 &(m_iProdWRTStdDBase[i]->GetPtr())[0] + j * nmodes,
219 1);
220 }
221 }
222 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
223 m_jac = pGeomData->GetJac(pCollExp);
224 m_wspSize = m_dim * m_nqe * m_numElmt;
225 }
226};
227
228/// Factory initialisation for the IProductWRTDerivBase_StdMat operators
229OperatorKey IProductWRTDerivBase_StdMat::m_typeArr[] = {
230 GetOperatorFactory().RegisterCreatorFunction(
231 OperatorKey(eSegment, eIProductWRTDerivBase, eStdMat, false),
232 IProductWRTDerivBase_StdMat::create, "IProductWRTDerivBase_StdMat_Seg"),
233 GetOperatorFactory().RegisterCreatorFunction(
234 OperatorKey(eTriangle, eIProductWRTDerivBase, eStdMat, false),
235 IProductWRTDerivBase_StdMat::create, "IProductWRTDerivBase_StdMat_Tri"),
236 GetOperatorFactory().RegisterCreatorFunction(
237 OperatorKey(eTriangle, eIProductWRTDerivBase, eStdMat, true),
238 IProductWRTDerivBase_StdMat::create,
239 "IProductWRTDerivBase_StdMat_NodalTri"),
240 GetOperatorFactory().RegisterCreatorFunction(
241 OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eStdMat, false),
242 IProductWRTDerivBase_StdMat::create,
243 "IProductWRTDerivBase_StdMat_Quad"),
244 GetOperatorFactory().RegisterCreatorFunction(
245 OperatorKey(eTetrahedron, eIProductWRTDerivBase, eStdMat, false),
246 IProductWRTDerivBase_StdMat::create, "IProductWRTDerivBase_StdMat_Tet"),
247 GetOperatorFactory().RegisterCreatorFunction(
248 OperatorKey(eTetrahedron, eIProductWRTDerivBase, eStdMat, true),
249 IProductWRTDerivBase_StdMat::create,
250 "IProductWRTDerivBase_StdMat_NodalTet"),
251 GetOperatorFactory().RegisterCreatorFunction(
252 OperatorKey(ePyramid, eIProductWRTDerivBase, eStdMat, false),
253 IProductWRTDerivBase_StdMat::create, "IProductWRTDerivBase_StdMat_Pyr"),
254 GetOperatorFactory().RegisterCreatorFunction(
255 OperatorKey(ePrism, eIProductWRTDerivBase, eStdMat, false),
256 IProductWRTDerivBase_StdMat::create,
257 "IProductWRTDerivBase_StdMat_Prism"),
258 GetOperatorFactory().RegisterCreatorFunction(
259 OperatorKey(ePrism, eIProductWRTDerivBase, eStdMat, true),
260 IProductWRTDerivBase_StdMat::create,
261 "IProductWRTDerivBase_StdMat_NodalPrism"),
262 GetOperatorFactory().RegisterCreatorFunction(
263 OperatorKey(eHexahedron, eIProductWRTDerivBase, eStdMat, false),
264 IProductWRTDerivBase_StdMat::create, "IProductWRTDerivBase_StdMat_Hex"),
265 GetOperatorFactory().RegisterCreatorFunction(
266 OperatorKey(ePyramid, eIProductWRTDerivBase, eSumFac, false),
267 IProductWRTDerivBase_StdMat::create,
268 "IProductWRTDerivBase_SumFac_Pyr")};
269
270/**
271 * @brief Inner product operator using operator using matrix free operators.
272 */
273class IProductWRTDerivBase_MatrixFree final : virtual public Operator,
274 MatrixFreeBase,
275 IProductWRTDerivBase_Helper
276{
277public:
278 OPERATOR_CREATE(IProductWRTDerivBase_MatrixFree)
279
280 ~IProductWRTDerivBase_MatrixFree() final = default;
281
282 void operator()(const Array<OneD, const NekDouble> &entry0,
283 Array<OneD, NekDouble> &entry1,
284 Array<OneD, NekDouble> &entry2,
285 Array<OneD, NekDouble> &entry3,
286 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
287 {
288 Array<OneD, NekDouble> output;
289 Array<OneD, Array<OneD, NekDouble>> input(m_coordim);
290
291 // copy into padded vector
292 switch (m_coordim)
293 {
294 case 1:
295 input[0] = entry0;
296 output = entry1;
297 break;
298 case 2:
299 input[0] = entry0;
300 input[1] = entry1;
301 output = entry2;
302 break;
303 case 3:
304 input[0] = entry0;
305 input[1] = entry1;
306 input[2] = entry2;
307 output = entry3;
308 break;
309 default:
310 NEKERROR(ErrorUtil::efatal, "coordim not valid");
311 break;
312 }
313
314 (*m_oper)(input, output);
315 }
316
317 void operator()([[maybe_unused]] int dir,
318 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
319 [[maybe_unused]] Array<OneD, NekDouble> &output,
320 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
321 {
322 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
323 }
324
325private:
326 std::shared_ptr<MatrixFree::IProductWRTDerivBase> m_oper;
327 int m_coordim;
328
329 IProductWRTDerivBase_MatrixFree(
330 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
331 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
332 : Operator(pCollExp, pGeomData, factors),
333 MatrixFreeBase(pCollExp[0]->GetStdExp()->GetTotPoints(),
334 pCollExp[0]->GetStdExp()->GetNcoeffs(),
335 pCollExp.size()),
336 IProductWRTDerivBase_Helper()
337 {
338 const auto dim = pCollExp[0]->GetStdExp()->GetShapeDimension();
339 m_coordim = pCollExp[0]->GetCoordim();
340
341 // Basis vector
342 std::vector<LibUtilities::BasisSharedPtr> basis(dim);
343 for (unsigned int i = 0; i < dim; ++i)
344 {
345 basis[i] = pCollExp[0]->GetBasis(i);
346 }
347
348 // Get shape type
349 auto shapeType = pCollExp[0]->GetStdExp()->DetShapeType();
350
351 // Generate operator string and create operator.
352 std::string op_string = "IProductWRTDerivBase";
353 op_string += MatrixFree::GetOpstring(shapeType, m_isDeformed);
354 auto oper = MatrixFree::GetOperatorFactory().CreateInstance(
355 op_string, basis, pCollExp.size());
356
357 // set up required copies for operations
358 oper->SetUpBdata(basis);
359 oper->SetUpDBdata(basis);
360 oper->SetUpZW(basis);
361
362 // Set Jacobian
363 oper->SetJac(pGeomData->GetJacInterLeave(pCollExp, m_nElmtPad));
364
365 // Set derivative factors
366 oper->SetDF(pGeomData->GetDerivFactorsInterLeave(pCollExp, m_nElmtPad));
367
368 m_oper =
369 std::dynamic_pointer_cast<MatrixFree::IProductWRTDerivBase>(oper);
370 ASSERTL0(m_oper, "Failed to cast pointer.");
371 }
372};
373
374/// Factory initialisation for the IProductWRTDerivBase_MatrixFree operators
375OperatorKey IProductWRTDerivBase_MatrixFree::m_typeArr[] = {
376 GetOperatorFactory().RegisterCreatorFunction(
377 OperatorKey(eSegment, eIProductWRTDerivBase, eMatrixFree, false),
378 IProductWRTDerivBase_MatrixFree::create,
379 "IProductWRTDerivBase_MatrixFree_Seg"),
380 GetOperatorFactory().RegisterCreatorFunction(
381 OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eMatrixFree, false),
382 IProductWRTDerivBase_MatrixFree::create,
383 "IProductWRTDerivBase_MatrixFree_Quad"),
384 GetOperatorFactory().RegisterCreatorFunction(
385 OperatorKey(eTriangle, eIProductWRTDerivBase, eMatrixFree, false),
386 IProductWRTDerivBase_MatrixFree::create,
387 "IProductWRTDerivBase_MatrixFree_Tri"),
388 GetOperatorFactory().RegisterCreatorFunction(
389 OperatorKey(eHexahedron, eIProductWRTDerivBase, eMatrixFree, false),
390 IProductWRTDerivBase_MatrixFree::create,
391 "IProductWRTDerivBase_MatrixFree_Hex"),
392 GetOperatorFactory().RegisterCreatorFunction(
393 OperatorKey(ePrism, eIProductWRTDerivBase, eMatrixFree, false),
394 IProductWRTDerivBase_MatrixFree::create,
395 "IProductWRTDerivBase_MatrixFree_Prism"),
396 GetOperatorFactory().RegisterCreatorFunction(
397 OperatorKey(ePyramid, eIProductWRTDerivBase, eMatrixFree, false),
398 IProductWRTDerivBase_MatrixFree::create,
399 "IProductWRTDerivBase_MatrixFree_Pyr"),
400 GetOperatorFactory().RegisterCreatorFunction(
401 OperatorKey(eTetrahedron, eIProductWRTDerivBase, eMatrixFree, false),
402 IProductWRTDerivBase_MatrixFree::create,
403 "IProductWRTDerivBase_MatrixFree_Tet")};
404
405/**
406 * @brief Inner product WRT deriv base operator using element-wise operation
407 */
408class IProductWRTDerivBase_IterPerExp final : virtual public Operator,
409 IProductWRTDerivBase_Helper
410{
411public:
412 OPERATOR_CREATE(IProductWRTDerivBase_IterPerExp)
413
414 ~IProductWRTDerivBase_IterPerExp() final = default;
415
416 void operator()(const Array<OneD, const NekDouble> &entry0,
417 Array<OneD, NekDouble> &entry1,
418 Array<OneD, NekDouble> &entry2,
419 Array<OneD, NekDouble> &entry3,
420 Array<OneD, NekDouble> &wsp) final
421 {
422 unsigned int nPhys = m_stdExp->GetTotPoints();
423 unsigned int ntot = m_numElmt * nPhys;
424 unsigned int nmodes = m_stdExp->GetNcoeffs();
425 unsigned int nmax = max(ntot, m_numElmt * nmodes);
426 Array<OneD, Array<OneD, const NekDouble>> in(3);
427 Array<OneD, NekDouble> output, tmp1;
428 Array<OneD, Array<OneD, NekDouble>> tmp(3);
429
430 in[0] = entry0;
431 in[1] = entry1;
432 in[2] = entry2;
433
434 output = (m_coordim == 3) ? entry3 : (m_coordim == 2) ? entry2 : entry1;
435
436 for (int i = 0; i < m_dim; ++i)
437 {
438 tmp[i] = wsp + i * nmax;
439 }
440
441 // calculate Iproduct WRT Std Deriv
442 // first component
443 if (m_isDeformed)
444 {
445 // calculate dx/dxi in[0] + dy/dxi in[2] + dz/dxi in[3]
446 for (int i = 0; i < m_dim; ++i)
447 {
448 Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
449 for (int j = 1; j < m_coordim; ++j)
450 {
451 Vmath::Vvtvp(ntot, m_derivFac[i + j * m_dim], 1, in[j], 1,
452 tmp[i], 1, tmp[i], 1);
453 }
454 }
455
456 Vmath::Vmul(ntot, m_jac, 1, tmp[0], 1, tmp[0], 1);
457 }
458 else
459 {
460 Array<OneD, NekDouble> t;
461 for (int e = 0; e < m_numElmt; ++e)
462 {
463 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
464 for (int i = 0; i < m_dim; ++i)
465 {
466 Vmath::Smul(m_nqe, m_derivFac[i][e], in[0] + e * m_nqe, 1,
467 t = tmp[i] + e * m_nqe, 1);
468 for (int j = 1; j < m_coordim; ++j)
469 {
470 Vmath::Svtvp(m_nqe, m_derivFac[i + j * m_dim][e],
471 in[j] + e * m_nqe, 1, tmp[i] + e * m_nqe,
472 1, t = tmp[i] + e * m_nqe, 1);
473 }
474 }
475
476 Vmath::Smul(m_nqe, m_jac[e], tmp[0] + e * m_nqe, 1,
477 t = tmp[0] + e * m_nqe, 1);
478 }
479 }
480
481 for (int n = 0; n < m_numElmt; ++n)
482 {
483 m_stdExp->IProductWRTDerivBase(0, tmp[0] + n * nPhys,
484 tmp1 = output + n * nmodes);
485 }
486
487 // other components
488 for (int i = 1; i < m_dim; ++i)
489 {
490 // multiply by Jacobian
491 if (m_isDeformed)
492 {
493 Vmath::Vmul(ntot, m_jac, 1, tmp[i], 1, tmp[i], 1);
494 }
495 else
496 {
497 Array<OneD, NekDouble> t;
498 for (int e = 0; e < m_numElmt; ++e)
499 {
500 Vmath::Smul(m_nqe, m_jac[e], tmp[i] + e * m_nqe, 1,
501 t = tmp[i] + e * m_nqe, 1);
502 }
503 }
504
505 for (int n = 0; n < m_numElmt; ++n)
506 {
507 m_stdExp->IProductWRTDerivBase(i, tmp[i] + n * nPhys, tmp[0]);
508 Vmath::Vadd(nmodes, tmp[0], 1, output + n * nmodes, 1,
509 tmp1 = output + n * nmodes, 1);
510 }
511 }
512 }
513
514 void operator()([[maybe_unused]] int dir,
515 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
516 [[maybe_unused]] Array<OneD, NekDouble> &output,
517 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
518 {
519 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
520 }
521
522protected:
523 Array<TwoD, const NekDouble> m_derivFac;
524 Array<OneD, const NekDouble> m_jac;
525 int m_dim;
526 int m_coordim;
527
528private:
529 IProductWRTDerivBase_IterPerExp(
530 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
531 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
532 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper()
533 {
534 m_dim = pCollExp[0]->GetShapeDimension();
535 m_coordim = pCollExp[0]->GetCoordim();
536
537 m_nqe = m_stdExp->GetTotPoints();
538
539 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
540 m_jac = pGeomData->GetJac(pCollExp);
541 m_wspSize = m_dim * m_nqe * m_numElmt;
542 }
543};
544
545/// Factory initialisation for the IProductWRTDerivBase_IterPerExp operators
546OperatorKey IProductWRTDerivBase_IterPerExp::m_typeArr[] = {
547 GetOperatorFactory().RegisterCreatorFunction(
548 OperatorKey(eSegment, eIProductWRTDerivBase, eIterPerExp, false),
549 IProductWRTDerivBase_IterPerExp::create,
550 "IProductWRTDerivBase_IterPerExp_Seg"),
551 GetOperatorFactory().RegisterCreatorFunction(
552 OperatorKey(eTriangle, eIProductWRTDerivBase, eIterPerExp, false),
553 IProductWRTDerivBase_IterPerExp::create,
554 "IProductWRTDerivBase_IterPerExp_Tri"),
555 GetOperatorFactory().RegisterCreatorFunction(
556 OperatorKey(eTriangle, eIProductWRTDerivBase, eIterPerExp, true),
557 IProductWRTDerivBase_IterPerExp::create,
558 "IProductWRTDerivBase_IterPerExp_NodalTri"),
559 GetOperatorFactory().RegisterCreatorFunction(
560 OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eIterPerExp, false),
561 IProductWRTDerivBase_IterPerExp::create,
562 "IProductWRTDerivBase_IterPerExp_Quad"),
563 GetOperatorFactory().RegisterCreatorFunction(
564 OperatorKey(eTetrahedron, eIProductWRTDerivBase, eIterPerExp, false),
565 IProductWRTDerivBase_IterPerExp::create,
566 "IProductWRTDerivBase_IterPerExp_Tet"),
567 GetOperatorFactory().RegisterCreatorFunction(
568 OperatorKey(eTetrahedron, eIProductWRTDerivBase, eIterPerExp, true),
569 IProductWRTDerivBase_IterPerExp::create,
570 "IProductWRTDerivBase_IterPerExp_NodalTet"),
571 GetOperatorFactory().RegisterCreatorFunction(
572 OperatorKey(ePyramid, eIProductWRTDerivBase, eIterPerExp, false),
573 IProductWRTDerivBase_IterPerExp::create,
574 "IProductWRTDerivBase_IterPerExp_Pyr"),
575 GetOperatorFactory().RegisterCreatorFunction(
576 OperatorKey(ePrism, eIProductWRTDerivBase, eIterPerExp, false),
577 IProductWRTDerivBase_IterPerExp::create,
578 "IProductWRTDerivBase_IterPerExp_Prism"),
579 GetOperatorFactory().RegisterCreatorFunction(
580 OperatorKey(ePrism, eIProductWRTDerivBase, eIterPerExp, true),
581 IProductWRTDerivBase_IterPerExp::create,
582 "IProductWRTDerivBase_IterPerExp_NodalPrism"),
583 GetOperatorFactory().RegisterCreatorFunction(
584 OperatorKey(eHexahedron, eIProductWRTDerivBase, eIterPerExp, false),
585 IProductWRTDerivBase_IterPerExp::create,
586 "IProductWRTDerivBase_IterPerExp_Hex")};
587
588/**
589 * @brief Inner product WRT deriv base operator using LocalRegions
590 * implementation.
591 */
592class IProductWRTDerivBase_NoCollection final : virtual public Operator,
593 IProductWRTDerivBase_Helper
594{
595public:
596 OPERATOR_CREATE(IProductWRTDerivBase_NoCollection)
597
598 ~IProductWRTDerivBase_NoCollection() final = default;
599
600 void operator()(const Array<OneD, const NekDouble> &entry0,
601 Array<OneD, NekDouble> &entry1,
602 Array<OneD, NekDouble> &entry2,
603 Array<OneD, NekDouble> &entry3,
604 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
605 {
606 unsigned int nmodes = m_expList[0]->GetNcoeffs();
607 unsigned int nPhys = m_expList[0]->GetTotPoints();
608 Array<OneD, NekDouble> tmp(nmodes), tmp1;
609
610 Array<OneD, Array<OneD, const NekDouble>> in(3);
611 Array<OneD, NekDouble> output;
612 in[0] = entry0;
613 in[1] = entry1;
614 in[2] = entry2;
615
616 output = (m_coordim == 3) ? entry3 : (m_coordim == 2) ? entry2 : entry1;
617
618 for (int n = 0; n < m_numElmt; ++n)
619 {
620 m_expList[n]->IProductWRTDerivBase(0, in[0] + n * nPhys,
621 tmp1 = output + n * nmodes);
622 }
623
624 for (int i = 1; i < m_dim; ++i)
625 {
626 for (int n = 0; n < m_numElmt; ++n)
627 {
628 m_expList[n]->IProductWRTDerivBase(i, in[i] + n * nPhys, tmp);
629
630 Vmath::Vadd(nmodes, tmp, 1, output + n * nmodes, 1,
631 tmp1 = output + n * nmodes, 1);
632 }
633 }
634 }
635
636 void operator()([[maybe_unused]] int dir,
637 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
638 [[maybe_unused]] Array<OneD, NekDouble> &output,
639 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
640 {
641 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
642 }
643
644protected:
645 int m_dim;
646 int m_coordim;
647 vector<StdRegions::StdExpansionSharedPtr> m_expList;
648
649private:
650 IProductWRTDerivBase_NoCollection(
651 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
652 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
653 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper()
654 {
655 m_expList = pCollExp;
656 m_dim = pCollExp[0]->GetNumBases();
657 m_coordim = pCollExp[0]->GetCoordim();
658 }
659};
660
661/// Factory initialisation for the IProductWRTDerivBase_NoCollection operators
662OperatorKey IProductWRTDerivBase_NoCollection::m_typeArr[] = {
663 GetOperatorFactory().RegisterCreatorFunction(
664 OperatorKey(eSegment, eIProductWRTDerivBase, eNoCollection, false),
665 IProductWRTDerivBase_NoCollection::create,
666 "IProductWRTDerivBase_NoCollection_Seg"),
667 GetOperatorFactory().RegisterCreatorFunction(
668 OperatorKey(eTriangle, eIProductWRTDerivBase, eNoCollection, false),
669 IProductWRTDerivBase_NoCollection::create,
670 "IProductWRTDerivBase_NoCollection_Tri"),
671 GetOperatorFactory().RegisterCreatorFunction(
672 OperatorKey(eTriangle, eIProductWRTDerivBase, eNoCollection, true),
673 IProductWRTDerivBase_NoCollection::create,
674 "IProductWRTDerivBase_NoCollection_NodalTri"),
675 GetOperatorFactory().RegisterCreatorFunction(
676 OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eNoCollection,
677 false),
678 IProductWRTDerivBase_NoCollection::create,
679 "IProductWRTDerivBase_NoCollection_Quad"),
680 GetOperatorFactory().RegisterCreatorFunction(
681 OperatorKey(eTetrahedron, eIProductWRTDerivBase, eNoCollection, false),
682 IProductWRTDerivBase_NoCollection::create,
683 "IProductWRTDerivBase_NoCollection_Tet"),
684 GetOperatorFactory().RegisterCreatorFunction(
685 OperatorKey(eTetrahedron, eIProductWRTDerivBase, eNoCollection, true),
686 IProductWRTDerivBase_NoCollection::create,
687 "IProductWRTDerivBase_NoCollection_NodalTet"),
688 GetOperatorFactory().RegisterCreatorFunction(
689 OperatorKey(ePyramid, eIProductWRTDerivBase, eNoCollection, false),
690 IProductWRTDerivBase_NoCollection::create,
691 "IProductWRTDerivBase_NoCollection_Pyr"),
692 GetOperatorFactory().RegisterCreatorFunction(
693 OperatorKey(ePrism, eIProductWRTDerivBase, eNoCollection, false),
694 IProductWRTDerivBase_NoCollection::create,
695 "IProductWRTDerivBase_NoCollection_Prism"),
696 GetOperatorFactory().RegisterCreatorFunction(
697 OperatorKey(ePrism, eIProductWRTDerivBase, eNoCollection, true),
698 IProductWRTDerivBase_NoCollection::create,
699 "IProductWRTDerivBase_NoCollection_NodalPrism"),
700 GetOperatorFactory().RegisterCreatorFunction(
701 OperatorKey(eHexahedron, eIProductWRTDerivBase, eNoCollection, false),
702 IProductWRTDerivBase_NoCollection::create,
703 "IProductWRTDerivBase_NoCollection_Hex")};
704
705/**
706 * @brief Inner product WRT deriv base operator using sum-factorisation
707 * (Segment)
708 */
709class IProductWRTDerivBase_SumFac_Seg final : virtual public Operator,
710 IProductWRTDerivBase_Helper
711{
712public:
713 OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Seg)
714
715 ~IProductWRTDerivBase_SumFac_Seg() final = default;
716
717 void operator()(const Array<OneD, const NekDouble> &entry0,
718 Array<OneD, NekDouble> &entry1,
719 Array<OneD, NekDouble> &entry2,
720 Array<OneD, NekDouble> &entry3,
721 Array<OneD, NekDouble> &wsp) final
722 {
723 Array<OneD, Array<OneD, const NekDouble>> in(3);
724 Array<OneD, NekDouble> output;
725
726 output = (m_coordim == 3) ? entry3 : (m_coordim == 2) ? entry2 : entry1;
727
728 in[0] = entry0;
729 in[1] = entry1;
730 in[2] = entry2;
731
732 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
733 if (m_isDeformed)
734 {
735 Vmath::Vmul(m_wspSize, m_derivFac[0], 1, in[0], 1, wsp, 1);
736 for (int i = 1; i < m_coordim; ++i)
737 {
738 Vmath::Vvtvp(m_wspSize, m_derivFac[i], 1, in[i], 1, wsp, 1, wsp,
739 1);
740 }
741 }
742 else
743 {
744 Array<OneD, NekDouble> t;
745 for (int e = 0; e < m_numElmt; ++e)
746 {
747 Vmath::Smul(m_nquad0, m_derivFac[0][e], in[0] + e * m_nquad0, 1,
748 t = wsp + e * m_nquad0, 1);
749 }
750
751 for (int i = 1; i < m_coordim; ++i)
752 {
753 for (int e = 0; e < m_numElmt; ++e)
754 {
755 Vmath::Svtvp(m_nquad0, m_derivFac[i][e],
756 in[i] + e * m_nquad0, 1, wsp + e * m_nquad0, 1,
757 t = wsp + e * m_nquad0, 1);
758 }
759 }
760 }
761
762 Vmath::Vmul(m_wspSize, m_jacWStdW, 1, wsp, 1, wsp, 1);
763
764 // out = B0*in;
765 Blas::Dgemm('T', 'N', m_nmodes0, m_numElmt, m_nquad0, 1.0,
766 m_derbase0.get(), m_nquad0, &wsp[0], m_nquad0, 0.0,
767 &output[0], m_nmodes0);
768 }
769
770 void operator()([[maybe_unused]] int dir,
771 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
772 [[maybe_unused]] Array<OneD, NekDouble> &output,
773 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
774 {
775 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
776 }
777
778protected:
779 const int m_nquad0;
780 const int m_nmodes0;
781 Array<OneD, const NekDouble> m_jacWStdW;
782 Array<OneD, const NekDouble> m_derbase0;
783 Array<TwoD, const NekDouble> m_derivFac;
784 int m_coordim;
785
786private:
787 IProductWRTDerivBase_SumFac_Seg(
788 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
789 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
790 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper(),
791 m_nquad0(m_stdExp->GetNumPoints(0)),
792 m_nmodes0(m_stdExp->GetBasisNumModes(0)),
793 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata())
794 {
795 m_coordim = pCollExp[0]->GetCoordim();
796 m_wspSize = m_numElmt * m_nquad0;
797 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
798 m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
799 }
800};
801
802/// Factory initialisation for the IProductWRTDerivBase_SumFac_Seg operator
803OperatorKey IProductWRTDerivBase_SumFac_Seg::m_type =
804 GetOperatorFactory().RegisterCreatorFunction(
805 OperatorKey(eSegment, eIProductWRTDerivBase, eSumFac, false),
806 IProductWRTDerivBase_SumFac_Seg::create,
807 "IProductWRTDerivBase_SumFac_Seg");
808
809/**
810 * @brief Inner product WRT deriv base operator using sum-factorisation (Quad)
811 */
812class IProductWRTDerivBase_SumFac_Quad final : virtual public Operator,
813 IProductWRTDerivBase_Helper
814{
815public:
816 OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Quad)
817
818 ~IProductWRTDerivBase_SumFac_Quad() final = default;
819
820 void operator()(const Array<OneD, const NekDouble> &entry0,
821 Array<OneD, NekDouble> &entry1,
822 Array<OneD, NekDouble> &entry2,
823 Array<OneD, NekDouble> &entry3,
824 Array<OneD, NekDouble> &wsp) final
825 {
826 unsigned int nPhys = m_stdExp->GetTotPoints();
827 unsigned int ntot = m_numElmt * nPhys;
828 unsigned int nmodes = m_stdExp->GetNcoeffs();
829 unsigned int nmax = max(ntot, m_numElmt * nmodes);
830 Array<OneD, Array<OneD, const NekDouble>> in(3);
831 Array<OneD, NekDouble> output, wsp1;
832 Array<OneD, Array<OneD, NekDouble>> tmp(2);
833
834 in[0] = entry0;
835 in[1] = entry1;
836 in[2] = entry2;
837
838 output = (m_coordim == 2) ? entry2 : entry3;
839
840 tmp[0] = wsp;
841 tmp[1] = wsp + nmax;
842 wsp1 = wsp + 2 * nmax;
843
844 if (m_isDeformed)
845 {
846 // calculate dx/dxi in[0] + dy/dxi in[1]
847 for (int i = 0; i < 2; ++i)
848 {
849 Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
850 for (int j = 1; j < m_coordim; ++j)
851 {
852 Vmath::Vvtvp(ntot, m_derivFac[i + 2 * j], 1, in[j], 1,
853 tmp[i], 1, tmp[i], 1);
854 }
855 }
856 }
857 else
858 {
859 Array<OneD, NekDouble> t;
860 for (int e = 0; e < m_numElmt; ++e)
861 {
862 // calculate dx/dxi in[0] + dy/dxi in[1]
863 for (int i = 0; i < 2; ++i)
864 {
865 Vmath::Smul(m_nqe, m_derivFac[i][e], in[0] + e * m_nqe, 1,
866 t = tmp[i] + e * m_nqe, 1);
867 for (int j = 1; j < m_coordim; ++j)
868 {
869 Vmath::Svtvp(m_nqe, m_derivFac[i + 2 * j][e],
870 in[j] + e * m_nqe, 1, tmp[i] + e * m_nqe,
871 1, t = tmp[i] + e * m_nqe, 1);
872 }
873 }
874 }
875 }
876 // Iproduct wrt derivative of base 1
877 QuadIProduct(false, m_colldir1, m_numElmt, m_nquad0, m_nquad1,
878 m_nmodes0, m_nmodes1, m_derbase0, m_base1, m_jacWStdW,
879 tmp[0], output, wsp1);
880
881 // Iproduct wrt derivative of base 1
882 QuadIProduct(m_colldir0, false, m_numElmt, m_nquad0, m_nquad1,
883 m_nmodes0, m_nmodes1, m_base0, m_derbase1, m_jacWStdW,
884 tmp[1], tmp[0], wsp1);
885
886 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
887 }
888
889 void operator()([[maybe_unused]] int dir,
890 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
891 [[maybe_unused]] Array<OneD, NekDouble> &output,
892 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
893 {
894 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
895 }
896
897protected:
898 const int m_nquad0;
899 const int m_nquad1;
900 const int m_nmodes0;
901 const int m_nmodes1;
902 const bool m_colldir0;
903 const bool m_colldir1;
904 int m_coordim;
905 Array<TwoD, const NekDouble> m_derivFac;
906 Array<OneD, const NekDouble> m_jacWStdW;
907 Array<OneD, const NekDouble> m_base0;
908 Array<OneD, const NekDouble> m_base1;
909 Array<OneD, const NekDouble> m_derbase0;
910 Array<OneD, const NekDouble> m_derbase1;
911
912private:
913 IProductWRTDerivBase_SumFac_Quad(
914 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
915 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
916 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper(),
917 m_nquad0(m_stdExp->GetNumPoints(0)),
918 m_nquad1(m_stdExp->GetNumPoints(1)),
919 m_nmodes0(m_stdExp->GetBasisNumModes(0)),
920 m_nmodes1(m_stdExp->GetBasisNumModes(1)),
921 m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
922 m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
923 m_base0(m_stdExp->GetBasis(0)->GetBdata()),
924 m_base1(m_stdExp->GetBasis(1)->GetBdata()),
925 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
926 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata())
927 {
928 m_coordim = pCollExp[0]->GetCoordim();
929 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
930 m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
931 m_wspSize =
932 4 * m_numElmt * (max(m_nquad0 * m_nquad1, m_nmodes0 * m_nmodes1));
933 }
934};
935
936/// Factory initialisation for the IProductWRTDerivBase_SumFac_Quad operator
937OperatorKey IProductWRTDerivBase_SumFac_Quad::m_type =
938 GetOperatorFactory().RegisterCreatorFunction(
939 OperatorKey(eQuadrilateral, eIProductWRTDerivBase, eSumFac, false),
940 IProductWRTDerivBase_SumFac_Quad::create,
941 "IProductWRTDerivBase_IterPerExp_Quad");
942
943/**
944 * @brief Inner product WRT deriv base operator using sum-factorisation (Tri)
945 */
946class IProductWRTDerivBase_SumFac_Tri final : virtual public Operator,
947 IProductWRTDerivBase_Helper
948{
949public:
950 OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Tri)
951
952 ~IProductWRTDerivBase_SumFac_Tri() final = default;
953
954 /**
955 * This method calculates:
956 *
957 * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) \f]
958 *
959 * which can be represented in terms of local cartesian
960 * derivaties as:
961 *
962 * \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
963 * d\phi/d\xi_1\, d\xi_1/dx),in[0]) + \f]
964 *
965 * \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
966 * d\phi/d\xi_1\, d\xi_1/dy),in[1]) + \f]
967 *
968 * where we note that
969 *
970 * \f[ d\phi/d\xi_0 = d\phi/d\eta_0\, d\eta_0/d\xi_0 =
971 * d\phi/d\eta_0 2/(1-\eta_1) \f]
972 *
973 * \f[ d\phi/d\xi_1 = d\phi/d\eta_1\, d\eta_1/d\xi_1 +
974 * d\phi/d\eta_1\, d\eta_1/d\xi_1 = d\phi/d\eta_0 (1+\eta_0)/(1-\eta_1)
975 * + d\phi/d\eta_1 \f]
976 *
977 * and so the full inner products are
978 *
979 * \f[ (d\phi/dx,in[0]) + (dphi/dy,in[1]) =
980 * (d\phi/d\eta_0, ((2/(1-\eta_1) (d\xi_0/dx in[0] + d\xi_0/dy in[1])
981 * + (1-\eta_0)/(1-\eta_1) (d\xi_1/dx in[0]+d\xi_1/dy in[1]))
982 * + (d\phi/d\eta_1, (d\xi_1/dx in[0] + d\xi_1/dy in[1])) \f]
983 *
984 */
985 void operator()(const Array<OneD, const NekDouble> &entry0,
986 Array<OneD, NekDouble> &entry1,
987 Array<OneD, NekDouble> &entry2,
988 Array<OneD, NekDouble> &entry3,
989 Array<OneD, NekDouble> &wsp) final
990 {
991 unsigned int nPhys = m_stdExp->GetTotPoints();
992 unsigned int ntot = m_numElmt * nPhys;
993 unsigned int nmodes = m_stdExp->GetNcoeffs();
994 unsigned int nmax = max(ntot, m_numElmt * nmodes);
995 Array<OneD, Array<OneD, const NekDouble>> in(3);
996 Array<OneD, NekDouble> output, wsp1;
997 Array<OneD, Array<OneD, NekDouble>> tmp(2);
998
999 in[0] = entry0;
1000 in[1] = entry1;
1001 in[2] = entry2;
1002
1003 output = (m_coordim == 2) ? entry2 : entry3;
1004
1005 tmp[0] = wsp;
1006 tmp[1] = wsp + nmax;
1007 wsp1 = wsp + 2 * nmax;
1008
1009 if (m_isDeformed)
1010 {
1011 for (int i = 0; i < 2; ++i)
1012 {
1013 Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
1014 for (int j = 1; j < m_coordim; ++j)
1015 {
1016 Vmath::Vvtvp(ntot, m_derivFac[i + 2 * j], 1, in[j], 1,
1017 tmp[i], 1, tmp[i], 1);
1018 }
1019 }
1020 }
1021 else
1022 {
1023 Array<OneD, NekDouble> t;
1024 for (int e = 0; e < m_numElmt; ++e)
1025 {
1026 // calculate dx/dxi in[0] + dy/dxi in[1]
1027 for (int i = 0; i < 2; ++i)
1028 {
1029 Vmath::Smul(m_nqe, m_derivFac[i][e], in[0] + e * m_nqe, 1,
1030 t = tmp[i] + e * m_nqe, 1);
1031 for (int j = 1; j < m_coordim; ++j)
1032 {
1033 Vmath::Svtvp(m_nqe, m_derivFac[i + 2 * j][e],
1034 in[j] + e * m_nqe, 1, tmp[i] + e * m_nqe,
1035 1, t = tmp[i] + e * m_nqe, 1);
1036 }
1037 }
1038 }
1039 }
1040
1041 // Multiply by factor: 2/(1-z1)
1042 for (int i = 0; i < m_numElmt; ++i)
1043 {
1044 // scale tmp[0] by geometric factor: 2/(1-z1)
1045 Vmath::Vmul(nPhys, &m_fac0[0], 1, tmp[0].get() + i * nPhys, 1,
1046 tmp[0].get() + i * nPhys, 1);
1047
1048 // scale tmp[1] by geometric factor (1+z0)/(1-z1)
1049 Vmath::Vvtvp(nPhys, &m_fac1[0], 1, tmp[1].get() + i * nPhys, 1,
1050 tmp[0].get() + i * nPhys, 1, tmp[0].get() + i * nPhys,
1051 1);
1052 }
1053
1054 // Iproduct wrt derivative of base 0
1055 TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1, m_nmodes0,
1056 m_nmodes1, m_derbase0, m_base1, m_jacWStdW, tmp[0], output,
1057 wsp1);
1058
1059 // Iproduct wrt derivative of base 1
1060 TriIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1, m_nmodes0,
1061 m_nmodes1, m_base0, m_derbase1, m_jacWStdW, tmp[1], tmp[0],
1062 wsp1);
1063
1064 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1065 }
1066
1067 void operator()([[maybe_unused]] int dir,
1068 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
1069 [[maybe_unused]] Array<OneD, NekDouble> &output,
1070 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
1071 {
1072 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1073 }
1074
1075protected:
1076 const int m_nquad0;
1077 const int m_nquad1;
1078 const int m_nmodes0;
1079 const int m_nmodes1;
1080 const bool m_colldir0;
1081 const bool m_colldir1;
1082 int m_coordim;
1083 Array<TwoD, const NekDouble> m_derivFac;
1084 Array<OneD, const NekDouble> m_jacWStdW;
1085 Array<OneD, const NekDouble> m_base0;
1086 Array<OneD, const NekDouble> m_base1;
1087 Array<OneD, const NekDouble> m_derbase0;
1088 Array<OneD, const NekDouble> m_derbase1;
1089 Array<OneD, NekDouble> m_fac0;
1090 Array<OneD, NekDouble> m_fac1;
1091 bool m_sortTopVertex;
1092
1093private:
1094 IProductWRTDerivBase_SumFac_Tri(
1095 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1096 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
1097 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper(),
1098 m_nquad0(m_stdExp->GetNumPoints(0)),
1099 m_nquad1(m_stdExp->GetNumPoints(1)),
1100 m_nmodes0(m_stdExp->GetBasisNumModes(0)),
1101 m_nmodes1(m_stdExp->GetBasisNumModes(1)),
1102 m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
1103 m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
1104 m_base0(m_stdExp->GetBasis(0)->GetBdata()),
1105 m_base1(m_stdExp->GetBasis(1)->GetBdata()),
1106 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1107 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata())
1108 {
1109 m_coordim = pCollExp[0]->GetCoordim();
1110 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1111 m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1112 m_wspSize =
1113 4 * m_numElmt * (max(m_nquad0 * m_nquad1, m_nmodes0 * m_nmodes1));
1114
1115 if (m_stdExp->GetBasis(0)->GetBasisType() == LibUtilities::eModified_A)
1116 {
1117 m_sortTopVertex = true;
1118 }
1119 else
1120 {
1121 m_sortTopVertex = false;
1122 }
1123
1124 const Array<OneD, const NekDouble> &z0 = m_stdExp->GetBasis(0)->GetZ();
1125 const Array<OneD, const NekDouble> &z1 = m_stdExp->GetBasis(1)->GetZ();
1126
1127 m_fac0 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1);
1128 // set up geometric factor: 2/(1-z1)
1129 for (int i = 0; i < m_nquad0; ++i)
1130 {
1131 for (int j = 0; j < m_nquad1; ++j)
1132 {
1133 m_fac0[i + j * m_nquad0] = 2.0 / (1 - z1[j]);
1134 }
1135 }
1136
1137 m_fac1 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1);
1138 // set up geometric factor: (1+z0)/(1-z1)
1139 for (int i = 0; i < m_nquad0; ++i)
1140 {
1141 for (int j = 0; j < m_nquad1; ++j)
1142 {
1143 m_fac1[i + j * m_nquad0] = (1 + z0[i]) / (1 - z1[j]);
1144 }
1145 }
1146 }
1147};
1148
1149/// Factory initialisation for the IProductWRTDerivBase_SumFac_Tri operator
1150OperatorKey IProductWRTDerivBase_SumFac_Tri::m_type =
1151 GetOperatorFactory().RegisterCreatorFunction(
1152 OperatorKey(eTriangle, eIProductWRTDerivBase, eSumFac, false),
1153 IProductWRTDerivBase_SumFac_Tri::create,
1154 "IProductWRTDerivBase_IterPerExp_Tri");
1155
1156/**
1157 * @brief Inner product WRT deriv base operator using sum-factorisation (Hex)
1158 */
1159class IProductWRTDerivBase_SumFac_Hex final : virtual public Operator,
1160 IProductWRTDerivBase_Helper
1161{
1162public:
1163 OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Hex)
1164
1165 ~IProductWRTDerivBase_SumFac_Hex() final = default;
1166
1167 void operator()(const Array<OneD, const NekDouble> &entry0,
1168 Array<OneD, NekDouble> &entry1,
1169 Array<OneD, NekDouble> &entry2,
1170 Array<OneD, NekDouble> &entry3,
1171 Array<OneD, NekDouble> &wsp) final
1172 {
1173 unsigned int nPhys = m_stdExp->GetTotPoints();
1174 unsigned int ntot = m_numElmt * nPhys;
1175 unsigned int nmodes = m_stdExp->GetNcoeffs();
1176 unsigned int nmax = max(ntot, m_numElmt * nmodes);
1177 Array<OneD, Array<OneD, const NekDouble>> in(3);
1178 Array<OneD, NekDouble> output, wsp1;
1179 Array<OneD, Array<OneD, NekDouble>> tmp(3);
1180
1181 in[0] = entry0;
1182 in[1] = entry1;
1183 in[2] = entry2;
1184
1185 output = entry3;
1186
1187 for (int i = 0; i < 3; ++i)
1188 {
1189 tmp[i] = wsp + i * nmax;
1190 }
1191
1192 if (m_isDeformed)
1193 {
1194 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1195 for (int i = 0; i < 3; ++i)
1196 {
1197 Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
1198 for (int j = 1; j < 3; ++j)
1199 {
1200 Vmath::Vvtvp(ntot, m_derivFac[i + 3 * j], 1, in[j], 1,
1201 tmp[i], 1, tmp[i], 1);
1202 }
1203 }
1204 }
1205 else
1206 {
1207 Array<OneD, NekDouble> t;
1208 for (int e = 0; e < m_numElmt; ++e)
1209 {
1210 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1211 for (int i = 0; i < 3; ++i)
1212 {
1213 Vmath::Smul(m_nqe, m_derivFac[i][e], in[0] + e * m_nqe, 1,
1214 t = tmp[i] + e * m_nqe, 1);
1215 for (int j = 1; j < 3; ++j)
1216 {
1217 Vmath::Svtvp(m_nqe, m_derivFac[i + 3 * j][e],
1218 in[j] + e * m_nqe, 1, tmp[i] + e * m_nqe,
1219 1, t = tmp[i] + e * m_nqe, 1);
1220 }
1221 }
1222 }
1223 }
1224
1225 wsp1 = wsp + 3 * nmax;
1226
1227 // calculate Iproduct WRT Std Deriv
1228 HexIProduct(false, m_colldir1, m_colldir2, m_numElmt, m_nquad0,
1229 m_nquad1, m_nquad2, m_nmodes0, m_nmodes1, m_nmodes2,
1230 m_derbase0, m_base1, m_base2, m_jacWStdW, tmp[0], output,
1231 wsp1);
1232
1233 HexIProduct(m_colldir0, false, m_colldir2, m_numElmt, m_nquad0,
1234 m_nquad1, m_nquad2, m_nmodes0, m_nmodes1, m_nmodes2,
1235 m_base0, m_derbase1, m_base2, m_jacWStdW, tmp[1], tmp[0],
1236 wsp1);
1237 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1238
1239 HexIProduct(m_colldir0, m_colldir1, false, m_numElmt, m_nquad0,
1240 m_nquad1, m_nquad2, m_nmodes0, m_nmodes1, m_nmodes2,
1241 m_base0, m_base1, m_derbase2, m_jacWStdW, tmp[2], tmp[0],
1242 wsp1);
1243 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1244 }
1245
1246 void operator()([[maybe_unused]] int dir,
1247 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
1248 [[maybe_unused]] Array<OneD, NekDouble> &output,
1249 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
1250 {
1251 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1252 }
1253
1254protected:
1255 const int m_nquad0;
1256 const int m_nquad1;
1257 const int m_nquad2;
1258 const int m_nmodes0;
1259 const int m_nmodes1;
1260 const int m_nmodes2;
1261 const bool m_colldir0;
1262 const bool m_colldir1;
1263 const bool m_colldir2;
1264 Array<OneD, const NekDouble> m_jacWStdW;
1265 Array<OneD, const NekDouble> m_base0;
1266 Array<OneD, const NekDouble> m_base1;
1267 Array<OneD, const NekDouble> m_base2;
1268 Array<OneD, const NekDouble> m_derbase0;
1269 Array<OneD, const NekDouble> m_derbase1;
1270 Array<OneD, const NekDouble> m_derbase2;
1271 Array<TwoD, const NekDouble> m_derivFac;
1272
1273private:
1274 IProductWRTDerivBase_SumFac_Hex(
1275 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1276 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
1277 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper(),
1278 m_nquad0(m_stdExp->GetNumPoints(0)),
1279 m_nquad1(m_stdExp->GetNumPoints(1)),
1280 m_nquad2(m_stdExp->GetNumPoints(2)),
1281 m_nmodes0(m_stdExp->GetBasisNumModes(0)),
1282 m_nmodes1(m_stdExp->GetBasisNumModes(1)),
1283 m_nmodes2(m_stdExp->GetBasisNumModes(2)),
1284 m_colldir0(m_stdExp->GetBasis(0)->Collocation()),
1285 m_colldir1(m_stdExp->GetBasis(1)->Collocation()),
1286 m_colldir2(m_stdExp->GetBasis(2)->Collocation()),
1287 m_base0(m_stdExp->GetBasis(0)->GetBdata()),
1288 m_base1(m_stdExp->GetBasis(1)->GetBdata()),
1289 m_base2(m_stdExp->GetBasis(2)->GetBdata()),
1290 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1291 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1292 m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1293
1294 {
1295 m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1296 m_wspSize = 6 * m_numElmt *
1297 (max(m_nquad0 * m_nquad1 * m_nquad2,
1298 m_nmodes0 * m_nmodes1 * m_nmodes2));
1299 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1300 }
1301};
1302
1303/// Factory initialisation for the IProductWRTDerivBase_SumFac_Hex operator
1304OperatorKey IProductWRTDerivBase_SumFac_Hex::m_type =
1305 GetOperatorFactory().RegisterCreatorFunction(
1306 OperatorKey(eHexahedron, eIProductWRTDerivBase, eSumFac, false),
1307 IProductWRTDerivBase_SumFac_Hex::create,
1308 "IProductWRTDerivBase_SumFac_Hex");
1309
1310/**
1311 * @brief Inner product WRT deriv base operator using sum-factorisation (Tet)
1312 */
1313class IProductWRTDerivBase_SumFac_Tet : virtual public Operator,
1314 IProductWRTDerivBase_Helper
1315{
1316public:
1317 OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Tet)
1318
1319 /**
1320 * This method calculates:
1321 *
1322 * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
1323 *
1324 * which can be represented in terms of local cartesian
1325 * derivaties as:
1326 *
1327 * \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
1328 * d\phi/d\xi_1\, d\xi_1/dx +
1329 * d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
1330 *
1331 * \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
1332 * d\phi/d\xi_1\, d\xi_1/dy +
1333 * d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
1334 *
1335 * \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
1336 * d\phi/d\xi_1\, d\xi_1/dz +
1337 * d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
1338 *
1339 * where we note that
1340 *
1341 * \f[ d\phi/d\xi_0 = d\phi/d\eta_0 4/((1-\eta_1)(1-\eta_2)) \f]
1342 *
1343 * \f[ d\phi/d\xi_1 = d\phi/d\eta_0 2(1+\eta_0)/((1-\eta_1)(1-\eta_2))
1344 * + d\phi/d\eta_1 2/(1-\eta_2) \f]
1345 *
1346 * \f[ d\phi/d\xi_2 = d\phi/d\eta_0 2(1+\eta_0)/((1-\eta_1)(1-\eta_2))
1347 * + d\phi/d\eta_1 (1+\eta_1)/(1-\eta_2) + d\phi/d\eta_2 \f]
1348 *
1349 * and so the full inner products are
1350 *
1351 * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
1352 *
1353 * \f[ (d\phi/d\eta_0, fac0 (tmp0 + fac1(tmp1 + tmp2)))
1354 * + (d\phi/d\eta_1, fac2 (tmp1 + fac3 tmp2))
1355 * + (d\phi/d\eta_2, tmp2) \f]
1356 *
1357 * where
1358 *
1359 * \f[ \begin{array}{lcl}
1360 * tmp0 &=& (d\xi_0/dx in[0] + d\xi_0/dy in[1] + d\xi_0/dz in[2]) \\
1361 * tmp1 &=& (d\xi_1/dx in[0] + d\xi_1/dy in[1] + d\xi_1/dz in[2]) \\
1362 * tmp2 &=& (d\xi_2/dx in[0] + d\xi_2/dy in[1] + d\xi_2/dz in[2])
1363 * \end{array} \f]
1364 *
1365 * \f[ \begin{array}{lcl}
1366 * fac0 &= & 4/((1-\eta_1)(1-\eta_2)) \\
1367 * fac1 &= & (1+\eta_0)/2 \\
1368 * fac2 &= & 2/(1-\eta_2) \\
1369 * fac3 &= & (1+\eta_1)/2 \end{array} \f]
1370 *
1371 */
1372 void operator()(const Array<OneD, const NekDouble> &entry0,
1373 Array<OneD, NekDouble> &entry1,
1374 Array<OneD, NekDouble> &entry2,
1375 Array<OneD, NekDouble> &entry3,
1376 Array<OneD, NekDouble> &wsp) final
1377 {
1378 unsigned int nPhys = m_stdExp->GetTotPoints();
1379 unsigned int ntot = m_numElmt * nPhys;
1380 unsigned int nmodes = m_stdExp->GetNcoeffs();
1381 unsigned int nmax = max(ntot, m_numElmt * nmodes);
1382 Array<OneD, Array<OneD, const NekDouble>> in(3);
1383 Array<OneD, NekDouble> output, wsp1;
1384 Array<OneD, Array<OneD, NekDouble>> tmp(3);
1385
1386 in[0] = entry0;
1387 in[1] = entry1;
1388 in[2] = entry2;
1389
1390 output = entry3;
1391
1392 for (int i = 0; i < 3; ++i)
1393 {
1394 tmp[i] = wsp + i * nmax;
1395 }
1396
1397 if (m_isDeformed)
1398 {
1399 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1400 for (int i = 0; i < 3; ++i)
1401 {
1402 Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
1403 for (int j = 1; j < 3; ++j)
1404 {
1405 Vmath::Vvtvp(ntot, m_derivFac[i + 3 * j], 1, in[j], 1,
1406 tmp[i], 1, tmp[i], 1);
1407 }
1408 }
1409 }
1410 else
1411 {
1412 Array<OneD, NekDouble> t;
1413 for (int e = 0; e < m_numElmt; ++e)
1414 {
1415 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1416 for (int i = 0; i < 3; ++i)
1417 {
1418 Vmath::Smul(m_nqe, m_derivFac[i][e], in[0] + e * m_nqe, 1,
1419 t = tmp[i] + e * m_nqe, 1);
1420 for (int j = 1; j < 3; ++j)
1421 {
1422 Vmath::Svtvp(m_nqe, m_derivFac[i + 3 * j][e],
1423 in[j] + e * m_nqe, 1, tmp[i] + e * m_nqe,
1424 1, t = tmp[i] + e * m_nqe, 1);
1425 }
1426 }
1427 }
1428 }
1429
1430 wsp1 = wsp + 3 * nmax;
1431
1432 // Sort into eta factors
1433 for (int i = 0; i < m_numElmt; ++i)
1434 {
1435 // add tmp[1] + tmp[2]
1436 Vmath::Vadd(nPhys, tmp[1].get() + i * nPhys, 1,
1437 tmp[2].get() + i * nPhys, 1, wsp1.get(), 1);
1438
1439 // scale wsp1 by fac1 and add to tmp0
1440 Vmath::Vvtvp(nPhys, &m_fac1[0], 1, wsp1.get(), 1,
1441 tmp[0].get() + i * nPhys, 1, tmp[0].get() + i * nPhys,
1442 1);
1443
1444 // scale tmp[0] by fac0
1445 Vmath::Vmul(nPhys, &m_fac0[0], 1, tmp[0].get() + i * nPhys, 1,
1446 tmp[0].get() + i * nPhys, 1);
1447
1448 // scale tmp[2] by fac3 and add to tmp1
1449 Vmath::Vvtvp(nPhys, &m_fac3[0], 1, tmp[2].get() + i * nPhys, 1,
1450 tmp[1].get() + i * nPhys, 1, tmp[1].get() + i * nPhys,
1451 1);
1452
1453 // scale tmp[1] by fac2
1454 Vmath::Vmul(nPhys, &m_fac2[0], 1, tmp[1].get() + i * nPhys, 1,
1455 tmp[1].get() + i * nPhys, 1);
1456 }
1457
1458 // calculate Iproduct WRT Std Deriv
1459 TetIProduct(m_sortTopEdge, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1460 m_nmodes0, m_nmodes1, m_nmodes2, m_derbase0, m_base1,
1461 m_base2, m_jacWStdW, tmp[0], output, wsp1);
1462
1463 TetIProduct(m_sortTopEdge, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1464 m_nmodes0, m_nmodes1, m_nmodes2, m_base0, m_derbase1,
1465 m_base2, m_jacWStdW, tmp[1], tmp[0], wsp1);
1466 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1467
1468 TetIProduct(m_sortTopEdge, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1469 m_nmodes0, m_nmodes1, m_nmodes2, m_base0, m_base1,
1470 m_derbase2, m_jacWStdW, tmp[2], tmp[0], wsp1);
1471 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1472 }
1473
1474 void operator()([[maybe_unused]] int dir,
1475 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
1476 [[maybe_unused]] Array<OneD, NekDouble> &output,
1477 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
1478 {
1479 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1480 }
1481
1482protected:
1483 const int m_nquad0;
1484 const int m_nquad1;
1485 const int m_nquad2;
1486 const int m_nmodes0;
1487 const int m_nmodes1;
1488 const int m_nmodes2;
1489 Array<OneD, const NekDouble> m_jacWStdW;
1490 Array<OneD, const NekDouble> m_base0;
1491 Array<OneD, const NekDouble> m_base1;
1492 Array<OneD, const NekDouble> m_base2;
1493 Array<OneD, const NekDouble> m_derbase0;
1494 Array<OneD, const NekDouble> m_derbase1;
1495 Array<OneD, const NekDouble> m_derbase2;
1496 Array<TwoD, const NekDouble> m_derivFac;
1497 Array<OneD, NekDouble> m_fac0;
1498 Array<OneD, NekDouble> m_fac1;
1499 Array<OneD, NekDouble> m_fac2;
1500 Array<OneD, NekDouble> m_fac3;
1501 bool m_sortTopEdge;
1502
1503private:
1504 IProductWRTDerivBase_SumFac_Tet(
1505 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1506 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
1507 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper(),
1508 m_nquad0(m_stdExp->GetNumPoints(0)),
1509 m_nquad1(m_stdExp->GetNumPoints(1)),
1510 m_nquad2(m_stdExp->GetNumPoints(2)),
1511 m_nmodes0(m_stdExp->GetBasisNumModes(0)),
1512 m_nmodes1(m_stdExp->GetBasisNumModes(1)),
1513 m_nmodes2(m_stdExp->GetBasisNumModes(2)),
1514 m_base0(m_stdExp->GetBasis(0)->GetBdata()),
1515 m_base1(m_stdExp->GetBasis(1)->GetBdata()),
1516 m_base2(m_stdExp->GetBasis(2)->GetBdata()),
1517 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1518 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1519 m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1520
1521 {
1522 m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1523 m_wspSize = 6 * m_numElmt *
1524 (max(m_nquad0 * m_nquad1 * m_nquad2,
1525 m_nmodes0 * m_nmodes1 * m_nmodes2));
1526 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1527
1528 const Array<OneD, const NekDouble> &z0 = m_stdExp->GetBasis(0)->GetZ();
1529 const Array<OneD, const NekDouble> &z1 = m_stdExp->GetBasis(1)->GetZ();
1530 const Array<OneD, const NekDouble> &z2 = m_stdExp->GetBasis(2)->GetZ();
1531
1532 m_fac0 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
1533 m_fac1 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
1534 m_fac2 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
1535 m_fac3 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
1536 // calculate 2.0/((1-eta_1)(1-eta_2))
1537 for (int i = 0; i < m_nquad0; ++i)
1538 {
1539 for (int j = 0; j < m_nquad1; ++j)
1540 {
1541 for (int k = 0; k < m_nquad2; ++k)
1542 {
1543 m_fac0[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
1544 4.0 / ((1 - z1[j]) * (1 - z2[k]));
1545 m_fac1[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
1546 (1 + z0[i]) * 0.5;
1547 m_fac2[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
1548 2.0 / (1 - z2[k]);
1549 m_fac3[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
1550 (1 + z1[j]) * 0.5;
1551 }
1552 }
1553 }
1554
1555 if (m_stdExp->GetBasis(0)->GetBasisType() == LibUtilities::eModified_A)
1556 {
1557 m_sortTopEdge = true;
1558 }
1559 else
1560 {
1561 m_sortTopEdge = false;
1562 }
1563 }
1564};
1565
1566/// Factory initialisation for the IProductWRTDerivBase_SumFac_Tet operator
1567OperatorKey IProductWRTDerivBase_SumFac_Tet::m_type =
1568 GetOperatorFactory().RegisterCreatorFunction(
1569 OperatorKey(eTetrahedron, eIProductWRTDerivBase, eSumFac, false),
1570 IProductWRTDerivBase_SumFac_Tet::create,
1571 "IProductWRTDerivBase_SumFac_Tet");
1572
1573/**
1574 * @brief Inner product WRT deriv base operator using sum-factorisation (Prism)
1575 */
1576class IProductWRTDerivBase_SumFac_Prism final : virtual public Operator,
1577 IProductWRTDerivBase_Helper
1578{
1579public:
1580 OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Prism)
1581
1582 ~IProductWRTDerivBase_SumFac_Prism() final = default;
1583
1584 /**
1585 * This method calculates:
1586 *
1587 * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
1588 *
1589 * which can be represented in terms of local cartesian
1590 * derivaties as:
1591 *
1592 * \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
1593 * d\phi/d\xi_1\, d\xi_1/dx +
1594 * d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
1595 *
1596 * \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
1597 * d\phi/d\xi_1\, d\xi_1/dy +
1598 * d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
1599 *
1600 * \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
1601 * d\phi/d\xi_1\, d\xi_1/dz +
1602 * d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
1603 *
1604 * where we note that
1605 *
1606 * \f[ d\phi/d\xi_0 =
1607 * d\phi/d\eta_0 d\eta_0/d\xi_0 = d\phi/d\eta_0 2/(1-\eta_2) \f]
1608 *
1609 * \f[ d\phi/d\xi_2 =
1610 * d\phi/d\eta_0 d\eta_0/d\xi_2 + d\phi/d\eta_2 d\eta_2/d\xi_2 =
1611 * d\phi/d\eta_0 (1+\eta_0)/(1-\eta_2) + d\phi/d\eta_2 \f]
1612 *
1613 *
1614 * and so the full inner products are
1615 *
1616 * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
1617 *
1618 * \f[ (d\phi/d\eta_0, ((2/(1-\eta_2) (d\xi_0/dx in[0] + d\xi_0/dy in[1]
1619 * + d\xi_0/dz in[2])
1620 * + (1-\eta_0)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1]
1621 * + d\xi_2/dz in[2] )) + \f]
1622 *
1623 * \f[ (d\phi/d\eta_1, (d\xi_1/dx in[0] + d\xi_1/dy in[1]
1624 * + d\xi_1/dz in[2])) + \f]
1625 *
1626 * \f[ (d\phi/d\eta_2, (d\xi_2/dx in[0] + d\xi_2/dy in[1]
1627 * + d\xi_2/dz in[2])) \f]
1628 *
1629 */
1630 void operator()(const Array<OneD, const NekDouble> &entry0,
1631 Array<OneD, NekDouble> &entry1,
1632 Array<OneD, NekDouble> &entry2,
1633 Array<OneD, NekDouble> &entry3,
1634 Array<OneD, NekDouble> &wsp) final
1635 {
1636 unsigned int nPhys = m_stdExp->GetTotPoints();
1637 unsigned int ntot = m_numElmt * nPhys;
1638 unsigned int nmodes = m_stdExp->GetNcoeffs();
1639 unsigned int nmax = max(ntot, m_numElmt * nmodes);
1640 Array<OneD, Array<OneD, const NekDouble>> in(3);
1641 Array<OneD, NekDouble> output, wsp1;
1642 Array<OneD, Array<OneD, NekDouble>> tmp(3);
1643
1644 in[0] = entry0;
1645 in[1] = entry1;
1646 in[2] = entry2;
1647
1648 output = entry3;
1649
1650 for (int i = 0; i < 3; ++i)
1651 {
1652 tmp[i] = wsp + i * nmax;
1653 }
1654
1655 if (m_isDeformed)
1656 {
1657 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1658 for (int i = 0; i < 3; ++i)
1659 {
1660 Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
1661 for (int j = 1; j < 3; ++j)
1662 {
1663 Vmath::Vvtvp(ntot, m_derivFac[i + 3 * j], 1, in[j], 1,
1664 tmp[i], 1, tmp[i], 1);
1665 }
1666 }
1667 }
1668 else
1669 {
1670 Array<OneD, NekDouble> t;
1671 for (int e = 0; e < m_numElmt; ++e)
1672 {
1673 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1674 for (int i = 0; i < 3; ++i)
1675 {
1676 Vmath::Smul(m_nqe, m_derivFac[i][e], in[0] + e * m_nqe, 1,
1677 t = tmp[i] + e * m_nqe, 1);
1678 for (int j = 1; j < 3; ++j)
1679 {
1680 Vmath::Svtvp(m_nqe, m_derivFac[i + 3 * j][e],
1681 in[j] + e * m_nqe, 1, tmp[i] + e * m_nqe,
1682 1, t = tmp[i] + e * m_nqe, 1);
1683 }
1684 }
1685 }
1686 }
1687
1688 wsp1 = wsp + 3 * nmax;
1689
1690 // Sort into eta factors
1691 for (int i = 0; i < m_numElmt; ++i)
1692 {
1693 // scale tmp[0] by fac0
1694 Vmath::Vmul(nPhys, &m_fac0[0], 1, tmp[0].get() + i * nPhys, 1,
1695 tmp[0].get() + i * nPhys, 1);
1696
1697 // scale tmp[2] by fac1 and add to tmp0
1698 Vmath::Vvtvp(nPhys, &m_fac1[0], 1, tmp[2].get() + i * nPhys, 1,
1699 tmp[0].get() + i * nPhys, 1, tmp[0].get() + i * nPhys,
1700 1);
1701 }
1702
1703 // calculate Iproduct WRT Std Deriv
1704 PrismIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1705 m_nmodes0, m_nmodes1, m_nmodes2, m_derbase0, m_base1,
1706 m_base2, m_jacWStdW, tmp[0], output, wsp1);
1707
1708 PrismIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1709 m_nmodes0, m_nmodes1, m_nmodes2, m_base0, m_derbase1,
1710 m_base2, m_jacWStdW, tmp[1], tmp[0], wsp1);
1711 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1712
1713 PrismIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1714 m_nmodes0, m_nmodes1, m_nmodes2, m_base0, m_base1,
1715 m_derbase2, m_jacWStdW, tmp[2], tmp[0], wsp1);
1716 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1717 }
1718
1719 void operator()([[maybe_unused]] int dir,
1720 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
1721 [[maybe_unused]] Array<OneD, NekDouble> &output,
1722 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
1723 {
1724 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1725 }
1726
1727protected:
1728 const int m_nquad0;
1729 const int m_nquad1;
1730 const int m_nquad2;
1731 const int m_nmodes0;
1732 const int m_nmodes1;
1733 const int m_nmodes2;
1734 Array<OneD, const NekDouble> m_jacWStdW;
1735 Array<OneD, const NekDouble> m_base0;
1736 Array<OneD, const NekDouble> m_base1;
1737 Array<OneD, const NekDouble> m_base2;
1738 Array<OneD, const NekDouble> m_derbase0;
1739 Array<OneD, const NekDouble> m_derbase1;
1740 Array<OneD, const NekDouble> m_derbase2;
1741 Array<TwoD, const NekDouble> m_derivFac;
1742 Array<OneD, NekDouble> m_fac0;
1743 Array<OneD, NekDouble> m_fac1;
1744 bool m_sortTopVertex;
1745
1746private:
1747 IProductWRTDerivBase_SumFac_Prism(
1748 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
1749 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
1750 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper(),
1751 m_nquad0(m_stdExp->GetNumPoints(0)),
1752 m_nquad1(m_stdExp->GetNumPoints(1)),
1753 m_nquad2(m_stdExp->GetNumPoints(2)),
1754 m_nmodes0(m_stdExp->GetBasisNumModes(0)),
1755 m_nmodes1(m_stdExp->GetBasisNumModes(1)),
1756 m_nmodes2(m_stdExp->GetBasisNumModes(2)),
1757 m_base0(m_stdExp->GetBasis(0)->GetBdata()),
1758 m_base1(m_stdExp->GetBasis(1)->GetBdata()),
1759 m_base2(m_stdExp->GetBasis(2)->GetBdata()),
1760 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
1761 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
1762 m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
1763
1764 {
1765 m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
1766 m_wspSize = 6 * m_numElmt *
1767 (max(m_nquad0 * m_nquad1 * m_nquad2,
1768 m_nmodes0 * m_nmodes1 * m_nmodes2));
1769 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
1770
1771 if (m_stdExp->GetBasis(0)->GetBasisType() == LibUtilities::eModified_A)
1772 {
1773 m_sortTopVertex = true;
1774 }
1775 else
1776 {
1777 m_sortTopVertex = false;
1778 }
1779
1780 const Array<OneD, const NekDouble> &z0 = m_stdExp->GetBasis(0)->GetZ();
1781 const Array<OneD, const NekDouble> &z2 = m_stdExp->GetBasis(2)->GetZ();
1782
1783 m_fac0 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
1784 m_fac1 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
1785
1786 for (int i = 0; i < m_nquad0; ++i)
1787 {
1788 for (int j = 0; j < m_nquad1; ++j)
1789 {
1790 for (int k = 0; k < m_nquad2; ++k)
1791 {
1792 // set up geometric factor: 2/(1-z1)
1793 m_fac0[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
1794 2.0 / (1 - z2[k]);
1795 // set up geometric factor: (1+z0)/(1-z1)
1796 m_fac1[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
1797 (1 + z0[i]) / (1 - z2[k]);
1798 }
1799 }
1800 }
1801 }
1802};
1803
1804/// Factory initialisation for the IProductWRTDerivBase_SumFac_Prism operator
1805OperatorKey IProductWRTDerivBase_SumFac_Prism::m_type =
1806 GetOperatorFactory().RegisterCreatorFunction(
1807 OperatorKey(ePrism, eIProductWRTDerivBase, eSumFac, false),
1808 IProductWRTDerivBase_SumFac_Prism::create,
1809 "IProductWRTDerivBase_SumFac_Prism");
1810
1811/**
1812 * @brief Inner product WRT deriv base operator using sum-factorisation (Pyr)
1813 */
1814class IProductWRTDerivBase_SumFac_Pyr final : virtual public Operator,
1815 IProductWRTDerivBase_Helper
1816{
1817public:
1818 OPERATOR_CREATE(IProductWRTDerivBase_SumFac_Pyr)
1819
1820 ~IProductWRTDerivBase_SumFac_Pyr() final = default;
1821
1822 /**
1823 * This method calculates:
1824 *
1825 * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) \f]
1826 *
1827 * which can be represented in terms of local cartesian
1828 * derivaties as:
1829 *
1830 * \f[ ((d\phi/d\xi_0\, d\xi_0/dx +
1831 * d\phi/d\xi_1\, d\xi_1/dx +
1832 * d\phi/d\xi_2\, d\xi_2/dx),in[0]) + \f]
1833 *
1834 * \f[ ((d\phi/d\xi_0\, d\xi_0/dy +
1835 * d\phi/d\xi_1\, d\xi_1/dy +
1836 * d\phi/d\xi_2\, d\xi_2/dy),in[1]) + \f]
1837 *
1838 * \f[ ((d\phi/d\xi_0\, d\xi_0/dz +
1839 * d\phi/d\xi_1\, d\xi_1/dz +
1840 * d\phi/d\xi_2\, d\xi_2/dz),in[2]) \, \f]
1841 *
1842 * where we note that
1843 *
1844 * \f[ d\phi/d\xi_0 =
1845 * d\phi/d\eta_0\, d\eta_0/d\xi_0 =
1846 * d\phi/d\eta_0\, 2/(1-\eta_2). \f]
1847 *
1848 * \f[ d\phi/d\xi_1 =
1849 * d\phi/d\eta_1\, d\eta_1/d\xi_1 =
1850 * d\phi/d\eta_1\, 2/(1-\eta_2) \f]
1851 *
1852 * \f[ d\phi/d\xi_2 =
1853 * d\phi/d\eta_0\, d\eta_0/d\xi_2 +
1854 * d\phi/d\eta_1\, d\eta_1/d\xi_2 +
1855 * d\phi/d\eta_2\, d\eta_2/d\xi_2 =
1856 * d\phi/d\eta_0 (1+\eta_0)/(1-\eta_2) +
1857 * d\phi/d\eta_1 (1+\eta_1)/(1-\eta_2) + d\phi/d\eta_2 \f]
1858 *
1859 * and so the full inner products are
1860 *
1861 * \f[ (d\phi/dx,in[0]) + (d\phi/dy,in[1]) + (d\phi/dz,in[2]) = \f]
1862 *
1863 * \f[ (d\phi/d\eta_0, ((2/(1-\eta_2) (d\xi_0/dx in[0] +
1864 * d\xi_0/dy in[1] +
1865 * (1-\eta_0)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1]
1866 * + d\xi_2/dz in[2] )) + \f]
1867 * \f[ (d\phi/d\eta_1, ((2/(1-\eta_2) (d\xi_1/dx in[0] +
1868 * d\xi_0/dy in[1] + d\xi_0/dz in[2]) +
1869 * (1-\eta_1)/(1-\eta_2) (d\xi_2/dx in[0] + d\xi_2/dy in[1] +
1870 * d\xi_2/dz in[2] )) \f]
1871 *
1872 * \f[ (d\phi/d\eta_2, (d\xi_2/dx in[0] + d\xi_2/dy in[1] +
1873 * d\xi_2/dz in[2])) \f]
1874 */
1875 void operator()(const Array<OneD, const NekDouble> &entry0,
1876 Array<OneD, NekDouble> &entry1,
1877 Array<OneD, NekDouble> &entry2,
1878 Array<OneD, NekDouble> &entry3,
1879 Array<OneD, NekDouble> &wsp) final
1880 {
1881 unsigned int nPhys = m_stdExp->GetTotPoints();
1882 unsigned int ntot = m_numElmt * nPhys;
1883 unsigned int nmodes = m_stdExp->GetNcoeffs();
1884 unsigned int nmax = max(ntot, m_numElmt * nmodes);
1885 Array<OneD, Array<OneD, const NekDouble>> in(3);
1886 Array<OneD, NekDouble> output, wsp1;
1887 Array<OneD, Array<OneD, NekDouble>> tmp(3);
1888
1889 in[0] = entry0;
1890 in[1] = entry1;
1891 in[2] = entry2;
1892
1893 output = entry3;
1894
1895 for (int i = 0; i < 3; ++i)
1896 {
1897 tmp[i] = wsp + i * nmax;
1898 }
1899
1900 if (m_isDeformed)
1901 {
1902 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1903 for (int i = 0; i < 3; ++i)
1904 {
1905 Vmath::Vmul(ntot, m_derivFac[i], 1, in[0], 1, tmp[i], 1);
1906 for (int j = 1; j < 3; ++j)
1907 {
1908 Vmath::Vvtvp(ntot, m_derivFac[i + 3 * j], 1, in[j], 1,
1909 tmp[i], 1, tmp[i], 1);
1910 }
1911 }
1912 }
1913 else
1914 {
1915 Array<OneD, NekDouble> t;
1916 for (int e = 0; e < m_numElmt; ++e)
1917 {
1918 // calculate dx/dxi in[0] + dy/dxi in[1] + dz/dxi in[2]
1919 for (int i = 0; i < 3; ++i)
1920 {
1921 Vmath::Smul(m_nqe, m_derivFac[i][e], in[0] + e * m_nqe, 1,
1922 t = tmp[i] + e * m_nqe, 1);
1923 for (int j = 1; j < 3; ++j)
1924 {
1925 Vmath::Svtvp(m_nqe, m_derivFac[i + 3 * j][e],
1926 in[j] + e * m_nqe, 1, tmp[i] + e * m_nqe,
1927 1, t = tmp[i] + e * m_nqe, 1);
1928 }
1929 }
1930 }
1931 }
1932
1933 wsp1 = wsp + 3 * nmax;
1934
1935 // Sort into eta factors
1936 for (int i = 0; i < m_numElmt; ++i)
1937 {
1938 // scale tmp[0] by fac0
1939 Vmath::Vmul(nPhys, &m_fac0[0], 1, tmp[0].get() + i * nPhys, 1,
1940 tmp[0].get() + i * nPhys, 1);
1941
1942 // scale tmp[2] by fac1 and add to tmp0
1943 Vmath::Vvtvp(nPhys, &m_fac1[0], 1, tmp[2].get() + i * nPhys, 1,
1944 tmp[0].get() + i * nPhys, 1, tmp[0].get() + i * nPhys,
1945 1);
1946
1947 // scale tmp[1] by fac0
1948 Vmath::Vmul(nPhys, &m_fac0[0], 1, tmp[1].get() + i * nPhys, 1,
1949 tmp[1].get() + i * nPhys, 1);
1950
1951 // scale tmp[2] by fac2 and add to tmp1
1952 Vmath::Vvtvp(nPhys, &m_fac2[0], 1, tmp[2].get() + i * nPhys, 1,
1953 tmp[1].get() + i * nPhys, 1, tmp[1].get() + i * nPhys,
1954 1);
1955 }
1956
1957 // calculate Iproduct WRT Std Deriv
1958 PyrIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1959 m_nmodes0, m_nmodes1, m_nmodes2, m_derbase0, m_base1,
1960 m_base2, m_jacWStdW, tmp[0], output, wsp1);
1961
1962 PyrIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1963 m_nmodes0, m_nmodes1, m_nmodes2, m_base0, m_derbase1,
1964 m_base2, m_jacWStdW, tmp[1], tmp[0], wsp1);
1965 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1966
1967 PyrIProduct(m_sortTopVertex, m_numElmt, m_nquad0, m_nquad1, m_nquad2,
1968 m_nmodes0, m_nmodes1, m_nmodes2, m_base0, m_base1,
1969 m_derbase2, m_jacWStdW, tmp[2], tmp[0], wsp1);
1970 Vmath::Vadd(m_numElmt * nmodes, tmp[0], 1, output, 1, output, 1);
1971 }
1972
1973 void operator()([[maybe_unused]] int dir,
1974 [[maybe_unused]] const Array<OneD, const NekDouble> &input,
1975 [[maybe_unused]] Array<OneD, NekDouble> &output,
1976 [[maybe_unused]] Array<OneD, NekDouble> &wsp) final
1977 {
1978 NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
1979 }
1980
1981protected:
1982 const int m_nquad0;
1983 const int m_nquad1;
1984 const int m_nquad2;
1985 const int m_nmodes0;
1986 const int m_nmodes1;
1987 const int m_nmodes2;
1988 Array<OneD, const NekDouble> m_jacWStdW;
1989 Array<OneD, const NekDouble> m_base0;
1990 Array<OneD, const NekDouble> m_base1;
1991 Array<OneD, const NekDouble> m_base2;
1992 Array<OneD, const NekDouble> m_derbase0;
1993 Array<OneD, const NekDouble> m_derbase1;
1994 Array<OneD, const NekDouble> m_derbase2;
1995 Array<TwoD, const NekDouble> m_derivFac;
1996 Array<OneD, NekDouble> m_fac0;
1997 Array<OneD, NekDouble> m_fac1;
1998 Array<OneD, NekDouble> m_fac2;
1999 bool m_sortTopVertex;
2000
2001private:
2002 IProductWRTDerivBase_SumFac_Pyr(
2003 vector<StdRegions::StdExpansionSharedPtr> pCollExp,
2004 CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
2005 : Operator(pCollExp, pGeomData, factors), IProductWRTDerivBase_Helper(),
2006 m_nquad0(m_stdExp->GetNumPoints(0)),
2007 m_nquad1(m_stdExp->GetNumPoints(1)),
2008 m_nquad2(m_stdExp->GetNumPoints(2)),
2009 m_nmodes0(m_stdExp->GetBasisNumModes(0)),
2010 m_nmodes1(m_stdExp->GetBasisNumModes(1)),
2011 m_nmodes2(m_stdExp->GetBasisNumModes(2)),
2012 m_base0(m_stdExp->GetBasis(0)->GetBdata()),
2013 m_base1(m_stdExp->GetBasis(1)->GetBdata()),
2014 m_base2(m_stdExp->GetBasis(2)->GetBdata()),
2015 m_derbase0(m_stdExp->GetBasis(0)->GetDbdata()),
2016 m_derbase1(m_stdExp->GetBasis(1)->GetDbdata()),
2017 m_derbase2(m_stdExp->GetBasis(2)->GetDbdata())
2018
2019 {
2020 m_jacWStdW = pGeomData->GetJacWithStdWeights(pCollExp);
2021 m_wspSize = 6 * m_numElmt *
2022 (max(m_nquad0 * m_nquad1 * m_nquad2,
2023 m_nmodes0 * m_nmodes1 * m_nmodes2));
2024 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
2025
2026 if (m_stdExp->GetBasis(0)->GetBasisType() == LibUtilities::eModified_A)
2027 {
2028 m_sortTopVertex = true;
2029 }
2030 else
2031 {
2032 m_sortTopVertex = false;
2033 }
2034
2035 const Array<OneD, const NekDouble> &z0 = m_stdExp->GetBasis(0)->GetZ();
2036 const Array<OneD, const NekDouble> &z1 = m_stdExp->GetBasis(1)->GetZ();
2037 const Array<OneD, const NekDouble> &z2 = m_stdExp->GetBasis(2)->GetZ();
2038
2039 m_fac0 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
2040 m_fac1 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
2041 m_fac2 = Array<OneD, NekDouble>(m_nquad0 * m_nquad1 * m_nquad2);
2042
2043 for (int i = 0; i < m_nquad0; ++i)
2044 {
2045 for (int j = 0; j < m_nquad1; ++j)
2046 {
2047 for (int k = 0; k < m_nquad2; ++k)
2048 {
2049 // set up geometric factor: 2/(1-z2)
2050 m_fac0[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
2051 2.0 / (1 - z2[k]);
2052 // set up geometric factor: (1+z0)/(1-z2)
2053 m_fac1[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
2054 (1 + z0[i]) / (1 - z2[k]);
2055 // set up geometric factor: (1+z1)/(1-z2)
2056 m_fac2[i + j * m_nquad0 + k * m_nquad0 * m_nquad1] =
2057 (1 + z1[j]) / (1 - z2[k]);
2058 }
2059 }
2060 }
2061 }
2062};
2063
2064/// Factory initialisation for the IProductWRTDerivBase_SumFac_Pyr operator
2065OperatorKey IProductWRTDerivBase_SumFac_Pyr::m_type =
2066 GetOperatorFactory().RegisterCreatorFunction(
2067 OperatorKey(ePyramid, eIProductWRTDerivBase, eSumFac, false),
2068 IProductWRTDerivBase_SumFac_Pyr::create,
2069 "IProductWRTDerivBase_SumFac_Pyr");
2070
2071} // namespace Nektar::Collections
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
NEKERROR
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
Definition: ErrorUtil.hpp:202
OPERATOR_CREATE
#define OPERATOR_CREATE(cname)
Definition: Operator.h:43
Nektar::Collections::IProductWRTDerivBase_Helper
Inner product deriv base help class to calculate the size of the collection that is given as an input...
Definition: IProductWRTDerivBase.cpp:61
Nektar::Collections::IProductWRTDerivBase_Helper::IProductWRTDerivBase_Helper
IProductWRTDerivBase_Helper()
Definition: IProductWRTDerivBase.cpp:63
Nektar::Collections::IProductWRTDerivBase_IterPerExp
Inner product WRT deriv base operator using element-wise operation.
Definition: IProductWRTDerivBase.cpp:410
Nektar::Collections::IProductWRTDerivBase_IterPerExp::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:523
Nektar::Collections::IProductWRTDerivBase_IterPerExp::IProductWRTDerivBase_IterPerExp
IProductWRTDerivBase_IterPerExp(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:529
Nektar::Collections::IProductWRTDerivBase_IterPerExp::m_dim
int m_dim
Definition: IProductWRTDerivBase.cpp:525
Nektar::Collections::IProductWRTDerivBase_IterPerExp::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:526
Nektar::Collections::IProductWRTDerivBase_IterPerExp::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:514
Nektar::Collections::IProductWRTDerivBase_IterPerExp::m_jac
Array< OneD, const NekDouble > m_jac
Definition: IProductWRTDerivBase.cpp:524
Nektar::Collections::IProductWRTDerivBase_IterPerExp::~IProductWRTDerivBase_IterPerExp
~IProductWRTDerivBase_IterPerExp() final=default
Nektar::Collections::IProductWRTDerivBase_MatrixFree
Inner product operator using operator using matrix free operators.
Definition: IProductWRTDerivBase.cpp:276
Nektar::Collections::IProductWRTDerivBase_MatrixFree::IProductWRTDerivBase_MatrixFree
IProductWRTDerivBase_MatrixFree(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:329
Nektar::Collections::IProductWRTDerivBase_MatrixFree::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:327
Nektar::Collections::IProductWRTDerivBase_MatrixFree::m_oper
std::shared_ptr< MatrixFree::IProductWRTDerivBase > m_oper
Definition: IProductWRTDerivBase.cpp:326
Nektar::Collections::IProductWRTDerivBase_MatrixFree::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:317
Nektar::Collections::IProductWRTDerivBase_MatrixFree::~IProductWRTDerivBase_MatrixFree
~IProductWRTDerivBase_MatrixFree() final=default
Nektar::Collections::IProductWRTDerivBase_NoCollection
Inner product WRT deriv base operator using LocalRegions implementation.
Definition: IProductWRTDerivBase.cpp:594
Nektar::Collections::IProductWRTDerivBase_NoCollection::m_expList
vector< StdRegions::StdExpansionSharedPtr > m_expList
Definition: IProductWRTDerivBase.cpp:647
Nektar::Collections::IProductWRTDerivBase_NoCollection::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:636
Nektar::Collections::IProductWRTDerivBase_NoCollection::IProductWRTDerivBase_NoCollection
IProductWRTDerivBase_NoCollection(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:650
Nektar::Collections::IProductWRTDerivBase_NoCollection::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:646
Nektar::Collections::IProductWRTDerivBase_NoCollection::~IProductWRTDerivBase_NoCollection
~IProductWRTDerivBase_NoCollection() final=default
Nektar::Collections::IProductWRTDerivBase_NoCollection::m_dim
int m_dim
Definition: IProductWRTDerivBase.cpp:645
Nektar::Collections::IProductWRTDerivBase_StdMat
Inner product WRT deriv base operator using standard matrix approach.
Definition: IProductWRTDerivBase.cpp:77
Nektar::Collections::IProductWRTDerivBase_StdMat::m_jac
Array< OneD, const NekDouble > m_jac
Definition: IProductWRTDerivBase.cpp:189
Nektar::Collections::IProductWRTDerivBase_StdMat::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:188
Nektar::Collections::IProductWRTDerivBase_StdMat::m_iProdWRTStdDBase
Array< OneD, DNekMatSharedPtr > m_iProdWRTStdDBase
Definition: IProductWRTDerivBase.cpp:187
Nektar::Collections::IProductWRTDerivBase_StdMat::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:178
Nektar::Collections::IProductWRTDerivBase_StdMat::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:191
Nektar::Collections::IProductWRTDerivBase_StdMat::~IProductWRTDerivBase_StdMat
~IProductWRTDerivBase_StdMat() final=default
Nektar::Collections::IProductWRTDerivBase_StdMat::m_dim
int m_dim
Definition: IProductWRTDerivBase.cpp:190
Nektar::Collections::IProductWRTDerivBase_StdMat::IProductWRTDerivBase_StdMat
IProductWRTDerivBase_StdMat(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:194
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex
Inner product WRT deriv base operator using sum-factorisation (Hex)
Definition: IProductWRTDerivBase.cpp:1161
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_colldir0
const bool m_colldir0
Definition: IProductWRTDerivBase.cpp:1261
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1269
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nquad2
const int m_nquad2
Definition: IProductWRTDerivBase.cpp:1257
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1246
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_colldir1
const bool m_colldir1
Definition: IProductWRTDerivBase.cpp:1262
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1264
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_base2
Array< OneD, const NekDouble > m_base2
Definition: IProductWRTDerivBase.cpp:1267
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::IProductWRTDerivBase_SumFac_Hex
IProductWRTDerivBase_SumFac_Hex(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1274
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1256
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1271
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_derbase2
Array< OneD, const NekDouble > m_derbase2
Definition: IProductWRTDerivBase.cpp:1270
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_colldir2
const bool m_colldir2
Definition: IProductWRTDerivBase.cpp:1263
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1258
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1268
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1255
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1266
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1265
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::~IProductWRTDerivBase_SumFac_Hex
~IProductWRTDerivBase_SumFac_Hex() final=default
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nmodes2
const int m_nmodes2
Definition: IProductWRTDerivBase.cpp:1260
Nektar::Collections::IProductWRTDerivBase_SumFac_Hex::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1259
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism
Inner product WRT deriv base operator using sum-factorisation (Prism)
Definition: IProductWRTDerivBase.cpp:1578
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1732
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_base2
Array< OneD, const NekDouble > m_base2
Definition: IProductWRTDerivBase.cpp:1737
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_derbase2
Array< OneD, const NekDouble > m_derbase2
Definition: IProductWRTDerivBase.cpp:1740
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1738
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1729
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1734
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1731
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::IProductWRTDerivBase_SumFac_Prism
IProductWRTDerivBase_SumFac_Prism(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1747
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1736
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nquad2
const int m_nquad2
Definition: IProductWRTDerivBase.cpp:1730
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1719
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: IProductWRTDerivBase.cpp:1742
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1739
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nmodes2
const int m_nmodes2
Definition: IProductWRTDerivBase.cpp:1733
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: IProductWRTDerivBase.cpp:1743
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1735
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1741
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_sortTopVertex
bool m_sortTopVertex
Definition: IProductWRTDerivBase.cpp:1744
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::~IProductWRTDerivBase_SumFac_Prism
~IProductWRTDerivBase_SumFac_Prism() final=default
Nektar::Collections::IProductWRTDerivBase_SumFac_Prism::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1728
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr
Inner product WRT deriv base operator using sum-factorisation (Pyr)
Definition: IProductWRTDerivBase.cpp:1816
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nquad2
const int m_nquad2
Definition: IProductWRTDerivBase.cpp:1984
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1973
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1995
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_base2
Array< OneD, const NekDouble > m_base2
Definition: IProductWRTDerivBase.cpp:1991
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1990
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1983
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_sortTopVertex
bool m_sortTopVertex
Definition: IProductWRTDerivBase.cpp:1999
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: IProductWRTDerivBase.cpp:1996
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1993
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::~IProductWRTDerivBase_SumFac_Pyr
~IProductWRTDerivBase_SumFac_Pyr() final=default
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nmodes2
const int m_nmodes2
Definition: IProductWRTDerivBase.cpp:1987
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: IProductWRTDerivBase.cpp:1997
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1988
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1992
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1985
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1989
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1986
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::IProductWRTDerivBase_SumFac_Pyr
IProductWRTDerivBase_SumFac_Pyr(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:2002
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_fac2
Array< OneD, NekDouble > m_fac2
Definition: IProductWRTDerivBase.cpp:1998
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1982
Nektar::Collections::IProductWRTDerivBase_SumFac_Pyr::m_derbase2
Array< OneD, const NekDouble > m_derbase2
Definition: IProductWRTDerivBase.cpp:1994
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad
Inner product WRT deriv base operator using sum-factorisation (Quad)
Definition: IProductWRTDerivBase.cpp:814
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:900
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:889
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:901
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::~IProductWRTDerivBase_SumFac_Quad
~IProductWRTDerivBase_SumFac_Quad() final=default
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:905
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:907
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:904
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:910
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:909
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_colldir1
const bool m_colldir1
Definition: IProductWRTDerivBase.cpp:903
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_colldir0
const bool m_colldir0
Definition: IProductWRTDerivBase.cpp:902
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:898
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:906
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::IProductWRTDerivBase_SumFac_Quad
IProductWRTDerivBase_SumFac_Quad(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:913
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:899
Nektar::Collections::IProductWRTDerivBase_SumFac_Quad::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:908
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg
Inner product WRT deriv base operator using sum-factorisation (Segment)
Definition: IProductWRTDerivBase.cpp:711
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:770
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:781
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:784
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:782
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:780
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:779
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::IProductWRTDerivBase_SumFac_Seg
IProductWRTDerivBase_SumFac_Seg(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:787
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::~IProductWRTDerivBase_SumFac_Seg
~IProductWRTDerivBase_SumFac_Seg() final=default
Nektar::Collections::IProductWRTDerivBase_SumFac_Seg::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:783
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet
Inner product WRT deriv base operator using sum-factorisation (Tet)
Definition: IProductWRTDerivBase.cpp:1315
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1474
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1483
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_fac2
Array< OneD, NekDouble > m_fac2
Definition: IProductWRTDerivBase.cpp:1499
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: IProductWRTDerivBase.cpp:1498
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_base2
Array< OneD, const NekDouble > m_base2
Definition: IProductWRTDerivBase.cpp:1492
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::IProductWRTDerivBase_SumFac_Tet
IProductWRTDerivBase_SumFac_Tet(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1504
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1486
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_sortTopEdge
bool m_sortTopEdge
Definition: IProductWRTDerivBase.cpp:1501
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1490
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nquad2
const int m_nquad2
Definition: IProductWRTDerivBase.cpp:1485
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_fac3
Array< OneD, NekDouble > m_fac3
Definition: IProductWRTDerivBase.cpp:1500
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: IProductWRTDerivBase.cpp:1497
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1489
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1496
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1494
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nmodes2
const int m_nmodes2
Definition: IProductWRTDerivBase.cpp:1488
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1491
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1493
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1487
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1484
Nektar::Collections::IProductWRTDerivBase_SumFac_Tet::m_derbase2
Array< OneD, const NekDouble > m_derbase2
Definition: IProductWRTDerivBase.cpp:1495
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri
Inner product WRT deriv base operator using sum-factorisation (Tri)
Definition: IProductWRTDerivBase.cpp:948
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_fac0
Array< OneD, NekDouble > m_fac0
Definition: IProductWRTDerivBase.cpp:1089
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_nquad1
const int m_nquad1
Definition: IProductWRTDerivBase.cpp:1077
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_base0
Array< OneD, const NekDouble > m_base0
Definition: IProductWRTDerivBase.cpp:1085
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_base1
Array< OneD, const NekDouble > m_base1
Definition: IProductWRTDerivBase.cpp:1086
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_coordim
int m_coordim
Definition: IProductWRTDerivBase.cpp:1082
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_colldir0
const bool m_colldir0
Definition: IProductWRTDerivBase.cpp:1080
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_jacWStdW
Array< OneD, const NekDouble > m_jacWStdW
Definition: IProductWRTDerivBase.cpp:1084
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_nmodes1
const int m_nmodes1
Definition: IProductWRTDerivBase.cpp:1079
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_fac1
Array< OneD, NekDouble > m_fac1
Definition: IProductWRTDerivBase.cpp:1090
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::operator()
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
Definition: IProductWRTDerivBase.cpp:1067
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_colldir1
const bool m_colldir1
Definition: IProductWRTDerivBase.cpp:1081
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_derivFac
Array< TwoD, const NekDouble > m_derivFac
Definition: IProductWRTDerivBase.cpp:1083
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_nmodes0
const int m_nmodes0
Definition: IProductWRTDerivBase.cpp:1078
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_derbase0
Array< OneD, const NekDouble > m_derbase0
Definition: IProductWRTDerivBase.cpp:1087
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_nquad0
const int m_nquad0
Definition: IProductWRTDerivBase.cpp:1076
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::IProductWRTDerivBase_SumFac_Tri
IProductWRTDerivBase_SumFac_Tri(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Definition: IProductWRTDerivBase.cpp:1094
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::~IProductWRTDerivBase_SumFac_Tri
~IProductWRTDerivBase_SumFac_Tri() final=default
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_sortTopVertex
bool m_sortTopVertex
Definition: IProductWRTDerivBase.cpp:1091
Nektar::Collections::IProductWRTDerivBase_SumFac_Tri::m_derbase1
Array< OneD, const NekDouble > m_derbase1
Definition: IProductWRTDerivBase.cpp:1088
Nektar::Collections::MatrixFreeBase::m_nElmtPad
unsigned int m_nElmtPad
size after padding
Definition: MatrixFreeBase.h:70
Nektar::Collections::Operator
Base class for operators on a collection of elements.
Definition: Operator.h:138
Nektar::Collections::Operator::m_stdExp
StdRegions::StdExpansionSharedPtr m_stdExp
Definition: Operator.h:217
Nektar::Collections::Operator::m_numElmt
unsigned int m_numElmt
number of elements that the operator is applied on
Definition: Operator.h:219
Nektar::Collections::Operator::m_outputSize
unsigned int m_outputSize
number of modes or quadrature points that are taken as output from an operator
Definition: Operator.h:227
Nektar::Collections::Operator::m_inputSize
unsigned int m_inputSize
number of modes or quadrature points that are passed as input to an operator
Definition: Operator.h:224
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: BasicUtils/NekFactory.hpp:197
Nektar::LibUtilities::NekFactory::CreateInstance
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Definition: BasicUtils/NekFactory.hpp:143
Nektar::MemoryManager::AllocateSharedPtr
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Definition: NekMemoryManager.hpp:166
Blas::Dgemm
static void Dgemm(const char &transa, const char &transb, const int &m, const int &n, const int &k, const double &alpha, const double *a, const int &lda, const double *b, const int &ldb, const double &beta, double *c, const int &ldc)
BLAS level 3: Matrix-matrix multiply C = A x B where op(A)[m x k], op(B)[k x n], C[m x n] DGEMM perfo...
Definition: Blas.hpp:383
Nektar::Collections::TetIProduct
void TetIProduct(bool sortTopEdge, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:496
Nektar::Collections::PyrIProduct
void PyrIProduct(bool sortTopVert, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:405
Nektar::Collections::TriIProduct
void TriIProduct(bool sortTopVertex, int numElmt, int nquad0, int nquad1, int nmodes0, int nmodes1, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:128
Nektar::Collections::OperatorKey
std::tuple< LibUtilities::ShapeType, OperatorType, ImplementationType, ExpansionIsNodal > OperatorKey
Key for describing an Operator.
Definition: Operator.h:120
Nektar::Collections::QuadIProduct
void QuadIProduct(bool colldir0, bool colldir1, int numElmt, int nquad0, int nquad1, int nmodes0, int nmodes1, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:46
Nektar::Collections::CoalescedGeomDataSharedPtr
std::shared_ptr< CoalescedGeomData > CoalescedGeomDataSharedPtr
Definition: CoalescedGeomData.h:88
Nektar::Collections::GetOperatorFactory
OperatorFactory & GetOperatorFactory()
Returns the singleton Operator factory object.
Definition: Operator.cpp:44
Nektar::Collections::PrismIProduct
void PrismIProduct(bool sortTopVert, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:338
Nektar::Collections::HexIProduct
void HexIProduct(bool colldir0, bool colldir1, bool colldir2, int numElmt, int nquad0, int nquad1, int nquad2, int nmodes0, int nmodes1, int nmodes2, const Array< OneD, const NekDouble > &base0, const Array< OneD, const NekDouble > &base1, const Array< OneD, const NekDouble > &base2, const Array< OneD, const NekDouble > &jac, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp)
Definition: IProduct.cpp:170
Nektar::LibUtilities::eModified_A
@ eModified_A
Principle Modified Functions .
Definition: BasisType.h:48
Nektar::StdRegions::FactorMap
ConstFactorMap FactorMap
Definition: StdRegions.hpp:434
Nektar::VarcoeffHashingTest::factors
StdRegions::ConstFactorMap factors
Definition: TestVarcoeffHashing.cpp:52
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Vmath::Vmul
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
Vmath::Svtvp
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Svtvp (scalar times vector plus vector): z = alpha*x + y.
Definition: Vmath.hpp:396
Vmath::Vvtvp
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
Vmath::Vadd
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
Vmath::Smul
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
Vmath::Zero
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
Vmath::Vcopy
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
std
STL namespace.
Generated on Wed Oct 30 2024 16:07:28 for Nektar++ by doxygen 1.9.4