39#include <MatrixFreeOps/Operator.hpp>
90 unsigned int nmodes =
m_expList[0]->GetNcoeffs();
98 m_expList[n]->GeneralMatrixOp(entry0 + n * nmodes,
99 tmp = entry1 + n * nmodes, mkey);
112 [[maybe_unused]]
int coll_phys_offset)
override
130 m_dim = pCollExp[0]->GetNumBases();
141 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_Seg"),
144 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_Tri"),
147 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_NodalTri"),
150 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_Quad"),
153 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_Tet"),
156 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_NodalTet"),
159 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_Pyr"),
162 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_Prism"),
165 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_NodalPrism"),
168 Helmholtz_NoCollection::create,
"Helmholtz_NoCollection_Hex")};
186 const int nCoeffs =
m_stdExp->GetNcoeffs();
187 const int nPhys =
m_stdExp->GetTotPoints();
190 "input array size is insufficient");
192 "output array size is insufficient");
201 dtmp[i - 1] = wsp + i * nPhys;
202 tmp[i - 1] = wsp + (i +
m_coordim) * nPhys;
207 m_stdExp->BwdTrans(input + i * nCoeffs, tmpphys);
210 m_stdExp->PhysDeriv(tmpphys, dtmp[0], dtmp[1], dtmp[2]);
223 m_stdExp->IProductWRTBase(tmpphys, t1 = output + i * nCoeffs);
225 t1 = output + i * nCoeffs, 1);
234 &dtmp[0][0], 1, &tmp[j][0], 1);
236 for (
int k = 1; k <
m_dim; ++k)
241 &dtmp[k][0], 1, &tmp[j][0], 1, &tmp[j][0], 1);
257 &tmpphys[0], 1, &tmpphys[0], 1);
260 for (
int j = 0; j <
m_dim; ++j)
263 1, &tmpphys[0], 1, &dtmp[j][0], 1);
274 &tmpphys[0], 1, &tmpphys[0], 1);
277 for (
int j = 0; j <
m_dim; ++j)
282 1, &tmpphys[0], 1, &dtmp[j][0], 1,
291 for (
int j = 0; j <
m_dim; ++j)
294 1, &tmp[0][0], 1, &dtmp[j][0], 1);
301 1, &tmp[k][0], 1, &dtmp[j][0], 1,
308 for (
int j = 0; j <
m_dim; ++j)
314 m_stdExp->IProductWRTDerivBase(j, dtmp[j], tmp[0]);
315 Vmath::Vadd(nCoeffs, tmp[0], 1, output + i * nCoeffs, 1,
316 t1 = output + i * nCoeffs, 1);
327 for (
int k = 1; k <
m_dim; ++k)
330 &dtmp[k][0], 1, &tmp[j][0], 1, &tmp[j][0],
347 &tmpphys[0], 1, &tmpphys[0], 1);
350 for (
int j = 0; j <
m_dim; ++j)
364 &tmpphys[0], 1, &tmpphys[0], 1);
367 for (
int j = 0; j <
m_dim; ++j)
370 &tmpphys[0], 1, &dtmp[j][0], 1,
379 for (
int j = 0; j <
m_dim; ++j)
387 &tmp[k][0], 1, &dtmp[j][0], 1,
394 for (
int j = 0; j <
m_dim; ++j)
399 m_stdExp->IProductWRTDerivBase(j, dtmp[j], tmp[0]);
400 Vmath::Vadd(nCoeffs, tmp[0], 1, output + i * nCoeffs, 1,
401 t1 = output + i * nCoeffs, 1);
422 [[maybe_unused]]
int coll_phys_offset)
override
436 "Constant factor not defined: " +
502 m_dim = pCollExp[0]->GetShapeDimension();
504 int nqtot =
m_stdExp->GetTotPoints();
506 m_derivFac = pGeomData->GetDerivFactors(pCollExp);
507 m_jac = pGeomData->GetJac(pCollExp);
521 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_Seg"),
524 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_Tri"),
527 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_NodalTri"),
530 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_Quad"),
533 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_Tet"),
536 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_NodalTet"),
539 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_Pyr"),
542 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_Prism"),
545 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_NodalPrism"),
548 Helmholtz_IterPerExp::create,
"Helmholtz_IterPerExp_Hex")};
577 (*m_oper)(input, output0);
593 [[maybe_unused]]
int coll_phys_offset)
override
606 "Constant factor not defined: " +
609 m_oper->SetLambda(x->second);
612 bool isConstVarDiff =
false;
614 diff[0] = diff[2] = diff[5] = 1.0;
617 if (xd00 !=
factors.end() && xd00->second != 1.0)
619 isConstVarDiff =
true;
620 diff[0] = xd00->second;
624 if (xd01 !=
factors.end() && xd01->second != 0.0)
626 isConstVarDiff =
true;
627 diff[1] = xd01->second;
631 if (xd11 !=
factors.end() && xd11->second != 1.0)
633 isConstVarDiff =
true;
634 diff[2] = xd11->second;
638 if (xd02 !=
factors.end() && xd02->second != 0.0)
640 isConstVarDiff =
true;
641 diff[3] = xd02->second;
645 if (xd12 !=
factors.end() && xd12->second != 0.0)
647 isConstVarDiff =
true;
648 diff[4] = xd12->second;
652 if (xd22 !=
factors.end() && xd22->second != 1.0)
654 isConstVarDiff =
true;
655 diff[5] = xd22->second;
660 m_oper->SetConstVarDiffusion(diff);
665 if (k !=
factors.end() && k->second != 0.0)
667 m_oper->SetVarDiffusion(diff);
672 std::shared_ptr<MatrixFree::Helmholtz>
m_oper;
681 pCollExp[0]->GetStdExp()->GetNcoeffs(),
688 const auto dim = pCollExp[0]->GetStdExp()->GetShapeDimension();
691 std::vector<LibUtilities::BasisSharedPtr> basis(dim);
692 for (
auto i = 0; i < dim; ++i)
694 basis[i] = pCollExp[0]->GetBasis(i);
698 auto shapeType = pCollExp[0]->GetStdExp()->DetShapeType();
701 std::string op_string =
"Helmholtz";
702 op_string += MatrixFree::GetOpstring(shapeType,
m_isDeformed);
707 oper->SetJac(pGeomData->GetJacInterLeave(pCollExp,
m_nElmtPad));
710 oper->SetDF(pGeomData->GetDerivFactorsInterLeave(pCollExp,
m_nElmtPad));
712 m_oper = std::dynamic_pointer_cast<MatrixFree::Helmholtz>(oper);
725 Helmholtz_MatrixFree::create,
"Helmholtz_MatrixFree_Quad"),
728 Helmholtz_MatrixFree::create,
"Helmholtz_MatrixFree_Tri"),
731 Helmholtz_MatrixFree::create,
"Helmholtz_MatrixFree_Hex"),
734 Helmholtz_MatrixFree::create,
"Helmholtz_MatrixFree_Prism"),
737 Helmholtz_MatrixFree::create,
"Helmholtz_MatrixFree_Pyr"),
740 Helmholtz_MatrixFree::create,
"Helmholtz_MatrixFree_Tet"),
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define OPERATOR_CREATE(cname)
Helmholtz help class to calculate the size of the collection that is given as an input and as an outp...
Helmholtz operator using LocalRegions implementation.
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
~Helmholtz_IterPerExp() final=default
StdRegions::FactorMap m_factors
Array< TwoD, const NekDouble > m_derivFac
const StdRegions::ConstFactorType m_factorCoeffDef[3][3]
Helmholtz_IterPerExp(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset) override
Check the validity of supplied constant factors.
Array< OneD, Array< OneD, NekDouble > > m_diff
Array< OneD, const NekDouble > m_jac
Helmholtz operator using matrix free operators.
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
~Helmholtz_MatrixFree() final=default
void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset) override
Check the validity of the supplied factor map.
std::shared_ptr< MatrixFree::Helmholtz > m_oper
StdRegions::FactorMap m_factors
Helmholtz_MatrixFree(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
Helmholtz operator using LocalRegions implementation.
Helmholtz_NoCollection(vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)
void operator()(int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final
StdRegions::FactorMap m_factors
~Helmholtz_NoCollection() final=default
vector< StdRegions::StdExpansionSharedPtr > m_expList
void CheckFactors(StdRegions::FactorMap factors, int coll_phys_offset) override
Check the validity of the supplied factor map.
bool m_isPadded
flag for padding
unsigned int m_nElmtPad
size after padding
Array< OneD, NekDouble > m_output
Array< OneD, NekDouble > m_input
padded input/output vectors
Base class for operators on a collection of elements.
StdRegions::StdExpansionSharedPtr m_stdExp
unsigned int m_numElmt
number of elements that the operator is applied on
unsigned int m_outputSize
number of modes or quadrature points that are taken as output from an operator
unsigned int m_inputSize
number of modes or quadrature points that are passed as input to an operator
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
std::tuple< LibUtilities::ShapeType, OperatorType, ImplementationType, ExpansionIsNodal > OperatorKey
Key for describing an Operator.
std::shared_ptr< CoalescedGeomData > CoalescedGeomDataSharedPtr
OperatorFactory & GetOperatorFactory()
Returns the singleton Operator factory object.
const NekPoint< NekDouble > origin
static FactorMap NullFactorMap
const char *const ConstFactorTypeMap[]
std::vector< double > d(NPUPPER *NPUPPER)
StdRegions::ConstFactorMap factors
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.
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.
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
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)