#include <FilterMovingBody.h>

Inheritance diagram for Nektar::FilterMovingBody:

Public Member Functions
void	UpdateForce (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &Aeroforces, const NekDouble &time)

void	UpdateMotion (const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, Array< OneD, NekDouble > &MotionVars, const NekDouble &time)

Public Member Functions inherited from Nektar::SolverUtils::Filter
SOLVER_UTILS_EXPORT	Filter (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< EquationSystem > &pEquation)

virtual SOLVER_UTILS_EXPORT	~Filter ()

SOLVER_UTILS_EXPORT void	Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)

SOLVER_UTILS_EXPORT void	Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)

SOLVER_UTILS_EXPORT void	Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)

SOLVER_UTILS_EXPORT bool	IsTimeDependent ()

SOLVER_UTILS_EXPORT std::string	SetupOutput (const std::string ext, const ParamMap &pParams)

SOLVER_UTILS_EXPORT std::string	SetupOutput (const std::string ext, const std::string inname)

SOLVER_UTILS_EXPORT void	SetUpdateOnInitialise (bool flag)

Static Public Member Functions
static SolverUtils::FilterSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)
	Creates an instance of this class. More...

Static Public Attributes
static std::string	className

Private Member Functions
	FilterMovingBody (const LibUtilities::SessionReaderSharedPtr &pSession, const std::shared_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)

	~FilterMovingBody () override=default

void	v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override

void	v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override

void	v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override

bool	v_IsTimeDependent () override

Private Attributes
std::vector< unsigned int >	m_boundaryRegionsIdList
	ID's of boundary regions where we want the forces. More...

std::vector< bool >	m_boundaryRegionIsInList
	Determines if a given Boundary Region is in m_boundaryRegionsIdList. More...

size_t	m_index_f

size_t	m_index_m

size_t	m_outputFrequency

size_t	m_outputPlane
	plane to take history point from if using a homogeneous1D expansion More...

bool	m_isHomogeneous1D

LibUtilities::BasisSharedPtr	m_homogeneousBasis

std::string	m_BoundaryString

size_t	m_planes
	number of planes for homogeneous1D expansion More...

Array< OneD, std::ofstream >	m_outputStream

std::string	m_outputFile_fce

std::string	m_outputFile_mot

Friends
class	MemoryManager< FilterMovingBody >

Additional Inherited Members
Public Types inherited from Nektar::SolverUtils::Filter
typedef std::map< std::string, std::string >	ParamMap

Protected Member Functions inherited from Nektar::SolverUtils::Filter
virtual void	v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0

virtual void	v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0

virtual void	v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)=0

virtual bool	v_IsTimeDependent ()=0

virtual SOLVER_UTILS_EXPORT std::string	v_SetupOutput (const std::string ext, const ParamMap &pParams)

virtual SOLVER_UTILS_EXPORT std::string	v_SetupOutput (const std::string ext, const std::string inname)

Protected Attributes inherited from Nektar::SolverUtils::Filter
LibUtilities::SessionReaderSharedPtr	m_session

const std::weak_ptr< EquationSystem >	m_equ

bool	m_updateOnInitialise = true

Detailed Description

Definition at line 50 of file FilterMovingBody.h.

Constructor & Destructor Documentation

◆ FilterMovingBody()

Nektar::FilterMovingBody::FilterMovingBody	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::shared_ptr< SolverUtils::EquationSystem > &	pEquation,
		const ParamMap &	pParams
	)

private

Definition at line 53 of file FilterMovingBody.cpp.

    : Filter(pSession, pEquation)
{
    // OutputFile
    // forces
    std::string ext;
    ext              = ".fce";
    m_outputFile_fce = Filter::SetupOutput(ext, pParams);
 
    // Motion
    ext              = ".mot";
    m_outputFile_mot = Filter::SetupOutput(ext, pParams);
 
    // OutputFrequency
    auto it = pParams.find("OutputFrequency");
    if (it == pParams.end())
    {
        m_outputFrequency = 1;
    }
    else
    {
        LibUtilities::Equation equ(m_session->GetInterpreter(), it->second);
        m_outputFrequency = round(equ.Evaluate());
    }
 
    pSession->MatchSolverInfo("Homogeneous", "1D", m_isHomogeneous1D, false);
    ASSERTL0(m_isHomogeneous1D, "Moving Body implemented just for 3D "
                                "Homogeneous 1D discetisations.");
 
    // Boundary (to calculate the forces)
    it = pParams.find("Boundary");
    ASSERTL0(it != pParams.end(), "Missing parameter 'Boundary'.");
    ASSERTL0(it->second.length() > 0, "Missing parameter 'Boundary'.");
    m_BoundaryString = it->second;
}

