Nektar::FilterMovingBody Class Reference

#include <FilterMovingBody.h>

Inheritance diagram for Nektar::FilterMovingBody:

Public Member Functions
	FilterMovingBody (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_ptr< SolverUtils::EquationSystem > &pEquation, const ParamMap &pParams)
	~FilterMovingBody () override
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	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)
void	v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time) override
bool	v_IsTimeDependent () override
Public Member Functions inherited from Nektar::SolverUtils::Filter
SOLVER_UTILS_EXPORT	Filter (const LibUtilities::SessionReaderSharedPtr &pSession, const std::weak_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 ()

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

Static Public Attributes
static std::string	className

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
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
Protected Attributes inherited from Nektar::SolverUtils::Filter
LibUtilities::SessionReaderSharedPtr	m_session
const std::weak_ptr< EquationSystem >	m_equ

Detailed Description

Definition at line 50 of file FilterMovingBody.h.

Constructor & Destructor Documentation

FilterMovingBody()

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

Definition at line 54 of file FilterMovingBody.cpp.

    : Filter(pSession, pEquation)
{
    // OutputFile
    auto it = pParams.find("OutputFile");
    if (it == pParams.end())
    {
        m_outputFile_fce = pSession->GetSessionName();
        m_outputFile_mot = pSession->GetSessionName();
    }
    else
    {
        ASSERTL0(it->second.length() > 0, "Missing parameter 'OutputFile'.");
 
        m_outputFile_fce = it->second;
        m_outputFile_mot = it->second;
    }
    if (!(m_outputFile_fce.length() >= 4 &&
          m_outputFile_fce.substr(m_outputFile_fce.length() - 4) == ".fce"))
    {
        m_outputFile_fce += ".fce";
    }
    if (!(m_outputFile_mot.length() >= 4 &&
          m_outputFile_mot.substr(m_outputFile_mot.length() - 4) == ".mot"))
    {
        m_outputFile_mot += ".mot";
    }
 
    // OutputFrequency
    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, and Nektar::SolverUtils::Filter::m_session.

~FilterMovingBody()

Nektar::FilterMovingBody::~FilterMovingBody ( )

override

Definition at line 111 of file FilterMovingBody.cpp.

{
}

Member Function Documentation

create()

static SolverUtils::FilterSharedPtr Nektar::FilterMovingBody::create ( const LibUtilities::SessionReaderSharedPtr &  pSession, const std::weak_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 220 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] << setprecision(6) << time;
 
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(6) << z_coords[i];
 
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << Fxp[i];
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << Fxv[i];
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << Fx[i];
 
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << Fyp[i];
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << Fyv[i];
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << Fy[i];
                m_outputStream[0] << 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] << setprecision(6) << time;
 
            m_outputStream[0].width(15);
            m_outputStream[0] << setprecision(6) << z_coords[0];
 
            m_outputStream[0].width(15);
            m_outputStream[0] << setprecision(8) << fces[2];
            m_outputStream[0].width(15);
            m_outputStream[0] << setprecision(8) << fces[4];
            m_outputStream[0].width(15);
            m_outputStream[0] << setprecision(8) << fces[0];
 
            m_outputStream[0].width(15);
            m_outputStream[0] << setprecision(8) << fces[3];
            m_outputStream[0].width(15);
            m_outputStream[0] << setprecision(8) << fces[5];
            m_outputStream[0].width(15);
            m_outputStream[0] << setprecision(8) << fces[1];
            m_outputStream[0] << endl;
 
            for (size_t i = 1; i < nstrips; i++)
            {
                vColComm->Recv(i, fces);
 
                m_outputStream[0].width(8);
                m_outputStream[0] << setprecision(6) << time;
 
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(6) << z_coords[i];
 
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << fces[2];
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << fces[4];
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << fces[0];
 
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << fces[3];
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << fces[5];
                m_outputStream[0].width(15);
                m_outputStream[0] << setprecision(8) << fces[1];
                m_outputStream[0] << 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 686 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] << setprecision(6) << time;
        m_outputStream[1].width(15);
        m_outputStream[1] << setprecision(6) << z_coords;
        m_outputStream[1].width(15);
        m_outputStream[1] << setprecision(8) << MotionVars[n];
        m_outputStream[1].width(15);
        m_outputStream[1] << setprecision(8) << MotionVars[npts + n];
        m_outputStream[1].width(15);
        m_outputStream[1] << setprecision(8) << MotionVars[2 * npts + n];
        m_outputStream[1].width(15);
        m_outputStream[1] << setprecision(8) << MotionVars[3 * npts + n];
        m_outputStream[1].width(15);
        m_outputStream[1] << setprecision(8) << MotionVars[4 * npts + n];
        m_outputStream[1].width(15);
        m_outputStream[1] << setprecision(8) << MotionVars[5 * npts + n];
        m_outputStream[1] << 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  )

overridevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 740 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  )

overridevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 118 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] << 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] << 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 ( )

overridevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 754 of file FilterMovingBody.cpp.

{
    return true;
}

v_Update()

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

inlineoverridevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 79 of file FilterMovingBody.h.

    {
    }

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 107 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 104 of file FilterMovingBody.h.

Referenced by v_Initialise().

m_BoundaryString

std::string Nektar::FilterMovingBody::m_BoundaryString

private

Definition at line 116 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

m_homogeneousBasis

LibUtilities::BasisSharedPtr Nektar::FilterMovingBody::m_homogeneousBasis

private

Definition at line 115 of file FilterMovingBody.h.

m_index_f

size_t Nektar::FilterMovingBody::m_index_f

private

Definition at line 108 of file FilterMovingBody.h.

Referenced by UpdateForce(), and v_Initialise().

m_index_m

size_t Nektar::FilterMovingBody::m_index_m

private

Definition at line 109 of file FilterMovingBody.h.

Referenced by UpdateMotion(), and v_Initialise().

m_isHomogeneous1D

bool Nektar::FilterMovingBody::m_isHomogeneous1D

private

Definition at line 114 of file FilterMovingBody.h.

Referenced by FilterMovingBody().

m_outputFile_fce

std::string Nektar::FilterMovingBody::m_outputFile_fce

private

Definition at line 120 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

m_outputFile_mot

std::string Nektar::FilterMovingBody::m_outputFile_mot

private

Definition at line 121 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

m_outputFrequency

size_t Nektar::FilterMovingBody::m_outputFrequency

private

Definition at line 110 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 113 of file FilterMovingBody.h.

m_outputStream

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

private

Definition at line 119 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 118 of file FilterMovingBody.h.