Helmholtz operator using LocalRegions implementation. More...

Inheritance diagram for Nektar::Collections::Helmholtz_IterPerExp:

Public Member Functions
	~Helmholtz_IterPerExp () final=default

void	operator() (const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp) final
	Perform operation. More...

void	operator() (int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp) final

void	UpdateFactors (StdRegions::FactorMap factors) override
	Check the validity of supplied constant factors. More...

void	UpdateVarcoeffs (StdRegions::VarCoeffMap &varcoeffs) override
	Check the validity of supplied variable coefficients. Note that the m_varcoeffs are not implemented for the IterPerExp operator. There exists a specialised implementation for variable diffusion in the above function UpdateFactors. More...

Public Member Functions inherited from Nektar::Collections::Operator
	Operator (std::vector< StdRegions::StdExpansionSharedPtr > pCollExp, std::shared_ptr< CoalescedGeomData > GeomData, StdRegions::FactorMap factors)
	Constructor. More...

virtual	~Operator ()=default

virtual COLLECTIONS_EXPORT void	operator() (const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output0, Array< OneD, NekDouble > &output1, Array< OneD, NekDouble > &output2, Array< OneD, NekDouble > &wsp=NullNekDouble1DArray)=0
	Perform operation. More...

virtual COLLECTIONS_EXPORT void	operator() (int dir, const Array< OneD, const NekDouble > &input, Array< OneD, NekDouble > &output, Array< OneD, NekDouble > &wsp=NullNekDouble1DArray)=0

virtual COLLECTIONS_EXPORT void	UpdateFactors (StdRegions::FactorMap factors)
	Update the supplied factor map. More...

virtual COLLECTIONS_EXPORT void	UpdateVarcoeffs (StdRegions::VarCoeffMap &varcoeffs)
	Update the supplied variable coefficients. More...

unsigned int	GetWspSize ()
	Get the size of the required workspace. More...

unsigned int	GetNumElmt ()
	Get number of elements. More...

StdRegions::StdExpansionSharedPtr	GetExpSharedPtr ()
	Get expansion pointer. More...

unsigned int	GetInputSize (bool defaultIn=true)

unsigned int	GetOutputSize (bool defaultOut=true)

Protected Attributes
Array< TwoD, const NekDouble >	m_derivFac

Array< OneD, const NekDouble >	m_jac

int	m_dim

int	m_coordim

StdRegions::FactorMap	m_factors

StdRegions::VarCoeffMap	m_varcoeffs

NekDouble	m_lambda

bool	m_HasVarCoeffDiff

Array< OneD, Array< OneD, NekDouble > >	m_diff

const StdRegions::ConstFactorType	m_factorCoeffDef [3][3]

Protected Attributes inherited from Nektar::Collections::Operator
bool	m_isDeformed

StdRegions::StdExpansionSharedPtr	m_stdExp

unsigned int	m_numElmt
	number of elements that the operator is applied on More...

unsigned int	m_nqe

unsigned int	m_wspSize

unsigned int	m_inputSize
	number of modes or quadrature points that are passed as input to an operator More...

unsigned int	m_outputSize
	number of modes or quadrature points that are taken as output from an operator More...

unsigned int	m_inputSizeOther
	Number of modes or quadrature points, opposite to m_inputSize. More...

unsigned int	m_outputSizeOther
	Number of modes or quadrature points, opposite to m_outputSize. More...

Private Member Functions
	Helmholtz_IterPerExp (vector< StdRegions::StdExpansionSharedPtr > pCollExp, CoalescedGeomDataSharedPtr pGeomData, StdRegions::FactorMap factors)

Private Member Functions inherited from Nektar::Collections::Helmholtz_Helper
	Helmholtz_Helper ()

Detailed Description

Helmholtz operator using LocalRegions implementation.

Definition at line 191 of file library/Collections/Helmholtz.cpp.

Constructor & Destructor Documentation

◆ ~Helmholtz_IterPerExp()

Nektar::Collections::Helmholtz_IterPerExp::~Helmholtz_IterPerExp ( )

finaldefault

◆ Helmholtz_IterPerExp()