References ASSERTL0, Nektar::LibUtilities::Equation::Evaluate(), m_BoundaryString, m_isHomogeneous1D, m_outputFile_fce, m_outputFile_mot, m_outputFrequency, Nektar::SolverUtils::Filter::m_session, and Nektar::SolverUtils::Filter::SetupOutput().

◆ ~FilterMovingBody()

Nektar::FilterMovingBody::~FilterMovingBody ( )

overrideprivatedefault

Member Function Documentation

◆ create()

static SolverUtils::FilterSharedPtr Nektar::FilterMovingBody::create	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::shared_ptr< SolverUtils::EquationSystem > &	pEquation,
		const ParamMap &	pParams
	)

inlinestatic

Creates an instance of this class.

Definition at line 56 of file FilterMovingBody.h.

    {
        SolverUtils::FilterSharedPtr p =
            MemoryManager<FilterMovingBody>::AllocateSharedPtr(
                pSession, pEquation, pParams);
        return p;
    }

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), and CellMLToNektar.cellml_metadata::p.

◆ UpdateForce()

void Nektar::FilterMovingBody::UpdateForce	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const Array< OneD, const MultiRegions::ExpListSharedPtr > &	pFields,
		Array< OneD, NekDouble > &	Aeroforces,
		const NekDouble &	time
	)

Definition at line 197 of file FilterMovingBody.cpp.

