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