Nektar::Collections::Helmholtz_IterPerExp::Helmholtz_IterPerExp	(	vector< StdRegions::StdExpansionSharedPtr >	pCollExp,
		CoalescedGeomDataSharedPtr	pGeomData,
		StdRegions::FactorMap	factors
	)

inlineprivate

Definition at line 528 of file library/Collections/Helmholtz.cpp.

        : Operator(pCollExp, pGeomData, factors), Helmholtz_Helper()
    {
        m_dim     = pCollExp[0]->GetShapeDimension();
        m_coordim = pCollExp[0]->GetCoordim();
        int nqtot = m_stdExp->GetTotPoints();
 
        m_derivFac = pGeomData->GetDerivFactors(pCollExp);
        m_jac      = pGeomData->GetJac(pCollExp);
        m_wspSize  = (2 * m_coordim + 1) * nqtot;
 
        m_lambda          = 1.0;
        m_HasVarCoeffDiff = false;
        m_factors         = StdRegions::NullFactorMap;
        m_varcoeffs       = StdRegions::NullVarCoeffMap;
        this->UpdateFactors(factors);
    }

References m_coordim, m_derivFac, m_dim, m_factors, m_HasVarCoeffDiff, m_jac, m_lambda, Nektar::Collections::Operator::m_stdExp, m_varcoeffs, Nektar::Collections::Operator::m_wspSize, Nektar::StdRegions::NullFactorMap, Nektar::StdRegions::NullVarCoeffMap, and UpdateFactors().

Member Function Documentation

◆ operator()() [1/2]

void Nektar::Collections::Helmholtz_IterPerExp::operator()	(	const Array< OneD, const NekDouble > &	input,
		Array< OneD, NekDouble > &	output0,
		Array< OneD, NekDouble > &	output1,
		Array< OneD, NekDouble > &	output2,
		Array< OneD, NekDouble > &	wsp
	)

inlinefinalvirtual

Perform operation.

Implements Nektar::Collections::Operator.