{
    size_t n, cnt, elmtid, nq, offset, boundary;
    size_t nt  = pFields[0]->GetNpoints();
    size_t dim = pFields.size() - 1;
 
    LocalRegions::ExpansionSharedPtr elmt;
    Array<OneD, int> BoundarytoElmtID;
    Array<OneD, int> BoundarytoTraceID;
    Array<OneD, MultiRegions::ExpListSharedPtr> BndExp;
 
    Array<OneD, const NekDouble> P(nt);
    Array<OneD, const NekDouble> U(nt);
    Array<OneD, const NekDouble> V(nt);
    Array<OneD, const NekDouble> W(nt);
 
    Array<OneD, Array<OneD, NekDouble>> gradU(dim);
    Array<OneD, Array<OneD, NekDouble>> gradV(dim);
    Array<OneD, Array<OneD, NekDouble>> gradW(dim);
 
    Array<OneD, Array<OneD, NekDouble>> fgradU(dim);
    Array<OneD, Array<OneD, NekDouble>> fgradV(dim);
    Array<OneD, Array<OneD, NekDouble>> fgradW(dim);
 
    LibUtilities::CommSharedPtr vComm     = pFields[0]->GetComm();
    LibUtilities::CommSharedPtr vRowComm  = vComm->GetRowComm();
    LibUtilities::CommSharedPtr vColComm  = vComm->GetColumnComm();
    LibUtilities::CommSharedPtr vCColComm = vColComm->GetColumnComm();
 
    // set up storage space for forces on all the planes (z-positions)
    // on each processors some of them will remain empty as we may
    // have just few planes per processor
    size_t Num_z_pos = pFields[0]->GetHomogeneousBasis()->GetNumModes();
    Array<OneD, NekDouble> Fx(Num_z_pos, 0.0);
    Array<OneD, NekDouble> Fxp(Num_z_pos, 0.0);
    Array<OneD, NekDouble> Fxv(Num_z_pos, 0.0);
    Array<OneD, NekDouble> Fy(Num_z_pos, 0.0);
    Array<OneD, NekDouble> Fyp(Num_z_pos, 0.0);
    Array<OneD, NekDouble> Fyv(Num_z_pos, 0.0);
 
    NekDouble rho = (pSession->DefinesParameter("rho"))
                        ? (pSession->GetParameter("rho"))
                        : 1;
    NekDouble mu  = rho * pSession->GetParameter("Kinvis");
 
    for (size_t i = 0; i < pFields.size(); ++i)
    {
        pFields[i]->SetWaveSpace(false);
        pFields[i]->BwdTrans(pFields[i]->GetCoeffs(), pFields[i]->UpdatePhys());
        pFields[i]->SetPhysState(true);
    }
 
    // Get the number of local planes on the process and their IDs
    // to properly locate the forces in the Fx, Fy etc. vectors.
    Array<OneD, unsigned int> ZIDs;
    ZIDs                = pFields[0]->GetZIDs();
    size_t local_planes = ZIDs.size();
 
    // Homogeneous 1D case  Compute forces on all WALL boundaries
    // This only has to be done on the zero (mean) Fourier mode.
    for (size_t plane = 0; plane < local_planes; plane++)
    {
        pFields[0]->GetPlane(plane)->GetBoundaryToElmtMap(BoundarytoElmtID,
                                                          BoundarytoTraceID);
        BndExp = pFields[0]->GetPlane(plane)->GetBndCondExpansions();
        StdRegions::StdExpansionSharedPtr bc;
 
        // loop over the types of boundary conditions
        for (cnt = n = 0; n < BndExp.size(); ++n)
        {
            if (m_boundaryRegionIsInList[n] == 1)
            {
                size_t nExp = BndExp[n]->GetExpSize();
                for (size_t i = 0; i < nExp; ++i, cnt++)
                {
                    // find element of this expansion.
                    elmtid = BoundarytoElmtID[cnt];
                    elmt   = pFields[0]->GetPlane(plane)->GetExp(elmtid);
                    nq     = elmt->GetTotPoints();
                    offset =
                        pFields[0]->GetPlane(plane)->GetPhys_Offset(elmtid);
 
                    // Initialise local arrays for the velocity
                    // gradients size of total number of quadrature
                    // points for each element (hence local).
                    for (size_t j = 0; j < dim; ++j)
                    {
                        gradU[j] = Array<OneD, NekDouble>(nq, 0.0);
                        gradV[j] = Array<OneD, NekDouble>(nq, 0.0);
                        gradW[j] = Array<OneD, NekDouble>(nq, 0.0);
                    }
 
                    // identify boundary of element
                    boundary = BoundarytoTraceID[cnt];
 
                    // Extract  fields
                    U = pFields[0]->GetPlane(plane)->GetPhys() + offset;
                    V = pFields[1]->GetPlane(plane)->GetPhys() + offset;
                    P = pFields[3]->GetPlane(plane)->GetPhys() + offset;
 
                    // compute the gradients
                    elmt->PhysDeriv(U, gradU[0], gradU[1]);
                    elmt->PhysDeriv(V, gradV[0], gradV[1]);
 
                    // Get face 1D expansion from element expansion
                    bc = BndExp[n]->GetExp(i)->as<LocalRegions::Expansion1D>();
 
                    // number of points on the boundary
                    size_t nbc = bc->GetTotPoints();
 
                    // several vectors for computing the forces
                    Array<OneD, NekDouble> Pb(nbc, 0.0);
 
                    for (size_t j = 0; j < dim; ++j)
                    {
                        fgradU[j] = Array<OneD, NekDouble>(nbc, 0.0);
                        fgradV[j] = Array<OneD, NekDouble>(nbc, 0.0);
                    }
 
                    Array<OneD, NekDouble> drag_t(nbc, 0.0);
                    Array<OneD, NekDouble> lift_t(nbc, 0.0);
                    Array<OneD, NekDouble> drag_p(nbc, 0.0);
                    Array<OneD, NekDouble> lift_p(nbc, 0.0);
                    Array<OneD, NekDouble> temp(nbc, 0.0);
                    Array<OneD, NekDouble> temp2(nbc, 0.0);
 
                    // identify boundary of element .
                    boundary = BoundarytoTraceID[cnt];
 
                    // extraction of the pressure and wss on the
                    // boundary of the element
                    elmt->GetTracePhysVals(boundary, bc, P, Pb);
 
                    for (size_t j = 0; j < dim; ++j)
                    {
                        elmt->GetTracePhysVals(boundary, bc, gradU[j],
                                               fgradU[j]);
                        elmt->GetTracePhysVals(boundary, bc, gradV[j],
                                               fgradV[j]);
                    }
 
                    // normals of the element
                    const Array<OneD, Array<OneD, NekDouble>> &normals =
                        elmt->GetTraceNormal(boundary);
 
                    //
                    // Compute viscous tractive forces on wall from
                    //
                    //  t_i  = - T_ij * n_j  (minus sign for force
                    //                        exerted BY fluid ON wall),
                    //
                    // where
                    //
                    //  T_ij = viscous stress tensor (here in Cartesian
                    //         coords)
                    //                          dU_i    dU_j
                    //       = RHO * KINVIS * ( ----  + ---- ) .
                    //                          dx_j    dx_i
 
                    // a) DRAG TERMS
                    //-rho*kinvis*(2*du/dx*nx+(du/dy+dv/dx)*ny)
 
                    Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, drag_t, 1);
                    Vmath::Vmul(nbc, drag_t, 1, normals[1], 1, drag_t, 1);
 
                    Vmath::Smul(nbc, 2.0, fgradU[0], 1, fgradU[0], 1);
                    Vmath::Vmul(nbc, fgradU[0], 1, normals[0], 1, temp2, 1);
                    Vmath::Smul(nbc, 0.5, fgradU[0], 1, fgradU[0], 1);
 
                    Vmath::Vadd(nbc, temp2, 1, drag_t, 1, drag_t, 1);
                    Vmath::Smul(nbc, -mu, drag_t, 1, drag_t, 1);
 
                    // zero temporary storage vector
                    Vmath::Zero(nbc, temp, 0.0);
                    Vmath::Zero(nbc, temp2, 0.0);
 
                    // b) LIFT TERMS
                    //-rho*kinvis*(2*dv/dy*ny+(du/dy+dv/dx)*nx)
 
                    Vmath::Vadd(nbc, fgradU[1], 1, fgradV[0], 1, lift_t, 1);
                    Vmath::Vmul(nbc, lift_t, 1, normals[0], 1, lift_t, 1);
 
                    Vmath::Smul(nbc, 2.0, fgradV[1], 1, fgradV[1], 1);
                    Vmath::Vmul(nbc, fgradV[1], 1, normals[1], 1, temp2, 1);
                    Vmath::Smul(nbc, -0.5, fgradV[1], 1, fgradV[1], 1);
 
                    Vmath::Vadd(nbc, temp2, 1, lift_t, 1, lift_t, 1);
                    Vmath::Smul(nbc, -mu, lift_t, 1, lift_t, 1);
 
                    // Compute normal tractive forces on all WALL
                    // boundaries
 
                    Vmath::Vvtvp(nbc, Pb, 1, normals[0], 1, drag_p, 1, drag_p,
                                 1);
                    Vmath::Vvtvp(nbc, Pb, 1, normals[1], 1, lift_p, 1, lift_p,
                                 1);
 
                    // integration over the boundary
                    Fxv[ZIDs[plane]] += bc->Integral(drag_t);
                    Fyv[ZIDs[plane]] += bc->Integral(lift_t);
 
                    Fxp[ZIDs[plane]] += bc->Integral(drag_p);
                    Fyp[ZIDs[plane]] += bc->Integral(lift_p);
                }
            }
            else
            {
                cnt += BndExp[n]->GetExpSize();
            }
        }
    }
 
    for (size_t i = 0; i < pFields.size(); ++i)
    {
        pFields[i]->SetWaveSpace(true);
        pFields[i]->BwdTrans(pFields[i]->GetCoeffs(), pFields[i]->UpdatePhys());
        pFields[i]->SetPhysState(false);
    }
 
    // In case we are using an hybrid parallelisation we need
    // to reduce the forces on the same plane coming from
    // different mesh partitions.
    // It is quite an expensive communication, therefore
    // we check to make sure it is actually required.
 
    if (vComm->GetRowComm()->GetSize() > 0)
    {
        // NOTE 1: We can eventually sum the viscous and pressure
        // component before doing the communication, thus
        // reducing by a factor 2 the communication.
        // NOTE 2: We may want to set up in the Comm class an AllReduce
        // routine wich can handle arrays more efficiently
        for (size_t plane = 0; plane < local_planes; plane++)
        {
            vRowComm->AllReduce(Fxp[ZIDs[plane]], LibUtilities::ReduceSum);
            vRowComm->AllReduce(Fxv[ZIDs[plane]], LibUtilities::ReduceSum);
            vRowComm->AllReduce(Fyp[ZIDs[plane]], LibUtilities::ReduceSum);
            vRowComm->AllReduce(Fyv[ZIDs[plane]], LibUtilities::ReduceSum);
        }
    }
 
    // At this point on rank (0,n) of the Mesh partion communicator we have
    // the total areo forces on the planes which are on the same column
    // communicator. Since the planes are scattered on different processors
    // some of the entries of the vector Fxp, Fxp etc. are still zero.
    // Now we need to reduce the values on a single vector on rank (0,0) of the
    // global communicator.
    if (!pSession->DefinesSolverInfo("HomoStrip"))
    {
        if (vComm->GetRowComm()->GetRank() == 0)
        {
            for (size_t z = 0; z < Num_z_pos; z++)
            {
                vColComm->AllReduce(Fxp[z], LibUtilities::ReduceSum);
                vColComm->AllReduce(Fxv[z], LibUtilities::ReduceSum);
                vColComm->AllReduce(Fyp[z], LibUtilities::ReduceSum);
                vColComm->AllReduce(Fyv[z], LibUtilities::ReduceSum);
            }
        }
    }
    else
    {
        if (vComm->GetRowComm()->GetRank() == 0)
        {
            for (size_t z = 0; z < Num_z_pos; z++)
            {
                vCColComm->AllReduce(Fxp[z], LibUtilities::ReduceSum);
                vCColComm->AllReduce(Fxv[z], LibUtilities::ReduceSum);
                vCColComm->AllReduce(Fyp[z], LibUtilities::ReduceSum);
                vCColComm->AllReduce(Fyv[z], LibUtilities::ReduceSum);
            }
        }
    }
 
    if (!pSession->DefinesSolverInfo("HomoStrip"))
    {
        // set the forces imparted on the cable's wall
        for (size_t plane = 0; plane < local_planes; plane++)
        {
            Aeroforces[plane] = Fxp[ZIDs[plane]] + Fxv[ZIDs[plane]];
            Aeroforces[plane + local_planes] =
                Fyp[ZIDs[plane]] + Fyv[ZIDs[plane]];
        }
 
        // Only output every m_outputFrequency.
        if ((m_index_f++) % m_outputFrequency)
        {
            return;
        }
 
        // At thi point in rank (0,0) we have the full vectors
        // containing Fxp,Fxv,Fyp and Fyv where different positions
        // in the vectors correspond to different planes.
        // Here we write it to file. We do it just on one porcess
 
        Array<OneD, NekDouble> z_coords(Num_z_pos, 0.0);
        Array<OneD, const NekDouble> pts =
            pFields[0]->GetHomogeneousBasis()->GetZ();
 
        NekDouble LZ;
        pSession->LoadParameter("LZ", LZ);
        Vmath::Smul(Num_z_pos, LZ / 2.0, pts, 1, z_coords, 1);
        Vmath::Sadd(Num_z_pos, LZ / 2.0, z_coords, 1, z_coords, 1);
        if (vComm->GetRank() == 0)
        {
 
            Vmath::Vadd(Num_z_pos, Fxp, 1, Fxv, 1, Fx, 1);
            Vmath::Vadd(Num_z_pos, Fyp, 1, Fyv, 1, Fy, 1);
 
            for (size_t i = 0; i < Num_z_pos; i++)
            {
                m_outputStream[0].width(8);
                m_outputStream[0] << std::setprecision(6) << time;
 
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(6) << z_coords[i];
 
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << Fxp[i];
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << Fxv[i];
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << Fx[i];
 
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << Fyp[i];
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << Fyv[i];
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << Fy[i];
                m_outputStream[0] << std::endl;
            }
        }
    }
    else
    {
        // average the forces over local planes for each processor
        Array<OneD, NekDouble> fces(6, 0.0);
        for (size_t plane = 0; plane < local_planes; plane++)
        {
            fces[0] += Fxp[ZIDs[plane]] + Fxv[ZIDs[plane]];
            fces[1] += Fyp[ZIDs[plane]] + Fyv[ZIDs[plane]];
            fces[2] += Fxp[ZIDs[plane]];
            fces[3] += Fyp[ZIDs[plane]];
            fces[4] += Fxv[ZIDs[plane]];
            fces[5] += Fyv[ZIDs[plane]];
        }
 
        fces[0] = fces[0] / local_planes;
        fces[1] = fces[1] / local_planes;
        fces[2] = fces[2] / local_planes;
        fces[3] = fces[3] / local_planes;
        fces[4] = fces[4] / local_planes;
        fces[5] = fces[5] / local_planes;
 
        // average the forces over communicators within each strip
        vCColComm->AllReduce(fces[0], LibUtilities::ReduceSum);
        vCColComm->AllReduce(fces[1], LibUtilities::ReduceSum);
        vCColComm->AllReduce(fces[2], LibUtilities::ReduceSum);
        vCColComm->AllReduce(fces[3], LibUtilities::ReduceSum);
        vCColComm->AllReduce(fces[4], LibUtilities::ReduceSum);
        vCColComm->AllReduce(fces[5], LibUtilities::ReduceSum);
 
        size_t npts = vComm->GetColumnComm()->GetColumnComm()->GetSize();
 
        fces[0] = fces[0] / npts;
        fces[1] = fces[1] / npts;
        fces[2] = fces[2] / npts;
        fces[3] = fces[3] / npts;
        fces[4] = fces[4] / npts;
        fces[5] = fces[5] / npts;
 
        for (size_t plane = 0; plane < local_planes; plane++)
        {
            Aeroforces[plane]                = fces[0];
            Aeroforces[plane + local_planes] = fces[1];
        }
 
        // Only output every m_outputFrequency.
        if ((m_index_f) % m_outputFrequency)
        {
            return;
        }
 
        size_t colrank = vColComm->GetRank();
        size_t nstrips;
 
        NekDouble DistStrip;
 
        pSession->LoadParameter("Strip_Z", nstrips);
        pSession->LoadParameter("DistStrip", DistStrip);
 
        Array<OneD, NekDouble> z_coords(nstrips);
        for (size_t i = 0; i < nstrips; i++)
        {
            z_coords[i] = i * DistStrip;
        }
 
        if (colrank == 0)
        {
            m_outputStream[0].width(8);
            m_outputStream[0] << std::setprecision(6) << time;
 
            m_outputStream[0].width(15);
            m_outputStream[0] << std::setprecision(6) << z_coords[0];
 
            m_outputStream[0].width(15);
            m_outputStream[0] << std::setprecision(8) << fces[2];
            m_outputStream[0].width(15);
            m_outputStream[0] << std::setprecision(8) << fces[4];
            m_outputStream[0].width(15);
            m_outputStream[0] << std::setprecision(8) << fces[0];
 
            m_outputStream[0].width(15);
            m_outputStream[0] << std::setprecision(8) << fces[3];
            m_outputStream[0].width(15);
            m_outputStream[0] << std::setprecision(8) << fces[5];
            m_outputStream[0].width(15);
            m_outputStream[0] << std::setprecision(8) << fces[1];
            m_outputStream[0] << std::endl;
 
            for (size_t i = 1; i < nstrips; i++)
            {
                vColComm->Recv(i, fces);
 
                m_outputStream[0].width(8);
                m_outputStream[0] << std::setprecision(6) << time;
 
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(6) << z_coords[i];
 
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << fces[2];
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << fces[4];
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << fces[0];
 
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << fces[3];
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << fces[5];
                m_outputStream[0].width(15);
                m_outputStream[0] << std::setprecision(8) << fces[1];
                m_outputStream[0] << std::endl;
            }
        }
        else
        {
            for (size_t i = 1; i < nstrips; i++)
            {
                if (colrank == i)
                {
                    vColComm->Send(0, fces);
                }
            }
        }
    }
}

