#include <RigidSolver.h>

Public Member Functions
void	InitObject (const LibUtilities::SessionReaderSharedPtr session, const MultiRegions::ExpListSharedPtr &pField, Array< OneD, NekDouble > pivot)

void	SetMovableDoFs (std::vector< bool > &moveDoFs)

void	SetInitialConditions (const LibUtilities::SessionReaderSharedPtr session, Array< OneD, NekDouble > MRFData)

void	UpdateFrameVelocity (Array< OneD, NekDouble > &aeroforce, const NekDouble &time, Array< OneD, NekDouble > &MRFData)
	Updates the forcing array with the current required forcing.

void	SetNewmarkBetaSolver (Array< OneD, NekDouble > &AddedMass)

void	GetFilterInfo (const LibUtilities::SessionReaderSharedPtr session, std::map< std::string, std::string > &vParams)

bool	HasFreeDoFs ()

void	SetOldFvis (Array< OneD, NekDouble > force)

Protected Member Functions
void	LoadParameters (const LibUtilities::SessionReaderSharedPtr session, const TiXmlElement *pSolver)

void	InitBodySolver (const LibUtilities::SessionReaderSharedPtr session, const TiXmlElement *pSolver, Array< OneD, NekDouble > pivot)

NekDouble	EvaluateExpression (const LibUtilities::SessionReaderSharedPtr session, std::string expression)

void	UpdatePrescribed (const NekDouble &time, std::map< int, NekDouble > &Dirs)

void	UpdateMRFData (Array< OneD, NekDouble > &MRFData)

void	SetInitialConditions (std::map< int, NekDouble > &Dirs)

void	SolveBodyMotion (Array< OneD, Array< OneD, NekDouble > > &bodyVel, const Array< OneD, NekDouble > &forcebody, std::map< int, NekDouble > &Dirs)

void	ParserFunctionToMap (const bool allowzero, const LibUtilities::SessionReaderSharedPtr session, const std::string &FuncName, const std::string &var, std::map< int, LibUtilities::EquationSharedPtr > &result, const int i)

void	CheckParameters ()

Private Member Functions
void	SolveInertialFrame (Array< OneD, Array< OneD, NekDouble > > &bodyVel, const Array< OneD, NekDouble > &forcebody, std::map< int, NekDouble > &Dirs)

void	SolveBodyFrame (Array< OneD, Array< OneD, NekDouble > > &bodyVel, const Array< OneD, NekDouble > &forcebody, std::map< int, NekDouble > &Dirs)

void	SolveRotOneFree (Array< OneD, Array< OneD, NekDouble > > &bodyVel, const Array< OneD, NekDouble > &forcebody, std::map< int, NekDouble > &Dirs)

Private Attributes
int	m_solveType

int	m_index

NekDouble	m_currentTime

NekDouble	m_timestep

NekDouble	m_beta

NekDouble	m_gamma

int	m_spacedim

bool	m_isRoot

bool	m_hasFreeMotion

bool	m_hasRotation

int	m_outputFrequency

std::ofstream	m_outputStream

std::set< int >	m_dirDoFs

Array< OneD, NekDouble >	m_inertialPosition

Array< OneD, Array< OneD, NekDouble > >	m_vel

std::map< int, LibUtilities::EquationSharedPtr >	m_extForceFunction

std::map< int, LibUtilities::EquationSharedPtr >	m_frameVelFunction

Array< OneD, NekDouble >	m_extForceXYZ

Array< OneD, NekDouble >	m_pivot

NekDouble	m_pivotdistance

Array< OneD, NekDouble >	m_oldFvis

NekDouble	m_mass

NekDouble	m_rotaionInertia

Array< OneD, NekDouble >	m_M

Array< OneD, NekDouble >	m_C

Array< OneD, NekDouble >	m_K

Newmark_BetaSolver	m_bodySolver

FrameTransform	m_frame

Detailed Description

Definition at line 106 of file RigidSolver.h.

Member Function Documentation

◆ CheckParameters()

void Nektar::RigidSolver::CheckParameters ( )

protected

Definition at line 83 of file RigidSolver.cpp.

{
    // Count free translational DoFs
    int nFreeTrans = 0;
    bool freeX     = false;
 
    for (int i = 0; i < m_spacedim; ++i)
    {
        if (m_dirDoFs.find(i) == m_dirDoFs.end())
        {
            ++nFreeTrans;
            if (i == 0)
            {
                freeX = true;
            }
        }
    }
    if (!m_hasRotation)
    {
        m_solveType = 0; // inertial frame
    }
    else if (nFreeTrans == 1 && freeX)
    {
        m_solveType = 1; // with given rotation and only x free
        if (m_spacedim >= 2)
        {
            m_dirDoFs.erase(1);
        }
    }
    else if ((nFreeTrans == 2 || nFreeTrans == 0) && m_hasRotation)
    {
        m_solveType = 2; // full planar rigid solver (rotation involved)
    }
    else
    {
        ASSERTL0(false, "Unsupported rigid motion: this solver supports only "
                        "(i) translation-only motion, "
                        "(ii) x-only translation with optional rotation, or "
                        "(iii) full planar rigid-body motion.");
    }
    if (m_hasRotation)
    {
        m_inertialPosition = Array<OneD, NekDouble>(2, 0.);
        // m_K should be zero
        NekDouble sum = 0;
        for (size_t i = 0; i < m_K.size(); ++i)
        {
            sum += fabs(m_K[i]);
        }
        for (size_t i = 0; i < m_C.size(); ++i)
        {
            sum += fabs(m_C[i]);
        }
        ASSERTL0(sum < NekConstants::kNekZeroTol,
                 "C, K should be zero [in body frame solver].");
    }
}

