#include <RigidSolver.h>

Public Member Functions
	Newmark_BetaSolver ()

	~Newmark_BetaSolver ()

void	SetNewmarkBeta (NekDouble beta, NekDouble gamma, NekDouble dt, Array< OneD, NekDouble > M, Array< OneD, NekDouble > C, Array< OneD, NekDouble > K, std::set< int > DirDoFs, int solveType)

void	SolvePrescribed (Array< OneD, Array< OneD, NekDouble > > u, std::map< int, NekDouble > motionPrescribed)

void	SolveFreeFixMat (Array< OneD, Array< OneD, NekDouble > > u, Array< OneD, NekDouble > force)

void	SolveFreeVarMat (Array< OneD, Array< OneD, NekDouble > > u, Array< OneD, NekDouble > force, const NekDouble mass)

void	SolveOneFree (Array< OneD, Array< OneD, NekDouble > > u, Array< OneD, NekDouble > force, const Array< OneD, NekDouble > theta, const NekDouble uy, const NekDouble mass)

Public Attributes
int	m_rows

int	m_motionDofs

std::vector< int >	m_index

Array< OneD, NekDouble >	m_coeffs

Array< OneD, Array< OneD, NekDouble > >	m_Matrix

Array< OneD, Array< OneD, NekDouble > >	m_M

Array< OneD, Array< OneD, NekDouble > >	m_C

Array< OneD, Array< OneD, NekDouble > >	m_K

Detailed Description

Definition at line 62 of file RigidSolver.h.

Constructor & Destructor Documentation

◆ Newmark_BetaSolver()

Nektar::Newmark_BetaSolver::Newmark_BetaSolver ( )

inline

Definition at line 65 of file RigidSolver.h.

65{};

◆ ~Newmark_BetaSolver()

Nektar::Newmark_BetaSolver::~Newmark_BetaSolver ( )

inline

Definition at line 66 of file RigidSolver.h.

66{};

Member Function Documentation

◆ SetNewmarkBeta()

void Nektar::Newmark_BetaSolver::SetNewmarkBeta	(	NekDouble	beta,
		NekDouble	gamma,
		NekDouble	dt,
		Array< OneD, NekDouble >	M,
		Array< OneD, NekDouble >	C,
		Array< OneD, NekDouble >	K,
		std::set< int >	DirDoFs,
		int	solveType
	)

Definition at line 1021 of file RigidSolver.cpp.