References Nektar::StdRegions::StdExpansion::GetTotPoints(), m_boundaryRegionIsInList, m_index_f, m_outputFrequency, m_outputStream, Nektar::LibUtilities::P, Nektar::LibUtilities::ReduceSum, Vmath::Sadd(), Vmath::Smul(), Vmath::Vadd(), Vmath::Vmul(), Vmath::Vvtvp(), Nektar::UnitTests::z(), and Vmath::Zero().

◆ UpdateMotion()

void Nektar::FilterMovingBody::UpdateMotion	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const Array< OneD, const MultiRegions::ExpListSharedPtr > &	pFields,
		Array< OneD, NekDouble > &	MotionVars,
		const NekDouble &	time
	)

Definition at line 663 of file FilterMovingBody.cpp.

{
    // Only output every m_outputFrequency.
    if ((m_index_m++) % m_outputFrequency)
    {
        return;
    }
 
    size_t npts;
    // Length of the cable
    NekDouble Length;
 
    if (!pSession->DefinesSolverInfo("HomoStrip"))
    {
        pSession->LoadParameter("LZ", Length);
        npts = m_session->GetParameter("HomModesZ");
    }
    else
    {
        pSession->LoadParameter("LC", Length);
        npts = m_session->GetParameter("HomStructModesZ");
    }
 
    NekDouble z_coords;
    for (size_t n = 0; n < npts; n++)
    {
        z_coords = Length / npts * n;
        m_outputStream[1].width(8);
        m_outputStream[1] << std::setprecision(6) << time;
        m_outputStream[1].width(15);
        m_outputStream[1] << std::setprecision(6) << z_coords;
        m_outputStream[1].width(15);
        m_outputStream[1] << std::setprecision(8) << MotionVars[n];
        m_outputStream[1].width(15);
        m_outputStream[1] << std::setprecision(8) << MotionVars[npts + n];
        m_outputStream[1].width(15);
        m_outputStream[1] << std::setprecision(8) << MotionVars[2 * npts + n];
        m_outputStream[1].width(15);
        m_outputStream[1] << std::setprecision(8) << MotionVars[3 * npts + n];
        m_outputStream[1].width(15);
        m_outputStream[1] << std::setprecision(8) << MotionVars[4 * npts + n];
        m_outputStream[1].width(15);
        m_outputStream[1] << std::setprecision(8) << MotionVars[5 * npts + n];
        m_outputStream[1] << std::endl;
    }
}