References ASSERTL0, Nektar::NekConstants::kNekZeroTol, m_C, m_dirDoFs, m_hasRotation, m_inertialPosition, m_K, m_solveType, and m_spacedim.

Referenced by InitObject().

◆ EvaluateExpression()

NekDouble Nektar::RigidSolver::EvaluateExpression	(	const LibUtilities::SessionReaderSharedPtr	session,
		std::string	expression
	)

protected

Definition at line 173 of file RigidSolver.cpp.

{
    NekDouble value = 0.;
    try
    {
        LibUtilities::Equation expession(session->GetInterpreter(), expression);
        value = expession.Evaluate();
    }
    catch (const std::runtime_error &)
    {
        NEKERROR(ErrorUtil::efatal, "Error evaluating expression" + expression);
    }
    return value;
}

References Nektar::ErrorUtil::efatal, Nektar::LibUtilities::Equation::Evaluate(), and NEKERROR.

Referenced by InitBodySolver(), and LoadParameters().

◆ GetFilterInfo()

void Nektar::RigidSolver::GetFilterInfo	(	const LibUtilities::SessionReaderSharedPtr	session,
		std::map< std::string, std::string > &	vParams
	)

Definition at line 982 of file RigidSolver.cpp.

{
    TiXmlElement *pForce      = session->GetElement("Nektar/RIGIDSOLVER");
    const TiXmlElement *param = pForce->FirstChildElement("BOUNDARY");
    ASSERTL0(param, "Body surface should be assigned");
 
    vParams["Boundary"]           = param->GetText();
    const TiXmlElement *pivotElmt = pForce->FirstChildElement("PIVOTPOINT");
    if (pivotElmt)
    {
        std::string pstr = pivotElmt->GetText();
        std::replace(pstr.begin(), pstr.end(), ',', ' ');
        vParams["MomentPoint"] = pstr;
    }
}

References ASSERTL0.

Referenced by Nektar::VCSFSI::InitialiseFilter().

◆ HasFreeDoFs()

bool Nektar::RigidSolver::HasFreeDoFs ( )

inline

Definition at line 122 of file RigidSolver.h.

    {
        return m_hasFreeMotion;
    };

References m_hasFreeMotion.

Referenced by Nektar::VCSFSI::InitialiseFilter(), and Nektar::VCSFSI::v_SolveSolid().

◆ InitBodySolver()

void Nektar::RigidSolver::InitBodySolver	(	const LibUtilities::SessionReaderSharedPtr	session,
		const TiXmlElement *	pSolver,
		Array< OneD, NekDouble >	pivot
	)

protected

Definition at line 329 of file RigidSolver.cpp.

{
    int NumDof = m_spacedim + 1;
    const TiXmlElement *mssgTag;
    std::string mssgStr;
    // allocate memory and initialise
    m_vel = Array<OneD, Array<OneD, NekDouble>>(3);
    for (size_t i = 0; i < 3; ++i)
    {
        m_vel[i] = Array<OneD, NekDouble>(NumDof, 0.);
    }
    // read prescribed motion DoFs
    m_dirDoFs.clear();
    for (int i = 0; i < NumDof; ++i)
    {
        m_dirDoFs.insert(i);
    }
    mssgTag = pSolver->FirstChildElement("MOTIONPRESCRIBED");
    if (mssgTag)
    {
        std::vector<std::string> values;
        mssgStr = mssgTag->GetText();
        ParseUtils::GenerateVector(mssgStr, values);
        ASSERTL0(values.size() == NumDof,
                 "MOTIONPRESCRIBED vector should be of size " +
                     std::to_string(NumDof));
        for (int i = 0; i < NumDof; ++i)
        {
            if (EvaluateExpression(session, values[i]) == 0)
            {
                m_dirDoFs.erase(i);
            }
        }
    }
    m_hasRotation =
        m_hasRotation || m_dirDoFs.find(m_spacedim) == m_dirDoFs.end();
    m_hasFreeMotion = m_dirDoFs.size() < NumDof;
    // read mass matrix
    m_M     = Array<OneD, NekDouble>(NumDof * NumDof, 0.);
    mssgTag = pSolver->FirstChildElement("MASS");
    ASSERTL0(m_dirDoFs.size() == NumDof || mssgTag, "Mass matrix is required.");
    if (mssgTag)
    {
        std::vector<std::string> values;
        mssgStr = mssgTag->GetText();
        ParseUtils::GenerateVector(mssgStr, values);
        ASSERTL0(values.size() == 1, "Mass matrix should be a scalar.");
        m_mass = EvaluateExpression(session, values[0]);
        for (int i = 0; i < m_spacedim; ++i)
        {
            m_M[i + i * NumDof] = m_mass;
        }
    }
    mssgTag = pSolver->FirstChildElement("ROTATIONINERTIA");
    if (mssgTag)
    {
        std::vector<std::string> values;
        mssgStr = mssgTag->GetText();
        ParseUtils::GenerateVector(mssgStr, values);
        ASSERTL0(values.size() == 1, "Inertia should be a scalar.");
        m_rotaionInertia         = EvaluateExpression(session, values[0]);
        m_M[NumDof * NumDof - 1] = m_rotaionInertia;
    }
    // read damping matrix
    m_C     = Array<OneD, NekDouble>(NumDof * NumDof, 0.);
    mssgTag = pSolver->FirstChildElement("DAMPING");
    if (mssgTag)
    {
        std::vector<std::string> values;
        mssgStr = mssgTag->GetText();
        ParseUtils::GenerateVector(mssgStr, values);
        ASSERTL0(values.size() == NumDof * NumDof,
                 "Damping matrix should be of size " + std::to_string(NumDof) +
                     "X" + std::to_string(NumDof));
        int count = 0;
        for (int i = 0; i < NumDof; ++i)
        {
            for (int j = 0; j < NumDof; ++j)
            {
                m_C[count] = EvaluateExpression(session, values[count]);
                ++count;
            }
        }
    }
    // read rigidity matrix
    m_K     = Array<OneD, NekDouble>(NumDof * NumDof, 0.);
    mssgTag = pSolver->FirstChildElement("RIGIDITY");
    if (mssgTag)
    {
        std::vector<std::string> values;
        mssgStr = mssgTag->GetText();
        ParseUtils::GenerateVector(mssgStr, values);
        ASSERTL0(values.size() == NumDof * NumDof,
                 "Rigidity matrix should be of size " + std::to_string(NumDof) +
                     "X" + std::to_string(NumDof));
        int count = 0;
        for (int i = 0; i < NumDof; ++i)
        {
            for (int j = 0; j < NumDof; ++j)
            {
                m_K[count] = EvaluateExpression(session, values[count]);
                ++count;
            }
        }
    }
    // read pivot point
    mssgTag = pSolver->FirstChildElement("PIVOTPOINT");
    m_pivot = Array<OneD, NekDouble>(m_spacedim, 0.);
    if (mssgTag)
    {
        std::vector<std::string> values;
        mssgStr = mssgTag->GetText();
        ParseUtils::GenerateVector(mssgStr, values);
        ASSERTL0(values.size() == m_spacedim,
                 "PIVOTPOINT vector should be of size " +
                     std::to_string(m_spacedim));
        for (int i = 0; i < m_spacedim; ++i)
        {
            m_pivot[i] = EvaluateExpression(session, values[i]);
        }
    }
    Vmath::Vcopy(m_spacedim, m_pivot, 1, pivot, 1);
    // read the distance between pivotpoint and masscenter
    mssgTag         = pSolver->FirstChildElement("PIVOTDISTANCE");
    m_pivotdistance = 0.;
    if (mssgTag)
    {
        std::vector<std::string> values;
        mssgStr = mssgTag->GetText();
        ParseUtils::GenerateVector(mssgStr, values);
        ASSERTL0(values.size() == 1, "PivotDistance should be a scalar.");
        m_pivotdistance = EvaluateExpression(session, values[0]);
    }
    // read Newmark Beta paramters
    m_timestep = session->GetParameter("TimeStep");
    m_beta     = 0.25;
    m_gamma    = 0.51;
    if (session->DefinesParameter("NewmarkBeta"))
    {
        m_beta = session->GetParameter("NewmarkBeta");
    }
    if (session->DefinesParameter("NewmarkGamma"))
    {
        m_gamma = session->GetParameter("NewmarkGamma");
    }
}