Definition at line 198 of file library/Collections/Helmholtz.cpp.

    {
        const int nCoeffs = m_stdExp->GetNcoeffs();
        const int nPhys   = m_stdExp->GetTotPoints();
 
        ASSERTL1(input.size() >= m_numElmt * nCoeffs,
                 "input array size is insufficient");
        ASSERTL1(output.size() >= m_numElmt * nCoeffs,
                 "output array size is insufficient");
 
        Array<OneD, NekDouble> tmpphys, t1;
        Array<OneD, Array<OneD, NekDouble>> dtmp(3);
        Array<OneD, Array<OneD, NekDouble>> tmp(3);
 
        tmpphys = wsp;
        for (int i = 1; i < m_coordim + 1; ++i)
        {
            dtmp[i - 1] = wsp + i * nPhys;
            tmp[i - 1]  = wsp + (i + m_coordim) * nPhys;
        }
 
        for (int i = 0; i < m_numElmt; ++i)
        {
            m_stdExp->BwdTrans(input + i * nCoeffs, tmpphys);
 
            // local derivative
            m_stdExp->PhysDeriv(tmpphys, dtmp[0], dtmp[1], dtmp[2]);
 
            // determine mass matrix term
            if (m_isDeformed)
            {
                Vmath::Vmul(nPhys, m_jac + i * nPhys, 1, tmpphys, 1, tmpphys,
                            1);
            }
            else
            {
                Vmath::Smul(nPhys, m_jac[i], tmpphys, 1, tmpphys, 1);
            }
 
            m_stdExp->IProductWRTBase(tmpphys, t1 = output + i * nCoeffs);
            Vmath::Smul(nCoeffs, m_lambda, output + i * nCoeffs, 1,
                        t1 = output + i * nCoeffs, 1);
 
            if (m_isDeformed)
            {
                // calculate full derivative
                for (int j = 0; j < m_coordim; ++j)
                {
                    Vmath::Vmul(nPhys,
                                m_derivFac[j * m_dim].origin() + i * nPhys, 1,
                                &dtmp[0][0], 1, &tmp[j][0], 1);
 
                    for (int k = 1; k < m_dim; ++k)
                    {
                        Vmath::Vvtvp(
                            nPhys,
                            m_derivFac[j * m_dim + k].origin() + i * nPhys, 1,
                            &dtmp[k][0], 1, &tmp[j][0], 1, &tmp[j][0], 1);
                    }
                }
 
                if (m_HasVarCoeffDiff)
                {
                    // calculate dtmp[i] = dx/dxi sum_j diff[0][j] tmp[j]
                    //                   + dy/dxi sum_j diff[1][j] tmp[j]
                    //                   + dz/dxi sum_j diff[2][j] tmp[j]
 
                    // First term
                    Vmath::Smul(nPhys, m_diff[0][0], &tmp[0][0], 1, &tmpphys[0],
                                1);
                    for (int l = 1; l < m_coordim; ++l)
                    {
                        Vmath::Svtvp(nPhys, m_diff[0][l], &tmp[l][0], 1,
                                     &tmpphys[0], 1, &tmpphys[0], 1);
                    }
 
                    for (int j = 0; j < m_dim; ++j)
                    {
                        Vmath::Vmul(nPhys, m_derivFac[j].origin() + i * nPhys,
                                    1, &tmpphys[0], 1, &dtmp[j][0], 1);
                    }
 
                    // Second and third terms
                    for (int k = 1; k < m_coordim; ++k)
                    {
                        Vmath::Smul(nPhys, m_diff[k][0], &tmp[0][0], 1,
                                    &tmpphys[0], 1);
                        for (int l = 1; l < m_coordim; ++l)
                        {
                            Vmath::Svtvp(nPhys, m_diff[k][l], &tmp[l][0], 1,
                                         &tmpphys[0], 1, &tmpphys[0], 1);
                        }
 
                        for (int j = 0; j < m_dim; ++j)
                        {
                            Vmath::Vvtvp(nPhys,
                                         m_derivFac[j + k * m_dim].origin() +
                                             i * nPhys,
                                         1, &tmpphys[0], 1, &dtmp[j][0], 1,
                                         &dtmp[j][0], 1);
                        }
                    }
                }
                else
                {
                    // calculate dx/dxi tmp[0] + dy/dxi tmp[1]
                    //                         + dz/dxi tmp[2]
                    for (int j = 0; j < m_dim; ++j)
                    {
                        Vmath::Vmul(nPhys, m_derivFac[j].origin() + i * nPhys,
                                    1, &tmp[0][0], 1, &dtmp[j][0], 1);
 
                        for (int k = 1; k < m_coordim; ++k)
                        {
                            Vmath::Vvtvp(nPhys,
                                         m_derivFac[j + k * m_dim].origin() +
                                             i * nPhys,
                                         1, &tmp[k][0], 1, &dtmp[j][0], 1,
                                         &dtmp[j][0], 1);
                        }
                    }
                }
 
                // calculate Iproduct WRT Std Deriv
                for (int j = 0; j < m_dim; ++j)
                {
                    // multiply by Jacobian
                    Vmath::Vmul(nPhys, m_jac + i * nPhys, 1, dtmp[j], 1,
                                dtmp[j], 1);
 
                    m_stdExp->IProductWRTDerivBase(j, dtmp[j], tmp[0]);
                    Vmath::Vadd(nCoeffs, tmp[0], 1, output + i * nCoeffs, 1,
                                t1 = output + i * nCoeffs, 1);
                }
            }
            else
            {
                // calculate full derivative
                for (int j = 0; j < m_coordim; ++j)
                {
                    Vmath::Smul(nPhys, m_derivFac[j * m_dim][i], &dtmp[0][0], 1,
                                &tmp[j][0], 1);
 
                    for (int k = 1; k < m_dim; ++k)
                    {
                        Vmath::Svtvp(nPhys, m_derivFac[j * m_dim + k][i],
                                     &dtmp[k][0], 1, &tmp[j][0], 1, &tmp[j][0],
                                     1);
                    }
                }
 
                if (m_HasVarCoeffDiff)
                {
                    // calculate dtmp[i] = dx/dxi sum_j diff[0][j] tmp[j]
                    //                   + dy/dxi sum_j diff[1][j] tmp[j]
                    //                   + dz/dxi sum_j diff[2][j] tmp[j]
 
                    // First term
                    Vmath::Smul(nPhys, m_diff[0][0], &tmp[0][0], 1, &tmpphys[0],
                                1);
                    for (int l = 1; l < m_coordim; ++l)
                    {
                        Vmath::Svtvp(nPhys, m_diff[0][l], &tmp[l][0], 1,
                                     &tmpphys[0], 1, &tmpphys[0], 1);
                    }
 
                    for (int j = 0; j < m_dim; ++j)
                    {
                        Vmath::Smul(nPhys, m_derivFac[j][i], &tmpphys[0], 1,
                                    &dtmp[j][0], 1);
                    }
 
                    // Second and third terms
                    for (int k = 1; k < m_coordim; ++k)
                    {
                        Vmath::Smul(nPhys, m_diff[k][0], &tmp[0][0], 1,
                                    &tmpphys[0], 1);
                        for (int l = 1; l < m_coordim; ++l)
                        {
                            Vmath::Svtvp(nPhys, m_diff[k][l], &tmp[l][0], 1,
                                         &tmpphys[0], 1, &tmpphys[0], 1);
                        }
 
                        for (int j = 0; j < m_dim; ++j)
                        {
                            Vmath::Svtvp(nPhys, m_derivFac[j + k * m_dim][i],
                                         &tmpphys[0], 1, &dtmp[j][0], 1,
                                         &dtmp[j][0], 1);
                        }
                    }
                }
                else
                {
                    // calculate dx/dxi tmp[0] + dy/dxi tmp[2]
                    //                         + dz/dxi tmp[3]
                    for (int j = 0; j < m_dim; ++j)
                    {
                        Vmath::Smul(nPhys, m_derivFac[j][i], &tmp[0][0], 1,
                                    &dtmp[j][0], 1);
 
                        for (int k = 1; k < m_coordim; ++k)
                        {
                            Vmath::Svtvp(nPhys, m_derivFac[j + k * m_dim][i],
                                         &tmp[k][0], 1, &dtmp[j][0], 1,
                                         &dtmp[j][0], 1);
                        }
                    }
                }
 
                // calculate Iproduct WRT Std Deriv
                for (int j = 0; j < m_dim; ++j)
                {
                    // multiply by Jacobian
                    Vmath::Smul(nPhys, m_jac[i], dtmp[j], 1, dtmp[j], 1);
 
                    m_stdExp->IProductWRTDerivBase(j, dtmp[j], tmp[0]);
                    Vmath::Vadd(nCoeffs, tmp[0], 1, output + i * nCoeffs, 1,
                                t1 = output + i * nCoeffs, 1);
                }
            }
        }
    }