References m_index_m, m_outputFrequency, m_outputStream, and Nektar::SolverUtils::Filter::m_session.

◆ v_Finalise()

void Nektar::FilterMovingBody::v_Finalise	(	const Array< OneD, const MultiRegions::ExpListSharedPtr > &	pFields,
		const NekDouble &	time
	)

overrideprivatevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 717 of file FilterMovingBody.cpp.

{
    if (pFields[0]->GetComm()->GetRank() == 0)
    {
        m_outputStream[0].close();
        m_outputStream[1].close();
    }
}

References m_outputStream.

◆ v_Initialise()

void Nektar::FilterMovingBody::v_Initialise	(	const Array< OneD, const MultiRegions::ExpListSharedPtr > &	pFields,
		const NekDouble &	time
	)

overrideprivatevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 95 of file FilterMovingBody.cpp.

{
    m_index_f      = 0;
    m_index_m      = 0;
    m_outputStream = Array<OneD, std::ofstream>(2);
    // Parse the boundary regions into a list.
    std::string::size_type FirstInd = m_BoundaryString.find_first_of('[') + 1;
    std::string::size_type LastInd  = m_BoundaryString.find_last_of(']') - 1;
 
    ASSERTL0(FirstInd <= LastInd,
             (std::string("Error reading boundary region definition:") +
              m_BoundaryString)
                 .c_str());
 
    std::string IndString =
        m_BoundaryString.substr(FirstInd, LastInd - FirstInd + 1);
 
    bool parseGood =
        ParseUtils::GenerateSeqVector(IndString, m_boundaryRegionsIdList);
 
    ASSERTL0(parseGood && !m_boundaryRegionsIdList.empty(),
             (std::string("Unable to read boundary regions index "
                          "range for FilterAeroForces: ") +
              IndString)
                 .c_str());
 
    // determine what boundary regions need to be considered
    size_t numBoundaryRegions = pFields[0]->GetBndConditions().size();
 
    m_boundaryRegionIsInList.insert(m_boundaryRegionIsInList.end(),
                                    numBoundaryRegions, 0);
 
    SpatialDomains::BoundaryConditions bcs(m_session, pFields[0]->GetGraph());
 
    const SpatialDomains::BoundaryRegionCollection &bregions =
        bcs.GetBoundaryRegions();
 
    size_t cnt = 0;
    for (auto &it : bregions)
    {
        if (std::find(m_boundaryRegionsIdList.begin(),
                      m_boundaryRegionsIdList.end(),
                      it.first) != m_boundaryRegionsIdList.end())
        {
            m_boundaryRegionIsInList[cnt] = 1;
        }
        ++cnt;
    }
 
    LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
 
    if (vComm->GetRank() == 0)
    {
        // Open output stream for cable forces
        m_outputStream[0].open(m_outputFile_fce.c_str());
        m_outputStream[0] << "#";
        m_outputStream[0].width(7);
        m_outputStream[0] << "Time";
        m_outputStream[0].width(15);
        m_outputStream[0] << "z";
        m_outputStream[0].width(15);
        m_outputStream[0] << "Fx (press)";
        m_outputStream[0].width(15);
        m_outputStream[0] << "Fx (visc)";
        m_outputStream[0].width(15);
        m_outputStream[0] << "Fx (tot)";
        m_outputStream[0].width(15);
        m_outputStream[0] << "Fy (press)";
        m_outputStream[0].width(15);
        m_outputStream[0] << "Fy (visc)";
        m_outputStream[0].width(15);
        m_outputStream[0] << "Fy (tot)";
        m_outputStream[0] << std::endl;
 
        // Open output stream for cable motions
        m_outputStream[1].open(m_outputFile_mot.c_str());
        m_outputStream[1] << "#";
        m_outputStream[1].width(7);
        m_outputStream[1] << "Time";
        m_outputStream[1].width(15);
        m_outputStream[1] << "z";
        m_outputStream[1].width(15);
        m_outputStream[1] << "Disp_x";
        m_outputStream[1].width(15);
        m_outputStream[1] << "Vel_x";
        m_outputStream[1].width(15);
        m_outputStream[1] << "Acel_x";
        m_outputStream[1].width(15);
        m_outputStream[1] << "Disp_y";
        m_outputStream[1].width(15);
        m_outputStream[1] << "Vel_y";
        m_outputStream[1].width(15);
        m_outputStream[1] << "Acel_y";
        m_outputStream[1] << std::endl;
    }
}

