Nektar::FilterMovingBody Class Reference

#include <FilterMovingBody.h>

Inheritance diagram for Nektar::FilterMovingBody:

Collaboration diagram for Nektar::FilterMovingBody:

Public Member Functions
	FilterMovingBody (const LibUtilities::SessionReaderSharedPtr &pSession, const ParamMap &pParams)
	~FilterMovingBody ()
virtual void	v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
virtual void	v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
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)
virtual void	v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
virtual bool	v_IsTimeDependent ()
Public Member Functions inherited from Nektar::SolverUtils::Filter
SOLVER_UTILS_EXPORT	Filter (const LibUtilities::SessionReaderSharedPtr &pSession)
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 ParamMap &pParams)
	Creates an instance of this class. More...

Static Public Attributes
static std::string	className

Private Attributes
LibUtilities::SessionReaderSharedPtr	m_session
vector< unsigned int >	m_boundaryRegionsIdList
	ID's of boundary regions where we want the forces. More...
vector< bool >	m_boundaryRegionIsInList
	Determines if a given Boundary Region is in m_boundaryRegionsIdList. More...
unsigned int	m_index_f
unsigned int	m_index_m
unsigned int	m_outputFrequency
unsigned int	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
int	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 Attributes inherited from Nektar::SolverUtils::Filter
LibUtilities::SessionReaderSharedPtr	m_session

Detailed Description

Definition at line 51 of file FilterMovingBody.h.

Constructor & Destructor Documentation

Nektar::FilterMovingBody::FilterMovingBody	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const ParamMap &	pParams
	)

Definition at line 53 of file FilterMovingBody.cpp.

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

    : Filter(pSession),
      m_session(pSession)
{
    ParamMap::const_iterator it;

    // OutputFile
    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, it->second);
        m_outputFrequency = floor(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;
}

Nektar::FilterMovingBody::~FilterMovingBody

(

)

Definition at line 113 of file FilterMovingBody.cpp.

{

}

Member Function Documentation

static SolverUtils::FilterSharedPtr Nektar::FilterMovingBody::create ( const LibUtilities::SessionReaderSharedPtr &  pSession, const ParamMap &  pParams  )

inlinestatic

Creates an instance of this class.

Definition at line 57 of file FilterMovingBody.h.

References Nektar::MemoryManager< DataType >::AllocateSharedPtr().

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

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 226 of file FilterMovingBody.cpp.

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