References ASSERTL0, EvaluateExpression(), Nektar::ParseUtils::GenerateVector(), m_beta, m_C, m_dirDoFs, m_gamma, m_hasFreeMotion, m_hasRotation, m_K, m_M, m_mass, m_pivot, m_pivotdistance, m_rotaionInertia, m_spacedim, m_timestep, m_vel, and Vmath::Vcopy().

Referenced by InitObject().

◆ InitObject()

void Nektar::RigidSolver::InitObject	(	const LibUtilities::SessionReaderSharedPtr	session,
		const MultiRegions::ExpListSharedPtr &	pField,
		Array< OneD, NekDouble >	pivot
	)

Definition at line 56 of file RigidSolver.cpp.

{
    bool isH1d, isH2d;
    session->MatchSolverInfo("Homogeneous", "1D", isH1d, false);
    session->MatchSolverInfo("Homogeneous", "2D", isH2d, false);
    bool singleMode, halfMode;
    session->MatchSolverInfo("ModeType", "SingleMode", singleMode, false);
    session->MatchSolverInfo("ModeType", "HalfMode", halfMode, false);
    if (singleMode || halfMode)
    {
        isH1d = false;
    }
    int expdim            = isH2d ? 1 : pField->GetGraph()->GetMeshDimension();
    m_spacedim            = expdim + (isH1d ? 1 : 0) + (isH2d ? 2 : 0);
    m_isRoot              = pField->GetComm()->TreatAsRankZero();
    m_index               = 1;
    m_currentTime         = -1.;
    m_extForceXYZ         = Array<OneD, NekDouble>(6, 0.0);
    m_oldFvis             = Array<OneD, NekDouble>(6, 0.0);
    TiXmlElement *pSolver = session->GetElement("Nektar/RIGIDSOLVER");
    LoadParameters(session, pSolver);
    InitBodySolver(session, pSolver, pivot);
    CheckParameters();
}

References CheckParameters(), InitBodySolver(), LoadParameters(), m_currentTime, m_extForceXYZ, m_index, m_isRoot, m_oldFvis, and m_spacedim.

Referenced by Nektar::PressDecompVCSFSI::v_InitObject(), and Nektar::VCSFSI::v_InitObject().

◆ LoadParameters()

void Nektar::RigidSolver::LoadParameters	(	const LibUtilities::SessionReaderSharedPtr	session,
		const TiXmlElement *	pSolver
	)

protected

Definition at line 205 of file RigidSolver.cpp.

