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