{
    int n, cnt, elmtid, nq, offset, boundary;
    int nt  = pFields[0]->GetNpoints();
    int dim = pFields.num_elements()-1;

    StdRegions::StdExpansionSharedPtr 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
    int 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(int i = 0; i < pFields.num_elements(); ++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();
    int local_planes = ZIDs.num_elements();

    // Homogeneous 1D case  Compute forces on all WALL boundaries
    // This only has to be done on the zero (mean) Fourier mode.
    for(int 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.num_elements(); ++n)
        {
            if(m_boundaryRegionIsInList[n] == 1)
            {
                for(int i = 0; i <  BndExp[n]->GetExpSize(); ++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(int 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
                    int nbc = bc->GetTotPoints();

                    // several vectors for computing the forces
                    Array<OneD, NekDouble> Pb(nbc,0.0);

                    for(int 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->GetEdgePhysVals(boundary,bc,P,Pb);

                    for(int j = 0; j < dim; ++j)
                    {
                        elmt->GetEdgePhysVals(boundary,bc,gradU[j],fgradU[j]);
                        elmt->GetEdgePhysVals(boundary,bc,gradV[j],fgradV[j]);
                    }

                    //normals of the element
                    const Array<OneD, Array<OneD, NekDouble> > &normals
                                            = elmt->GetEdgeNormal(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(int i = 0; i < pFields.num_elements(); ++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(int 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(int 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(int 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(int 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(int 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(int 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);

        int 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(int 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;
        }

        int colrank = vColComm->GetRank();
        int nstrips;

        NekDouble DistStrip;

        pSession->LoadParameter("Strip_Z", nstrips);
        pSession->LoadParameter("DistStrip", DistStrip);

        Array<OneD, NekDouble> z_coords(nstrips);
        for(int 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(int 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(int i = 1; i < nstrips; i++)
            {
                if(colrank == i)
                {
                    vColComm->Send(0, fces);
                }
            }
        }
    }
}

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 697 of file FilterMovingBody.cpp.

References m_index_m, m_outputFrequency, m_outputStream, m_session, and npts.

{
    // Only output every m_outputFrequency.
    if ((m_index_m++) % m_outputFrequency)
    {
        return;
    }

    int 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(int 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;
    }
}

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

virtual

Implements Nektar::SolverUtils::Filter.

Definition at line 752 of file FilterMovingBody.cpp.

References m_outputStream.

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

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

virtual

Implements Nektar::SolverUtils::Filter.

Definition at line 122 of file FilterMovingBody.cpp.

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 m_session.

{
    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.c_str(),
                                                   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
    int cnt;

    unsigned int numBoundaryRegions
                    = pFields[0]->GetBndConditions().num_elements();

    m_boundaryRegionIsInList.insert(m_boundaryRegionIsInList.end(),
                                    numBoundaryRegions, 0);

    SpatialDomains::BoundaryConditions bcs(m_session,pFields[0]->GetGraph());

    const SpatialDomains::BoundaryRegionCollection &bregions
                    = bcs.GetBoundaryRegions();

    SpatialDomains::BoundaryRegionCollection::const_iterator it;

    for (cnt = 0, it = bregions.begin(); it != bregions.end();
            ++it, cnt++)
    {
        if ( std::find(m_boundaryRegionsIdList.begin(),
                       m_boundaryRegionsIdList.end(), it->first) !=
                m_boundaryRegionsIdList.end() )
        {
            m_boundaryRegionIsInList[cnt] = 1;
        }
    }

    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;
    }
}

bool Nektar::FilterMovingBody::v_IsTimeDependent ( )

virtual

Implements Nektar::SolverUtils::Filter.

Definition at line 767 of file FilterMovingBody.cpp.

{
    return true;
}

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

inlinevirtual

Implements Nektar::SolverUtils::Filter.

Definition at line 76 of file FilterMovingBody.h.

                                   {}

Friends And Related Function Documentation

friend class MemoryManager< FilterMovingBody >

friend

Definition at line 54 of file FilterMovingBody.h.

Member Data Documentation

std::string Nektar::FilterMovingBody::className

static

Initial value:

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

Definition at line 65 of file FilterMovingBody.h.

vector<bool> Nektar::FilterMovingBody::m_boundaryRegionIsInList

private

Determines if a given Boundary Region is in m_boundaryRegionsIdList.

Definition at line 105 of file FilterMovingBody.h.

Referenced by UpdateForce(), and v_Initialise().

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

private

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

Definition at line 102 of file FilterMovingBody.h.

Referenced by v_Initialise().

std::string Nektar::FilterMovingBody::m_BoundaryString

private

Definition at line 114 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

LibUtilities::BasisSharedPtr Nektar::FilterMovingBody::m_homogeneousBasis

private

Definition at line 113 of file FilterMovingBody.h.

unsigned int Nektar::FilterMovingBody::m_index_f

private

Definition at line 106 of file FilterMovingBody.h.

Referenced by UpdateForce(), and v_Initialise().

unsigned int Nektar::FilterMovingBody::m_index_m

private

Definition at line 107 of file FilterMovingBody.h.

Referenced by UpdateMotion(), and v_Initialise().

bool Nektar::FilterMovingBody::m_isHomogeneous1D

private

Definition at line 112 of file FilterMovingBody.h.

Referenced by FilterMovingBody().

std::string Nektar::FilterMovingBody::m_outputFile_fce

private

Definition at line 118 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

std::string Nektar::FilterMovingBody::m_outputFile_mot

private

Definition at line 119 of file FilterMovingBody.h.

Referenced by FilterMovingBody(), and v_Initialise().

unsigned int Nektar::FilterMovingBody::m_outputFrequency

private

Definition at line 108 of file FilterMovingBody.h.

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

unsigned int Nektar::FilterMovingBody::m_outputPlane

private

plane to take history point from if using a homogeneous1D expansion

Definition at line 111 of file FilterMovingBody.h.

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

private

Definition at line 117 of file FilterMovingBody.h.

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

int Nektar::FilterMovingBody::m_planes

private

number of planes for homogeneous1D expansion

Definition at line 116 of file FilterMovingBody.h.

LibUtilities::SessionReaderSharedPtr Nektar::FilterMovingBody::m_session

private

Definition at line 99 of file FilterMovingBody.h.

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