{
    const TiXmlElement *funcNameElmt;
    // load frame velocity
    funcNameElmt = pSolver->FirstChildElement("FRAMEVELOCITY");
    if (funcNameElmt)
    {
        std::string FuncName = funcNameElmt->GetText();
        ASSERTL0(session->DefinesFunction(FuncName),
                 "Function '" + FuncName + "' is not defined in the session.");
        // linear velocity
        for (int i = 0; i < m_spacedim; ++i)
        {
            std::string var = session->GetVariable(i);
            ParserFunctionToMap(false, session, FuncName, var,
                                m_frameVelFunction, i);
        }
        // linear displacement and acceleration
        std::vector<std::string> linearDispVar = {"X", "Y", "Z"};
        std::vector<std::string> linearAcceVar = {"A_x", "A_y", "A_z"};
        for (int i = 0; i < m_spacedim; ++i)
        {
            ParserFunctionToMap(true, session, FuncName, linearDispVar[i],
                                m_frameVelFunction, i + 6);
            ParserFunctionToMap(true, session, FuncName, linearAcceVar[i],
                                m_frameVelFunction, i + 12);
        }
        // angular velocities
        std::vector<std::string> angularVar = {"Omega_x", "Omega_y", "Omega_z"};
        m_hasRotation                       = false;
        for (int i = 0; i < 3; ++i)
        {
            std::string var = angularVar[i];
            ParserFunctionToMap(false, session, FuncName, var,
                                m_frameVelFunction, i + 3);
            m_hasRotation = m_hasRotation || m_frameVelFunction.find(i + 3) !=
                                                 m_frameVelFunction.end();
        }
        // angular displacement and acceleration
        std::vector<std::string> angularDispVar  = {"Theta_x", "Theta_y",
                                                    "Theta_z"};
        std::vector<std::string> angularAccepVar = {"DOmega_x", "DOmega_y",
                                                    "DOmega_z"};
        for (int i = 0; i < 3; ++i)
        {
            ParserFunctionToMap(true, session, FuncName, angularDispVar[i],
                                m_frameVelFunction, i + 3 + 6);
            ParserFunctionToMap(true, session, FuncName, angularAccepVar[i],
                                m_frameVelFunction, i + 3 + 12);
        }
    }
 
    // load external force
    funcNameElmt = pSolver->FirstChildElement("EXTERNALFORCE");
    if (funcNameElmt)
    {
        std::string FuncName = funcNameElmt->GetText();
        if (session->DefinesFunction(FuncName))
        {
            std::vector<std::string> forceVar = {"Fx", "Fy", "Fz"};
            for (int i = 0; i < m_spacedim; ++i)
            {
                std::string var = forceVar[i];
                ParserFunctionToMap(false, session, FuncName, var,
                                    m_extForceFunction, i);
            }
            std::vector<std::string> momentVar = {"Mx", "My", "Mz"};
            for (int i = 0; i < 3; ++i)
            {
                std::string var = momentVar[i];
                ParserFunctionToMap(false, session, FuncName, var,
                                    m_extForceFunction, i + 3);
            }
        }
    }
 
    // OutputFile
    const TiXmlElement *mssgTag = pSolver->FirstChildElement("OutputFile");
    std::string filename;
    if (mssgTag)
    {
        filename = mssgTag->GetText();
    }
    else
    {
        filename = session->GetSessionName();
    }
    if (!(filename.length() >= 4 &&
          filename.substr(filename.length() - 4) == ".mrf"))
    {
        filename += ".mrf";
    }
    if (m_isRoot)
    {
        m_outputStream.open(filename.c_str());
        if (m_spacedim == 2)
        {
            m_outputStream
                << "Variables = t, x, ux, ax, y, uy, ay, theta, omega, domega"
                << std::endl;
        }
        else if (m_spacedim == 3)
        {
            m_outputStream << "Variables = t, x, ux, ax, y, uy, ay, z, uz, az, "
                              "theta, omega, domega"
                           << std::endl;
        }
    }
 
    // output frequency
    m_outputFrequency = 1;
    mssgTag           = pSolver->FirstChildElement("OutputFrequency");
    if (mssgTag)
    {
        std::vector<std::string> values;
        std::string mssgStr = mssgTag->GetText();
        m_outputFrequency   = round(EvaluateExpression(session, mssgStr));
    }
    ASSERTL0(m_outputFrequency > 0,
             "OutputFrequency should be greater than zero.");
}

References ASSERTL0, EvaluateExpression(), m_extForceFunction, m_frameVelFunction, m_hasRotation, m_isRoot, m_outputFrequency, m_outputStream, m_spacedim, and ParserFunctionToMap().

Referenced by InitObject().

◆ ParserFunctionToMap()

void Nektar::RigidSolver::ParserFunctionToMap	(	const bool	allowzero,
		const LibUtilities::SessionReaderSharedPtr	session,
		const std::string &	FuncName,
		const std::string &	var,
		std::map< int, LibUtilities::EquationSharedPtr > &	result,
		const int	i
	)

protected

Definition at line 189 of file RigidSolver.cpp.

{
    if (session->DefinesFunction(FuncName, var))
    {
        LibUtilities::EquationSharedPtr equ =
            session->GetFunction(FuncName, var);
        if (allowzero || "0" != equ->GetExpression())
        {
            result[i] = equ;
        }
    }
}

Referenced by LoadParameters().

◆ SetInitialConditions() [1/2]

void Nektar::RigidSolver::SetInitialConditions	(	const LibUtilities::SessionReaderSharedPtr	session,
		Array< OneD, NekDouble >	MRFData
	)

Definition at line 860 of file RigidSolver.cpp.

