FilterMovingBody.cpp

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////////
//
// File: FilterMovingBody.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Output moving body forces during time-stepping.
//
///////////////////////////////////////////////////////////////////////////////
 
#include <IncNavierStokesSolver/Filters/FilterMovingBody.h>
#include <LibUtilities/BasicUtils/ParseUtils.h>
#include <LibUtilities/Memory/NekMemoryManager.hpp>
#include <LocalRegions/Expansion1D.h>
#include <LocalRegions/Expansion2D.h>
#include <LocalRegions/Expansion3D.h>
#include <iomanip>
 
using namespace std;
 
namespace Nektar
{
 
std::string FilterMovingBody::className =
    SolverUtils::GetFilterFactory().RegisterCreatorFunction(
        "MovingBody", FilterMovingBody::create, "Moving Body Filter");
/**
 *
 */
FilterMovingBody::FilterMovingBody(
    const LibUtilities::SessionReaderSharedPtr &pSession,
    const std::weak_ptr<SolverUtils::EquationSystem> &pEquation,
    const ParamMap &pParams)
    : 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;
}
 
/**
 *
 */
FilterMovingBody::~FilterMovingBody()
{
}
 
/**
 *
 */
void FilterMovingBody::v_Initialise(
    const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
    const NekDouble &time)
{
    boost::ignore_unused(time);
 
    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;
    }
}
 
/**
 *
 */
void FilterMovingBody::UpdateForce(
    const LibUtilities::SessionReaderSharedPtr &pSession,
    const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
    Array<OneD, NekDouble> &Aeroforces, const NekDouble &time)
{
    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);
                }
            }
        }
    }
}
 
/**
 *
 */
void FilterMovingBody::UpdateMotion(
    const LibUtilities::SessionReaderSharedPtr &pSession,
    const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
    Array<OneD, NekDouble> &MotionVars, const NekDouble &time)
{
    boost::ignore_unused(pFields);
 
    // 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;
    }
}
 
/**
 *
 */
void FilterMovingBody::v_Finalise(
    const Array<OneD, const MultiRegions::ExpListSharedPtr> &pFields,
    const NekDouble &time)
{
    boost::ignore_unused(time);
 
    if (pFields[0]->GetComm()->GetRank() == 0)
    {
        m_outputStream[0].close();
        m_outputStream[1].close();
    }
}
 
/**
 *
 */
bool FilterMovingBody::v_IsTimeDependent()
{
    return true;
}
} // namespace Nektar