References ASSERTL1, m_coordim, m_derivFac, m_diff, m_dim, m_HasVarCoeffDiff, Nektar::Collections::Operator::m_isDeformed, m_jac, m_lambda, Nektar::Collections::Operator::m_numElmt, Nektar::Collections::Operator::m_stdExp, Nektar::MovementTests::origin, Vmath::Smul(), Vmath::Svtvp(), Vmath::Vadd(), Vmath::Vmul(), and Vmath::Vvtvp().

◆ operator()() [2/2]

void Nektar::Collections::Helmholtz_IterPerExp::operator()	(	int	dir,
		const Array< OneD, const NekDouble > &	input,
		Array< OneD, NekDouble > &	output,
		Array< OneD, NekDouble > &	wsp
	)

inlinefinalvirtual

Implements Nektar::Collections::Operator.

Definition at line 425 of file library/Collections/Helmholtz.cpp.

    {
        NEKERROR(ErrorUtil::efatal, "Not valid for this operator.");
    }

References Nektar::ErrorUtil::efatal, and NEKERROR.

◆ UpdateFactors()

void Nektar::Collections::Helmholtz_IterPerExp::UpdateFactors ( StdRegions::FactorMap factors )

inlineoverridevirtual

Check the validity of supplied constant factors.

Parameters

factors Map of factors

Reimplemented from Nektar::Collections::Operator.