{
    NekDouble time = 0.;
    std::map<std::string, std::string> fieldMetaDataMap;
    std::vector<std::string> strFrameData = {
        "X",   "Y",   "Z",   "Theta_x",  "Theta_y",  "Theta_z",
        "U",   "V",   "W",   "Omega_x",  "Omega_y",  "Omega_z",
        "A_x", "A_y", "A_z", "DOmega_x", "DOmega_y", "DOmega_z"};
    std::map<std::string, NekDouble> fileData;
    if (session->DefinesFunction("InitialConditions"))
    {
        for (int i = 0; i < session->GetVariables().size(); ++i)
        {
            if (session->GetFunctionType("InitialConditions",
                                         session->GetVariable(i)) ==
                LibUtilities::eFunctionTypeFile)
            {
                std::string filename = session->GetFunctionFilename(
                    "InitialConditions", session->GetVariable(i));
                fs::path pfilename(filename);
                // redefine path for parallel file which is in directory
                if (fs::is_directory(pfilename))
                {
                    fs::path metafile("Info.xml");
                    fs::path fullpath = pfilename / metafile;
                    filename          = LibUtilities::PortablePath(fullpath);
                }
                LibUtilities::FieldIOSharedPtr fld =
                    LibUtilities::FieldIO::CreateForFile(session, filename);
                fld->ImportFieldMetaData(filename, fieldMetaDataMap);
 
                // check to see if time is defined
                if (fieldMetaDataMap != LibUtilities::NullFieldMetaDataMap)
                {
                    if (fieldMetaDataMap.find("Time") != fieldMetaDataMap.end())
                    {
                        time = std::stod(fieldMetaDataMap["Time"]);
                    }
                    fileData.clear();
                    for (auto &var : strFrameData)
                    {
                        if (fieldMetaDataMap.find(var) !=
                            fieldMetaDataMap.end())
                        {
                            fileData[var] = std::stod(fieldMetaDataMap[var]);
                        }
                    }
                    if (fileData.size() == strFrameData.size())
                    {
                        break;
                    }
                }
            }
        }
    }
    if (session->DefinesCmdLineArgument("set-start-time"))
    {
        time = std::stod(
            session->GetCmdLineArgument<std::string>("set-start-time"));
    }
    if (fileData.size() == strFrameData.size())
    {
        int NumDofm1 = m_vel[0].size() - 1;
        for (int i = 0; i < m_spacedim; ++i)
        {
            m_vel[0][i] = fileData[strFrameData[i]];
            m_vel[1][i] = fileData[strFrameData[i + 6]];
            m_vel[2][i] = fileData[strFrameData[i + 12]];
        }
        m_vel[0][NumDofm1] = fileData[strFrameData[5]];
        m_vel[1][NumDofm1] = fileData[strFrameData[11]];
        m_vel[2][NumDofm1] = fileData[strFrameData[17]];
        if (m_hasRotation)
        {
            m_inertialPosition[0] = m_vel[0][0];
            m_inertialPosition[1] = m_vel[0][1];
        }
    }
    std::map<int, NekDouble> Dirs;
    UpdatePrescribed(time, Dirs);
    SetInitialConditions(Dirs);
    UpdateMRFData(MRFData);
    // output initial status for rigid body
    if (m_isRoot)
    {
        Array<OneD, Array<OneD, NekDouble>> tmpVel;
        if (m_hasRotation)
        {
            tmpVel = Array<OneD, Array<OneD, NekDouble>>(m_vel.size());
            for (size_t i = 0; i < m_vel.size(); ++i)
            {
                tmpVel[i] = Array<OneD, NekDouble>(m_vel[i].size());
                Vmath::Vcopy(m_vel[i].size(), m_vel[i], 1, tmpVel[i], 1);
            }
            // transform u, du to body frame
            Array<OneD, NekDouble> angle(3, 0.);
            angle[2] = tmpVel[0][m_spacedim];
            m_frame.SetAngle(angle);
            m_frame.BodyToInerital(m_spacedim, tmpVel[1], tmpVel[1]);
            m_frame.BodyToInerital(m_spacedim, tmpVel[2], tmpVel[2]);
            tmpVel[2][0] -= tmpVel[1][m_spacedim] * tmpVel[1][1];
            tmpVel[2][1] += tmpVel[1][m_spacedim] * tmpVel[1][0];
            tmpVel[0][0] = m_inertialPosition[0];
            tmpVel[0][1] = m_inertialPosition[1];
        }
        else
        {
            tmpVel = m_vel;
        }
        m_outputStream << boost::format("%25.19e") % time << " ";
        for (size_t i = 0; i < tmpVel[0].size(); ++i)
        {
            m_outputStream << boost::format("%25.19e") % tmpVel[0][i] << " "
                           << boost::format("%25.19e") % tmpVel[1][i] << " "
                           << boost::format("%25.19e") % tmpVel[2][i] << " ";
        }
        m_outputStream << std::endl;
    }
}

References Nektar::FrameTransform::BodyToInerital(), Nektar::LibUtilities::FieldIO::CreateForFile(), Nektar::LibUtilities::eFunctionTypeFile, m_frame, m_hasRotation, m_inertialPosition, m_isRoot, m_outputStream, m_spacedim, m_vel, Nektar::LibUtilities::NullFieldMetaDataMap, Nektar::LibUtilities::PortablePath(), Nektar::FrameTransform::SetAngle(), SetInitialConditions(), UpdateMRFData(), UpdatePrescribed(), and Vmath::Vcopy().

Referenced by SetInitialConditions(), Nektar::PressDecompVCSFSI::v_DoInitialise(), and Nektar::VCSFSI::v_DoInitialise().

◆ SetInitialConditions() [2/2]

void Nektar::RigidSolver::SetInitialConditions ( std::map< int, NekDouble > & Dirs )

protected

Definition at line 1000 of file RigidSolver.cpp.

{
    for (auto it : Dirs)
    {
        int NumDof  = m_vel[0].size();
        int NumDof2 = NumDof << 1;
        if (it.first < m_vel[0].size())
        {
            m_vel[1][it.first] = it.second;
        }
        else if (it.first < NumDof2)
        {
            m_vel[0][it.first - NumDof] = it.second;
        }
        else
        {
            m_vel[2][it.first - NumDof2] = it.second;
        }
    }
}

References m_vel.

◆ SetMovableDoFs()

void Nektar::RigidSolver::SetMovableDoFs ( std::vector< bool > & moveDoFs )

Definition at line 141 of file RigidSolver.cpp.