{
    m_coeffs    = Array<OneD, NekDouble>(5, 0.);
    m_coeffs[0] = 1. / (gamma * dt);
    m_coeffs[1] = 1. / gamma - 1.;
    m_coeffs[2] = beta * dt / gamma;
    m_coeffs[3] = dt * (1. - beta / gamma);
    m_coeffs[4] = (0.5 - beta / gamma) * dt * dt;
 
    m_rows = sqrt(M.size());
    m_index.resize(m_rows, -1);
    m_motionDofs = 0;
    for (int i = 0; i < m_rows; ++i)
    {
        if (DirDoFs.find(i) == DirDoFs.end())
        {
            m_index[m_motionDofs++] = i;
        }
    }
    for (int i = 0, count = m_motionDofs; i < m_rows; ++i)
    {
        if (DirDoFs.find(i) != DirDoFs.end())
        {
            m_index[count++] = i;
        }
    }
    if (1 == solveType)
    {
        ASSERTL0(m_motionDofs == 2,
                 "2 Dofs if body is free only in x direction.");
    }
    if (m_motionDofs)
    {
        Array<OneD, NekDouble> temp;
        m_M = Array<OneD, Array<OneD, NekDouble>>(m_motionDofs);
        m_C = Array<OneD, Array<OneD, NekDouble>>(m_motionDofs);
        m_K = Array<OneD, Array<OneD, NekDouble>>(m_motionDofs);
        DNekMatSharedPtr inverseMatrix =
            MemoryManager<DNekMat>::AllocateSharedPtr(m_motionDofs,
                                                      m_motionDofs, 0.0, eFULL);
        for (int i = 0; i < m_motionDofs; ++i)
        {
            m_M[i]     = Array<OneD, NekDouble>(m_rows, 0.);
            m_C[i]     = Array<OneD, NekDouble>(m_rows, 0.);
            m_K[i]     = Array<OneD, NekDouble>(m_rows, 0.);
            int offset = m_index[i] * m_rows;
            for (int j = 0; j < m_rows; ++j)
            {
                int ind   = offset + m_index[j];
                m_M[i][j] = M[ind];
                m_C[i][j] = C[ind];
                m_K[i][j] = K[ind];
                NekDouble value =
                    m_coeffs[0] * M[ind] + C[ind] + m_coeffs[2] * K[ind];
                if (j < m_motionDofs)
                {
                    inverseMatrix->SetValue(i, j, value);
                }
            }
        }
 
        m_Matrix = Array<OneD, Array<OneD, NekDouble>>(m_motionDofs);
        if (solveType == 0)
        {
            inverseMatrix->Invert();
        }
        for (int i = 0; i < m_motionDofs; ++i)
        {
            m_Matrix[i] = Array<OneD, NekDouble>(m_motionDofs, 0.);
            for (int j = 0; j < m_rows; ++j)
            {
                if (j < m_motionDofs)
                {
                    m_Matrix[i][j] = inverseMatrix->GetValue(i, j);
                }
            }
        }
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::LibUtilities::beta, Nektar::eFULL, m_C, m_coeffs, m_index, m_K, m_M, m_Matrix, m_motionDofs, m_rows, and tinysimd::sqrt().

Referenced by Nektar::RigidSolver::SetNewmarkBetaSolver().

◆ SolveFreeFixMat()

void Nektar::Newmark_BetaSolver::SolveFreeFixMat	(	Array< OneD, Array< OneD, NekDouble > >	u,
		Array< OneD, NekDouble >	force
	)

Definition at line 1146 of file RigidSolver.cpp.

{
    if (m_motionDofs)
    {
        Array<OneD, NekDouble> bm(m_motionDofs, 0.);
        Array<OneD, NekDouble> bk(m_motionDofs, 0.);
        for (int j = 0; j < m_motionDofs; ++j)
        {
            int j1 = m_index[j];
            bm[j]  = m_coeffs[0] * u[1][j1] + m_coeffs[1] * u[2][j1];
            bk[j]  = u[0][j1] + m_coeffs[3] * u[1][j1] + m_coeffs[4] * u[2][j1];
        }
        Array<OneD, NekDouble> rhs(m_motionDofs, 0.);
        for (int i = 0; i < m_motionDofs; ++i)
        {
            rhs[i] = force[m_index[i]];
            for (int j = 0; j < m_motionDofs; ++j)
            {
                rhs[i] += m_M[i][j] * bm[j] - m_K[i][j] * bk[j];
            }
            for (int j = m_motionDofs; j < m_rows; ++j)
            {
                int j1 = m_index[j];
                rhs[i] -= m_M[i][j] * u[2][j1] + m_C[i][j] * u[1][j1] +
                          m_K[i][j] * u[0][j1];
            }
        }
        for (int j = 0; j < m_motionDofs; ++j)
        {
            int j1   = m_index[j];
            u[1][j1] = Vmath::Dot(m_motionDofs, m_Matrix[j], 1, rhs, 1);
            u[0][j1] = m_coeffs[2] * u[1][j1] + bk[j];
            u[2][j1] = m_coeffs[0] * u[1][j1] - bm[j];
        }
    }
}

References Vmath::Dot(), m_C, m_coeffs, m_index, m_K, m_M, m_Matrix, m_motionDofs, and m_rows.

Referenced by Nektar::RigidSolver::SolveInertialFrame().

◆ SolveFreeVarMat()

void Nektar::Newmark_BetaSolver::SolveFreeVarMat	(	Array< OneD, Array< OneD, NekDouble > >	u,
		Array< OneD, NekDouble >	force,
		const NekDouble	mass
	)

Definition at line 784 of file RigidSolver.cpp.

{
    Array<OneD, NekDouble> bm(m_motionDofs, 0.);
    Array<OneD, NekDouble> bk(m_motionDofs, 0.);
    double *dMatrix = new double[m_motionDofs * m_motionDofs];
    int *ipiv       = new int[m_motionDofs];
    double *drhs    = new double[m_motionDofs];
    int info;
 
    for (int j = 0; j < m_motionDofs; ++j)
    {
        int j1 = m_index[j];
        bm[j]  = m_coeffs[0] * u[1][j1] + m_coeffs[1] * u[2][j1];
        bk[j]  = u[0][j1] + m_coeffs[3] * u[1][j1] + m_coeffs[4] * u[2][j1];
    }
    Array<OneD, NekDouble> rhs(m_motionDofs, 0.);
    for (int i = 0; i < m_motionDofs; ++i)
    {
        rhs[i] = force[m_index[i]];
        for (int j = 0; j < m_motionDofs; ++j)
        {
            rhs[i] += m_M[i][j] * bm[j] - m_K[i][j] * bk[j];
        }
        for (int j = m_motionDofs; j < m_rows; ++j)
        {
            int j1 = m_index[j];
            rhs[i] -= m_M[i][j] * u[2][j1] + m_C[i][j] * u[1][j1] +
                      m_K[i][j] * u[0][j1];
        }
    }
    for (int i = 0; i < m_motionDofs; ++i)
    {
        drhs[i] = rhs[i];
        for (int j = 0; j < m_motionDofs; ++j)
        {
            dMatrix[j * m_motionDofs + i] = m_Matrix[i][j];
        }
    }
    dMatrix[1] += mass * u[1][m_rows - 1];
    dMatrix[m_motionDofs] += -mass * u[1][m_rows - 1];
 
    Lapack::DoSgetrf(m_motionDofs, m_motionDofs, dMatrix, m_motionDofs, ipiv,
                     info);
    Lapack::Dgetrs('N', m_motionDofs, 1, dMatrix, m_motionDofs, ipiv, drhs,
                   m_motionDofs, info);
    for (int j = 0; j < m_motionDofs; ++j)
    {
        int j1   = m_index[j];
        u[1][j1] = drhs[j];
    }
 
    for (int j = 0; j < m_motionDofs; ++j)
    {
        int j1   = m_index[j];
        u[0][j1] = m_coeffs[2] * u[1][j1] + bk[j];
        u[2][j1] = m_coeffs[0] * u[1][j1] - bm[j];
    }
 
    delete[] dMatrix;
    delete[] drhs;
    delete[] ipiv;
}

References Lapack::Dgetrs(), Lapack::DoSgetrf(), m_C, m_coeffs, m_index, m_K, m_M, m_Matrix, m_motionDofs, and m_rows.

Referenced by Nektar::RigidSolver::SolveBodyFrame().

◆ SolveOneFree()

void Nektar::Newmark_BetaSolver::SolveOneFree	(	Array< OneD, Array< OneD, NekDouble > >	u,
		Array< OneD, NekDouble >	force,
		const Array< OneD, NekDouble >	theta,
		const NekDouble	uy,
		const NekDouble	mass
	)

e_x \cdot M [du0, du1, du2] + m (e_x \times Omega) \cdot u = F \cdot e_x e_x = (c, -s, 0) (c M00 - s M10, c M01 - s M11, c M02 - s M 12) [du0, du1, du2]^T + (-s Omega, -c Omega, 0) [u0, u1, u2]^T = c F0 - s F1 equation 2, e_y \cdot (u0, u1) = uy = s * u0 + c * u1

Definition at line 1192 of file RigidSolver.cpp.

{
    int iOmega  = m_rows - 1;
    NekDouble c = cos(theta[2]), s = sin(theta[2]);
    NekDouble C00, C01, C10, C11, M0, M1, M2, F0, F1;
    M0  = c * m_M[0][0] - s * m_M[1][0];
    M1  = c * m_M[0][1] - s * m_M[1][1];
    M2  = c * m_M[0][iOmega] - s * m_M[1][iOmega];
    C00 = m_coeffs[0] * M0 - s * u[1][iOmega] * mass;
    C01 = m_coeffs[0] * M1 - c * u[1][iOmega] * mass;
    C10 = s;
    C11 = c;
    F1  = uy;
    // F0
    Array<OneD, NekDouble> bm(m_motionDofs, 0.);
    Array<OneD, NekDouble> bk(m_motionDofs, 0.);
    for (int j = 0; j < m_motionDofs; ++j)
    {
        bm[j] = m_coeffs[0] * u[1][j] + m_coeffs[1] * u[2][j];
        bk[j] = u[0][j] + m_coeffs[3] * u[1][j] + m_coeffs[4] * u[2][j];
    }
    F0 = c * force[0] - s * force[1];
    F0 = F0 + M0 * bm[0] + M1 * bm[1] - M2 * u[2][iOmega];
    // solve
    NekDouble det = 1. / (C00 * C11 - C01 * C10);
    u[1][0]       = det * (C11 * F0 - C01 * F1);
    u[1][1]       = det * (-C10 * F0 + C00 * F1);
    for (int j = 0; j < m_motionDofs; ++j)
    {
        u[0][j] = m_coeffs[2] * u[1][j] + bk[j];
        u[2][j] = m_coeffs[0] * u[1][j] - bm[j];
    }
}

References m_coeffs, m_M, m_motionDofs, and m_rows.

Referenced by Nektar::RigidSolver::SolveRotOneFree().

◆ SolvePrescribed()

void Nektar::Newmark_BetaSolver::SolvePrescribed	(	Array< OneD, Array< OneD, NekDouble > >	u,
		std::map< int, NekDouble >	motionPrescribed
	)

Definition at line 1105 of file RigidSolver.cpp.

{
    for (int i = 0; i < m_rows; ++i)
    {
        if (motionPrescribed.find(i) != motionPrescribed.end())
        {
            int i0       = i + m_rows;
            int i2       = i0 + m_rows;
            NekDouble bm = 0., bk = 0.;
            if (motionPrescribed.find(i2) == motionPrescribed.end())
            {
                bm = m_coeffs[0] * u[1][i] + m_coeffs[1] * u[2][i];
            }
            if (motionPrescribed.find(i0) == motionPrescribed.end())
            {
                bk = u[0][i] + m_coeffs[3] * u[1][i] + m_coeffs[4] * u[2][i];
            }
 
            u[1][i] = motionPrescribed[i];
            if (motionPrescribed.find(i2) == motionPrescribed.end())
            {
                u[2][i] = m_coeffs[0] * u[1][i] - bm;
            }
            else
            {
                u[2][i] = motionPrescribed[i2];
            }
            if (motionPrescribed.find(i0) == motionPrescribed.end())
            {
                u[0][i] = m_coeffs[2] * u[1][i] + bk;
            }
            else
            {
                u[0][i] = motionPrescribed[i0];
            }
        }
    }
}