References ASSERTL0, Nektar::StdRegions::find(), Nektar::ParseUtils::GenerateSeqVector(), Nektar::SpatialDomains::BoundaryConditions::GetBoundaryRegions(), m_boundaryRegionIsInList, m_boundaryRegionsIdList, m_BoundaryString, m_index_f, m_index_m, m_outputFile_fce, m_outputFile_mot, m_outputStream, and Nektar::SolverUtils::Filter::m_session.

◆ v_IsTimeDependent()

bool Nektar::FilterMovingBody::v_IsTimeDependent ( )

overrideprivatevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 731 of file FilterMovingBody.cpp.

{
    return true;
}

◆ v_Update()

void Nektar::FilterMovingBody::v_Update	(	const Array< OneD, const MultiRegions::ExpListSharedPtr > &	pFields,
		const NekDouble &	time
	)

inlineoverrideprivatevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 111 of file FilterMovingBody.h.

115 {

116 }

Friends And Related Function Documentation

◆ MemoryManager< FilterMovingBody >

friend class MemoryManager< FilterMovingBody >

friend

Definition at line 48 of file FilterMovingBody.h.

Member Data Documentation

◆ className

std::string Nektar::FilterMovingBody::className

static

Initial value:

=
    SolverUtils::GetFilterFactory().RegisterCreatorFunction(
        "MovingBody", FilterMovingBody::create, "Moving Body Filter")