{
    for (const auto &it : m_frameVelFunction)
    {
        if (it.first < 6)
        {
            moveDoFs[it.first] = true;
        }
    }
    for (int i = 0; i < m_spacedim; ++i)
    {
        if (m_dirDoFs.find(i) == m_dirDoFs.end())
        {
            moveDoFs[i] = true;
        }
    }
    if (m_dirDoFs.find(m_spacedim) == m_dirDoFs.end())
    {
        moveDoFs[5] = true;
    }
    if (moveDoFs[5])
    {
        // since m_frameVelFunction is defined in inerital frame
        // motions in 2 directions occur if there is rotation
        if (moveDoFs[0] || moveDoFs[1])
        {
            moveDoFs[0] = true;
            moveDoFs[1] = true;
        }
    }
}

References m_dirDoFs, m_frameVelFunction, and m_spacedim.

Referenced by Nektar::PressDecompVCSFSI::v_InitObject(), and Nektar::VCSFSI::v_InitObject().

◆ SetNewmarkBetaSolver()

void Nektar::RigidSolver::SetNewmarkBetaSolver ( Array< OneD, NekDouble > & AddedMass )

Definition at line 849 of file RigidSolver.cpp.

{
    int NumDof = m_spacedim + 1;
    if (AddedMass.size() >= NumDof * NumDof)
    {
        Vmath::Vadd(NumDof * NumDof, AddedMass, 1, m_M, 1, m_M, 1);
    }
    m_bodySolver.SetNewmarkBeta(m_beta, m_gamma, m_timestep, m_M, m_C, m_K,
                                m_dirDoFs, m_solveType);
}

References m_beta, m_bodySolver, m_C, m_dirDoFs, m_gamma, m_K, m_M, m_solveType, m_spacedim, m_timestep, Nektar::Newmark_BetaSolver::SetNewmarkBeta(), and Vmath::Vadd().

Referenced by Nektar::PressDecompVCSFSI::v_DoInitialise(), and Nektar::VCSFSI::v_DoInitialise().

◆ SetOldFvis()

void Nektar::RigidSolver::SetOldFvis ( Array< OneD, NekDouble > force )

Definition at line 522 of file RigidSolver.cpp.

{
    for (int i = 0; i < 6; ++i)
    {
        m_oldFvis[i] = force[6 + i];
    }
}

References m_oldFvis.

Referenced by Nektar::PressDecompVCSFSI::v_DoInitialise(), and Nektar::VCSFSI::v_DoInitialise().

◆ SolveBodyFrame()

void Nektar::RigidSolver::SolveBodyFrame	(	Array< OneD, Array< OneD, NekDouble > > &	bodyVel,
		const Array< OneD, NekDouble > &	forcebody,
		std::map< int, NekDouble > &	Dirs
	)

private

Definition at line 716 of file RigidSolver.cpp.

{
    m_bodySolver.SolvePrescribed(bodyVel, Dirs); // at most 1 rotation
    Array<OneD, NekDouble> force(6, 0.), angle(3, 0.);
    if (!m_hasFreeMotion)
    {
        angle[2] = bodyVel[0][m_spacedim];
        m_frame.SetAngle(angle);
        m_frame.IneritalToBody(m_spacedim, bodyVel[1], bodyVel[1]);
        m_frame.IneritalToBody(m_spacedim, bodyVel[2], bodyVel[2]);
        bodyVel[2][0] += bodyVel[1][m_spacedim] * bodyVel[1][1];
        bodyVel[2][1] -= bodyVel[1][m_spacedim] * bodyVel[1][0];
        return;
    }
    if (m_dirDoFs.find(m_spacedim) != m_dirDoFs.end())
    {
        // known rotation
        angle[2] = bodyVel[0][m_spacedim];
        m_frame.SetAngle(angle);
        m_frame.IneritalToBody(m_spacedim, m_extForceXYZ, force);
        for (int i = 0; i < m_spacedim; ++i)
        {
            force[i] = forcebody[i] + force[i];
        }
        force[0] += m_mass * bodyVel[1][m_spacedim] * bodyVel[1][m_spacedim] *
                    m_pivotdistance;
        force[1] -= m_mass * bodyVel[2][m_spacedim] * m_pivotdistance;
        m_bodySolver.SolveFreeVarMat(bodyVel, force, m_mass);
    }
    else
    {
        // all free
        Array<OneD, Array<OneD, NekDouble>> tmpbodyVel(bodyVel.size());
        for (size_t i = 0; i < bodyVel.size(); ++i)
        {
            tmpbodyVel[i] = Array<OneD, NekDouble>(bodyVel[i].size());
            Vmath::Vcopy(bodyVel[i].size(), bodyVel[i], 1, tmpbodyVel[i], 1);
        }
        for (int iter = 0; iter < 2; ++iter)
        {
            if (iter > 0)
            {
                for (size_t i = 0; i < bodyVel.size(); ++i)
                {
                    Vmath::Vcopy(bodyVel[i].size() - 1, bodyVel[i], 1,
                                 tmpbodyVel[i], 1);
                }
            }
            angle[2] = tmpbodyVel[0][m_spacedim];
            m_frame.SetAngle(angle);
            m_frame.IneritalToBody(3, m_extForceXYZ, force);
            for (int i = 0; i < m_spacedim; ++i)
            {
                force[i] = forcebody[i] + force[i];
            }
            force[m_spacedim] = forcebody[5] + m_extForceXYZ[5];
            m_bodySolver.SolveFreeVarMat(tmpbodyVel, force, m_mass);
        }
        // copy final results
        for (size_t i = 0; i < bodyVel.size(); ++i)
        {
            Vmath::Vcopy(bodyVel[i].size(), tmpbodyVel[i], 1, bodyVel[i], 1);
        }
    }
}

