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