Definition at line 67 of file FilterMovingBody.h.

◆ m_boundaryRegionIsInList

std::vector<bool> Nektar::FilterMovingBody::m_boundaryRegionIsInList

private

Determines if a given Boundary Region is in m_boundaryRegionsIdList.

Definition at line 84 of file FilterMovingBody.h.

Referenced by UpdateForce(), and v_Initialise().

◆ m_boundaryRegionsIdList

std::vector<unsigned int> Nektar::FilterMovingBody::m_boundaryRegionsIdList

private

ID's of boundary regions where we want the forces.

Definition at line 81 of file FilterMovingBody.h.

Referenced by v_Initialise().

◆ m_BoundaryString

std::string Nektar::FilterMovingBody::m_BoundaryString

private

Definition at line 93 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

◆ m_homogeneousBasis

LibUtilities::BasisSharedPtr Nektar::FilterMovingBody::m_homogeneousBasis

private

Definition at line 92 of file FilterMovingBody.h.

◆ m_index_f

size_t Nektar::FilterMovingBody::m_index_f

private

Definition at line 85 of file FilterMovingBody.h.

Referenced by UpdateForce(), and v_Initialise().

◆ m_index_m

size_t Nektar::FilterMovingBody::m_index_m

private

Definition at line 86 of file FilterMovingBody.h.