References Nektar::FrameTransform::IneritalToBody(), m_bodySolver, m_dirDoFs, m_extForceXYZ, m_frame, m_hasFreeMotion, m_mass, m_pivotdistance, m_spacedim, Nektar::FrameTransform::SetAngle(), Nektar::Newmark_BetaSolver::SolveFreeVarMat(), Nektar::Newmark_BetaSolver::SolvePrescribed(), and Vmath::Vcopy().

Referenced by SolveBodyMotion().

◆ SolveBodyMotion()

void Nektar::RigidSolver::SolveBodyMotion	(	Array< OneD, Array< OneD, NekDouble > > &	bodyVel,
		const Array< OneD, NekDouble > &	forcebody,
		std::map< int, NekDouble > &	Dirs
	)

protected

Definition at line 647 of file RigidSolver.cpp.

{
    if (0 == m_solveType)
    {
        SolveInertialFrame(bodyVel, forcebody, Dirs);
    }
    else if (1 == m_solveType)
    {
        SolveRotOneFree(bodyVel, forcebody, Dirs);
    }
    else if (2 == m_solveType)
    {
        SolveBodyFrame(bodyVel, forcebody, Dirs);
    }
    else
    {
        ASSERTL0(false, "Unsupported rigid solver type.");
    }
}

References ASSERTL0, m_solveType, SolveBodyFrame(), SolveInertialFrame(), and SolveRotOneFree().

Referenced by UpdateFrameVelocity().

◆ SolveInertialFrame()

void Nektar::RigidSolver::SolveInertialFrame	(	Array< OneD, Array< OneD, NekDouble > > &	bodyVel,
		const Array< OneD, NekDouble > &	forcebody,
		std::map< int, NekDouble > &	Dirs
	)

private

Definition at line 669 of file RigidSolver.cpp.

{
    // only translational motion or prescribed tranlation with rotational
    // solve in absolute frame and transform to body frame
    m_bodySolver.SolvePrescribed(bodyVel, Dirs);
    if (m_hasFreeMotion)
    {
        Array<OneD, NekDouble> force(6, 0.);
        for (int i = 0; i < m_spacedim; ++i)
        {
            force[i] = forcebody[i] + m_extForceXYZ[i];
        }
        force[m_spacedim] = forcebody[5] + m_extForceXYZ[5];
        m_bodySolver.SolveFreeFixMat(bodyVel, force);
    }
}

References m_bodySolver, m_extForceXYZ, m_hasFreeMotion, m_spacedim, Nektar::Newmark_BetaSolver::SolveFreeFixMat(), and Nektar::Newmark_BetaSolver::SolvePrescribed().

Referenced by SolveBodyMotion().

◆ SolveRotOneFree()

void Nektar::RigidSolver::SolveRotOneFree	(	Array< OneD, Array< OneD, NekDouble > > &	bodyVel,
		const Array< OneD, NekDouble > &	forcebody,
		std::map< int, NekDouble > &	Dirs
	)

private

Definition at line 689 of file RigidSolver.cpp.

{
    // one direction free
    NekDouble uy = 0.;
    if (Dirs.find(1) != Dirs.end())
    {
        uy = Dirs[1];
        Dirs.erase(1);
    }
    m_bodySolver.SolvePrescribed(bodyVel, Dirs);
    Array<OneD, NekDouble> force(6, 0.), angle(3, 0.);
    angle[2] = bodyVel[0][m_spacedim];
    m_frame.SetAngle(angle);
    m_frame.IneritalToBody(3, m_extForceXYZ, force);
    for (int i = 0; i < m_spacedim; ++i)
    {
        force[i] = forcebody[i] + force[i];
    }
    force[0] += m_mass * bodyVel[1][m_spacedim] * bodyVel[1][m_spacedim] *
                m_pivotdistance;
    force[1] -= m_mass * bodyVel[2][m_spacedim] * m_pivotdistance;
    m_bodySolver.SolveOneFree(bodyVel, force, angle, uy, m_mass);
}

References Nektar::FrameTransform::IneritalToBody(), m_bodySolver, m_extForceXYZ, m_frame, m_mass, m_pivotdistance, m_spacedim, Nektar::FrameTransform::SetAngle(), Nektar::Newmark_BetaSolver::SolveOneFree(), and Nektar::Newmark_BetaSolver::SolvePrescribed().

Referenced by SolveBodyMotion().

◆ UpdateFrameVelocity()

void Nektar::RigidSolver::UpdateFrameVelocity	(	Array< OneD, NekDouble > &	aeroforce,
		const NekDouble &	time,
		Array< OneD, NekDouble > &	MRFData
	)

Updates the forcing array with the current required forcing.

Parameters

pFields
time

Definition at line 535 of file RigidSolver.cpp.