Definition at line 438 of file library/Collections/Helmholtz.cpp.

    {
        // If match previous factors, nothing to do.
        if (m_factors == factors)
        {
            return;
        }
 
        m_factors = factors;
 
        // Check Lambda constant of Helmholtz operator
        auto x = factors.find(StdRegions::eFactorLambda);
        ASSERTL1(
            x != factors.end(),
            "Constant factor not defined: " +
                std::string(
                    StdRegions::ConstFactorTypeMap[StdRegions::eFactorLambda]));
        m_lambda = x->second;
 
        // If varcoeffs not supplied, nothing else to do.
        m_HasVarCoeffDiff = false;
        auto d            = factors.find(StdRegions::eFactorCoeffD00);
        if (d == factors.end())
        {
            return;
        }
 
        m_diff = Array<OneD, Array<OneD, NekDouble>>(m_coordim);
        for (int i = 0; i < m_coordim; ++i)
        {
            m_diff[i] = Array<OneD, NekDouble>(m_coordim, 0.0);
        }
 
        for (int i = 0; i < m_coordim; ++i)
        {
            d = factors.find(m_factorCoeffDef[i][i]);
            if (d != factors.end())
            {
                m_diff[i][i] = d->second;
            }
            else
            {
                m_diff[i][i] = 1.0;
            }
 
            for (int j = i + 1; j < m_coordim; ++j)
            {
                d = factors.find(m_factorCoeffDef[i][j]);
                if (d != factors.end())
                {
                    m_diff[i][j] = m_diff[j][i] = d->second;
                }
            }
        }
        m_HasVarCoeffDiff = true;
    }

References ASSERTL1, Nektar::StdRegions::ConstFactorTypeMap, Nektar::UnitTests::d(), Nektar::StdRegions::eFactorCoeffD00, Nektar::StdRegions::eFactorLambda, Nektar::VarcoeffHashingTest::factors, m_coordim, m_diff, m_factorCoeffDef, m_factors, m_HasVarCoeffDiff, and m_lambda.

Referenced by Helmholtz_IterPerExp().

◆ UpdateVarcoeffs()

void Nektar::Collections::Helmholtz_IterPerExp::UpdateVarcoeffs ( StdRegions::VarCoeffMap & varcoeffs )

inlineoverridevirtual

Check the validity of supplied variable coefficients. Note that the m_varcoeffs are not implemented for the IterPerExp operator. There exists a specialised implementation for variable diffusion in the above function UpdateFactors.

Parameters

varcoeffs Map of varcoeffs

Reimplemented from Nektar::Collections::Operator.

Definition at line 504 of file library/Collections/Helmholtz.cpp.

    {
        m_varcoeffs = varcoeffs;
    }

References m_varcoeffs.

Member Data Documentation

◆ m_coordim

int Nektar::Collections::Helmholtz_IterPerExp::m_coordim

protected

Definition at line 513 of file library/Collections/Helmholtz.cpp.

Referenced by Helmholtz_IterPerExp(), operator()(), and UpdateFactors().

◆ m_derivFac

Array<TwoD, const NekDouble> Nektar::Collections::Helmholtz_IterPerExp::m_derivFac

protected

Definition at line 510 of file library/Collections/Helmholtz.cpp.

Referenced by Helmholtz_IterPerExp(), and operator()().

◆ m_diff

Array<OneD, Array<OneD, NekDouble> > Nektar::Collections::Helmholtz_IterPerExp::m_diff

protected

Definition at line 518 of file library/Collections/Helmholtz.cpp.

Referenced by operator()(), and UpdateFactors().

◆ m_dim

int Nektar::Collections::Helmholtz_IterPerExp::m_dim

protected

Definition at line 512 of file library/Collections/Helmholtz.cpp.

Referenced by Helmholtz_IterPerExp(), and operator()().

◆ m_factorCoeffDef

const StdRegions::ConstFactorType Nektar::Collections::Helmholtz_IterPerExp::m_factorCoeffDef[3][3]

protected

Initial value:

= {
        {StdRegions::eFactorCoeffD00, StdRegions::eFactorCoeffD01,
         StdRegions::eFactorCoeffD02},
        {StdRegions::eFactorCoeffD01, StdRegions::eFactorCoeffD11,
         StdRegions::eFactorCoeffD12},
        {StdRegions::eFactorCoeffD02, StdRegions::eFactorCoeffD12,
         StdRegions::eFactorCoeffD22}}