Referenced by UpdateMotion(), and v_Initialise().

◆ m_isHomogeneous1D

bool Nektar::FilterMovingBody::m_isHomogeneous1D

private

Definition at line 91 of file FilterMovingBody.h.

Referenced by FilterMovingBody().

◆ m_outputFile_fce

std::string Nektar::FilterMovingBody::m_outputFile_fce

private

Definition at line 97 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

◆ m_outputFile_mot

std::string Nektar::FilterMovingBody::m_outputFile_mot

private

Definition at line 98 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

◆ m_outputFrequency

size_t Nektar::FilterMovingBody::m_outputFrequency

private

Definition at line 87 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), UpdateForce(), and UpdateMotion().

◆ m_outputPlane

size_t Nektar::FilterMovingBody::m_outputPlane

private

plane to take history point from if using a homogeneous1D expansion

Definition at line 90 of file FilterMovingBody.h.

◆ m_outputStream

Array<OneD, std::ofstream> Nektar::FilterMovingBody::m_outputStream

private

Definition at line 96 of file FilterMovingBody.h.

Referenced by UpdateForce(), UpdateMotion(), v_Finalise(), and v_Initialise().

◆ m_planes

size_t Nektar::FilterMovingBody::m_planes

private

number of planes for homogeneous1D expansion

Definition at line 95 of file FilterMovingBody.h.

Public Member Functions

Static Public Member Functions

Static Public Attributes

Private Member Functions

Private Attributes

Friends

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ FilterMovingBody()

◆ ~FilterMovingBody()

Member Function Documentation

◆ create()

◆ UpdateForce()

◆ UpdateMotion()

◆ v_Finalise()

◆ v_Initialise()

◆ v_IsTimeDependent()

◆ v_Update()

Friends And Related Function Documentation

◆ MemoryManager< FilterMovingBody >

Member Data Documentation

◆ className

◆ m_boundaryRegionIsInList

◆ m_boundaryRegionsIdList

◆ m_BoundaryString

◆ m_homogeneousBasis

◆ m_index_f

◆ m_index_m

◆ m_isHomogeneous1D

◆ m_outputFile_fce

◆ m_outputFile_mot

◆ m_outputFrequency

◆ m_outputPlane

◆ m_outputStream

◆ m_planes