{
    if (m_currentTime >= time)
    {
        return;
    }
 
    for (int i = 0; i < 6; ++i)
    {
        aeroforce[i] +=
            2. * aeroforce[6 + i] - m_oldFvis[i]; // extrapolate viscous
        m_oldFvis[i] = aeroforce[6 + i];
    }
    m_currentTime = time;
    std::map<int, NekDouble> Dirs;
    UpdatePrescribed(time, Dirs);
    // compute the velocites whoes functions are provided in inertial frame
    if (m_hasFreeMotion)
    {
        for (auto it : m_extForceFunction)
        {
            m_extForceXYZ[it.first] = it.second->Evaluate(0., 0., 0., time);
        }
    }
    SolveBodyMotion(m_vel, aeroforce, Dirs);
    Array<OneD, Array<OneD, NekDouble>> tmpVel;
    if (m_hasRotation)
    {
        // update inertial position of the body
        tmpVel = Array<OneD, Array<OneD, NekDouble>>(m_vel.size());
        for (size_t i = 0; i < m_vel.size(); ++i)
        {
            tmpVel[i] = Array<OneD, NekDouble>(m_vel[i].size());
            Vmath::Vcopy(m_vel[i].size(), m_vel[i], 1, tmpVel[i], 1);
        }
        Array<OneD, NekDouble> angle(3, 0.);
        angle[2] = tmpVel[0][m_spacedim];
        m_frame.SetAngle(angle);
        m_frame.BodyToInerital(m_spacedim, tmpVel[1], tmpVel[1]);
        m_frame.BodyToInerital(m_spacedim, tmpVel[2], tmpVel[2]);
        tmpVel[2][0] -= tmpVel[1][m_spacedim] * tmpVel[1][1];
        tmpVel[2][1] += tmpVel[1][m_spacedim] * tmpVel[1][0];
        m_inertialPosition[0] +=
            m_timestep * (tmpVel[1][0] + 0.5 * m_timestep * tmpVel[2][0]);
        m_inertialPosition[1] +=
            m_timestep * (tmpVel[1][1] + 0.5 * m_timestep * tmpVel[2][1]);
        tmpVel[0][0] = m_inertialPosition[0];
        tmpVel[0][1] = m_inertialPosition[1];
    }
    if (m_isRoot && m_index % m_outputFrequency == 0)
    {
        if (!m_hasRotation)
        {
            tmpVel = m_vel;
        }
        m_outputStream << boost::format("%25.19e") % time << " ";
        for (size_t i = 0; i < tmpVel[0].size(); ++i)
        {
            m_outputStream << boost::format("%25.19e") % tmpVel[0][i] << " "
                           << boost::format("%25.19e") % tmpVel[1][i] << " "
                           << boost::format("%25.19e") % tmpVel[2][i] << " ";
        }
        m_outputStream << std::endl;
    }
    // set displacements at the current time step
    UpdateMRFData(MRFData);
    ++m_index;
}

References Nektar::FrameTransform::BodyToInerital(), m_currentTime, m_extForceFunction, m_extForceXYZ, m_frame, m_hasFreeMotion, m_hasRotation, m_index, m_inertialPosition, m_isRoot, m_oldFvis, m_outputFrequency, m_outputStream, m_spacedim, m_timestep, m_vel, Nektar::FrameTransform::SetAngle(), SolveBodyMotion(), UpdateMRFData(), UpdatePrescribed(), and Vmath::Vcopy().

Referenced by Nektar::VCSFSI::v_SolveSolid().

◆ UpdateMRFData()

void Nektar::RigidSolver::UpdateMRFData ( Array< OneD, NekDouble > & MRFData )

protected

MRFData: X, Y, Z, Theta_x, Theta_y, Theta_z, [inertial frame 0-5] U, V, W, Omega_x, Omega_y, Omega_z, [body frame 6-11] A_x, A_y, A_z, DOmega_x, DOmega_y, DOmega_z, [body frame 12-17] pivot_x, pivot_y, pivot_z, [body frame]

Definition at line 606 of file RigidSolver.cpp.

{
    /// MRFData:
    /// X, Y, Z, Theta_x, Theta_y, Theta_z, [inertial frame 0-5]
    /// U, V, W, Omega_x, Omega_y, Omega_z, [body frame 6-11]
    /// A_x, A_y, A_z, DOmega_x, DOmega_y, DOmega_z, [body frame 12-17]
    /// pivot_x, pivot_y, pivot_z, [body frame]
    for (int i = 0; i < 18; ++i)
    {
        MRFData[i] = 0.0;
    }
 
    if (m_hasRotation)
    {
        MRFData[0] = m_inertialPosition[0];
        MRFData[1] = m_inertialPosition[1];
    }
    else
    {
        MRFData[0] = m_vel[0][0];
        MRFData[1] = m_vel[0][1];
    }
 
    if (m_hasRotation)
    {
        MRFData[5] = m_vel[0][m_spacedim];
    }
 
    for (int i = 0; i < m_spacedim; ++i)
    {
        MRFData[i + 6]  = m_vel[1][i];
        MRFData[i + 12] = m_vel[2][i];
    }
 
    if (m_hasRotation)
    {
        MRFData[11] = m_vel[1][m_spacedim];
        MRFData[17] = m_vel[2][m_spacedim];
    }
}

References m_hasRotation, m_inertialPosition, m_spacedim, and m_vel.

Referenced by SetInitialConditions(), and UpdateFrameVelocity().

◆ UpdatePrescribed()

void Nektar::RigidSolver::UpdatePrescribed	(	const NekDouble &	time,
		std::map< int, NekDouble > &	Dirs
	)

protected

Definition at line 478 of file RigidSolver.cpp.

{
    int NumDof = m_spacedim + 1;
    for (auto it : m_frameVelFunction)
    {
        if (it.first < 3)
        {
            Dirs[it.first] = it.second->Evaluate(0., 0., 0., time);
        }
        else if (it.first == 5)
        {
            Dirs[m_spacedim] = it.second->Evaluate(0., 0., 0., time);
        }
        else if (it.first < 9)
        {
            Dirs[NumDof + it.first - 6] = it.second->Evaluate(0., 0., 0., time);
        }
        else if (it.first == 11)
        {
            Dirs[NumDof + m_spacedim] = it.second->Evaluate(0., 0., 0., time);
        }
        else if (it.first < 15)
        {
            Dirs[(NumDof << 1) + it.first - 12] =
                it.second->Evaluate(0., 0., 0., time);
        }
        else if (it.first == 17)
        {
            Dirs[(NumDof << 1) + m_spacedim] =
                it.second->Evaluate(0., 0., 0., time);
        }
    }
    for (auto i : m_dirDoFs)
    {
        if (Dirs.find(i) == Dirs.end())
        {
            Dirs[i]                 = 0.;
            Dirs[i + NumDof]        = 0.;
            Dirs[i + (NumDof << 1)] = 0.;
        }
    }
}