Nektar::BetaPressureArea Class Reference

#include <BetaPressureArea.h>

Inheritance diagram for Nektar::BetaPressureArea:

Public Member Functions
	BetaPressureArea (Array< OneD, MultiRegions::ExpListSharedPtr > pVessel, const LibUtilities::SessionReaderSharedPtr pSession)
virtual	~BetaPressureArea ()
Public Member Functions inherited from Nektar::PulseWavePressureArea
	PulseWavePressureArea (Array< OneD, MultiRegions::ExpListSharedPtr > &pVessel, const LibUtilities::SessionReaderSharedPtr &pSession)
virtual	~PulseWavePressureArea ()
void	GetPressure (NekDouble &P, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &dAUdx, const NekDouble &gamma=0, const NekDouble &alpha=0.5)
void	GetC (NekDouble &c, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
void	GetW1 (NekDouble &W1, const NekDouble &u, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
void	GetW2 (NekDouble &W2, const NekDouble &u, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
void	GetAFromChars (NekDouble &A, const NekDouble &W1, const NekDouble &W2, const NekDouble &beta, const NekDouble &A0, const NekDouble &alpha=0.5)
void	GetUFromChars (NekDouble &u, const NekDouble &W1, const NekDouble &W2)
void	GetCharIntegral (NekDouble &I, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
void	GetJacobianInverse (NekMatrix< NekDouble > &invJ, const Array< OneD, NekDouble > &Au, const Array< OneD, NekDouble > &uu, const Array< OneD, NekDouble > &beta, const Array< OneD, NekDouble > &A0, const Array< OneD, NekDouble > &alpha, const std::string &type)

Static Public Member Functions
static PulseWavePressureAreaSharedPtr	create (Array< OneD, MultiRegions::ExpListSharedPtr > &pVessel, const LibUtilities::SessionReaderSharedPtr &pSession)

Static Public Attributes
static std::string	className

Protected Member Functions
virtual void	v_GetPressure (NekDouble &P, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &dAUdx, const NekDouble &gamma=0, const NekDouble &alpha=0.5)
virtual void	v_GetC (NekDouble &c, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
virtual void	v_GetW1 (NekDouble &W1, const NekDouble &u, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
virtual void	v_GetW2 (NekDouble &W2, const NekDouble &u, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
virtual void	v_GetAFromChars (NekDouble &A, const NekDouble &W1, const NekDouble &W2, const NekDouble &beta, const NekDouble &A0, const NekDouble &alpha=0.5)
virtual void	v_GetUFromChars (NekDouble &u, const NekDouble &W1, const NekDouble &W2)
virtual void	v_GetCharIntegral (NekDouble &I, const NekDouble &beta, const NekDouble &A, const NekDouble &A0, const NekDouble &alpha=0.5)
virtual void	v_GetJacobianInverse (NekMatrix< NekDouble > &invJ, const Array< OneD, NekDouble > &Au, const Array< OneD, NekDouble > &uu, const Array< OneD, NekDouble > &beta, const Array< OneD, NekDouble > &A0, const Array< OneD, NekDouble > &alpha, const std::string &type)

Additional Inherited Members
Protected Attributes inherited from Nektar::PulseWavePressureArea
Array< OneD, MultiRegions::ExpListSharedPtr >	m_vessels
LibUtilities::SessionReaderSharedPtr	m_session
NekDouble	m_PExt
NekDouble	m_rho

Detailed Description

Definition at line 50 of file BetaPressureArea.h.

Constructor & Destructor Documentation

BetaPressureArea()

Nektar::BetaPressureArea::BetaPressureArea	(	Array< OneD, MultiRegions::ExpListSharedPtr >	pVessel,
		const LibUtilities::SessionReaderSharedPtr	pSession
	)

Definition at line 48 of file BetaPressureArea.cpp.

    : PulseWavePressureArea(pVessel, pSession)
{
}

~BetaPressureArea()

Nektar::BetaPressureArea::~BetaPressureArea ( )

virtual

Definition at line 55 of file BetaPressureArea.cpp.

{
}

Member Function Documentation

create()

static PulseWavePressureAreaSharedPtr Nektar::BetaPressureArea::create ( Array< OneD, MultiRegions::ExpListSharedPtr > &  pVessel, const LibUtilities::SessionReaderSharedPtr &  pSession  )

inlinestatic

Definition at line 54 of file BetaPressureArea.h.

        {
            return MemoryManager<BetaPressureArea>::AllocateSharedPtr(pVessel,
                                                                      pSession);
        }

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

v_GetAFromChars()

void Nektar::BetaPressureArea::v_GetAFromChars ( NekDouble &  A, const NekDouble &  W1, const NekDouble &  W2, const NekDouble &  beta, const NekDouble &  A0, const NekDouble &  alpha = 0.5  )

protectedvirtual

Implements Nektar::PulseWavePressureArea.

Definition at line 92 of file BetaPressureArea.cpp.

{
    A = pow((W1 - W2) * sqrt(2 * m_rho / beta) / 8 + sqrt(sqrt(A0)), 4);
}

References Nektar::PulseWavePressureArea::m_rho, and tinysimd::sqrt().

v_GetC()

void Nektar::BetaPressureArea::v_GetC ( NekDouble &  c, const NekDouble &  beta, const NekDouble &  A, const NekDouble &  A0, const NekDouble &  alpha = 0.5  )

protectedvirtual

Implements Nektar::PulseWavePressureArea.

Definition at line 66 of file BetaPressureArea.cpp.

{
    c = sqrt(beta / (2 * m_rho)) * sqrt(sqrt(A)); // Elastic
}

References Nektar::PulseWavePressureArea::m_rho, and tinysimd::sqrt().

v_GetCharIntegral()

void Nektar::BetaPressureArea::v_GetCharIntegral ( NekDouble &  I, const NekDouble &  beta, const NekDouble &  A, const NekDouble &  A0, const NekDouble &  alpha = 0.5  )

protectedvirtual

Implements Nektar::PulseWavePressureArea.

Definition at line 105 of file BetaPressureArea.cpp.

{
    NekDouble c  = 0.0;
    NekDouble c0 = 0.0;
 
    GetC(c,  beta, A,  A0);
    GetC(c0, beta, A0, A0);
 
    I = 4 * (c - c0);
}

References Nektar::PulseWavePressureArea::GetC().

v_GetJacobianInverse()

void Nektar::BetaPressureArea::v_GetJacobianInverse ( NekMatrix< NekDouble > &  invJ, const Array< OneD, NekDouble > &  Au, const Array< OneD, NekDouble > &  uu, const Array< OneD, NekDouble > &  beta, const Array< OneD, NekDouble > &  A0, const Array< OneD, NekDouble > &  alpha, const std::string &  type  )

protectedvirtual

Implements Nektar::PulseWavePressureArea.

Definition at line 117 of file BetaPressureArea.cpp.

{
    /*
    In the interest of speed, the inverse of the Jacobians for bifurcations,
    merges and junctions for the beta law have been calculated analytically.
    This can be done for other laws too, or the general formulation can be used
    instead.
    */
 
    NekDouble k = 0.0;
 
    if (type == "Bifurcation")
    {
        NekMatrix<NekDouble> J(6, 6);
        Array<OneD, NekDouble> c(3);
        Array<OneD, NekDouble> K(3);
 
        for (int i = 0; i < 3; ++i)
        {
            GetC(c[i], beta[i], Au[i], A0[i], alpha[i]);
        }
 
        k = c[0] * Au[1] * c[2] + Au[0] * c[2] * c[1] + Au[2] * c[0] * c[1];
        K[0] = (c[0] - uu[0]) * k;
        K[1] = (c[1] + uu[1]) * k;
        K[2] = (c[2] + uu[2]) * k;
 
        invJ.SetValue(0, 0, (-c[1] * uu[0] * c[2] * Au[0] + Au[2] * c[1] * c[0]
                                   * c[0] + Au[1] * c[0] * c[0] * c[2]) / K[0]);
        invJ.SetValue(0, 1, Au[1] * (c[1] - uu[1]) * c[0] * c[2] / K[0]);
        invJ.SetValue(0, 2, Au[2] * (c[2] - uu[2]) * c[0] * c[1] / K[0]);
        invJ.SetValue(0, 3, c[0] * c[1] * c[2] / K[0]);
        invJ.SetValue(0, 4, -0.5 * c[0] * Au[1] * c[2] / K[0]);
        invJ.SetValue(0, 5, -0.5 * Au[2] * c[0] * c[1] / K[0]);
 
        invJ.SetValue(1, 0, Au[0] * (c[0] + uu[0]) * c[1] * c[2] / K[1]);
        invJ.SetValue(1, 1, (c[0] * uu[1] * c[2] * Au[1] + Au[2] * c[0] * c[1]
                                   * c[1] + c[2] * c[1] * c[1] * Au[0]) / K[1]);
        invJ.SetValue(1, 2, -Au[2] * (c[2] - uu[2]) * c[0] * c[1] / K[1]);
        invJ.SetValue(1, 3, -c[0] * c[1] * c[2] / K[1]);
        invJ.SetValue(1, 4, -0.5 * (c[0] * Au[2] + Au[0] * c[2]) * c[1] / K[1]);
        invJ.SetValue(1, 5, 0.5 * Au[2] * c[0] * c[1] / K[1]);
 
        invJ.SetValue(2, 0, Au[0] * (c[0] + uu[0]) * c[1] * c[2] / K[2]);
        invJ.SetValue(2, 1, -Au[1] * (c[1] - uu[1]) * c[0] * c[2] / K[2]);
        invJ.SetValue(2, 2, (c[0] * c[1] * uu[2] * Au[2] + c[0] * Au[1] * c[2]
                                   * c[2] + c[1] * c[2] * c[2] * Au[0]) / K[2]);
        invJ.SetValue(2, 3, -c[0] * c[1] * c[2] / K[2]);
        invJ.SetValue(2, 4, 0.5 * c[0] * Au[1] * c[2] / K[2]);
        invJ.SetValue(2, 5, -0.5 * (Au[1] * c[0] + c[1] * Au[0]) * c[2] / K[2]);
 
        invJ.SetValue(3, 0, Au[0] * (Au[0] * c[2] * c[1] - uu[0] * c[2] * Au[1]
                                                - uu[0] * c[1] * Au[2]) / K[0]);
        invJ.SetValue(3, 1, -Au[0] * Au[1] * (c[1] - uu[1]) * c[2] / K[0]);
        invJ.SetValue(3, 2, -Au[0] * Au[2] * (c[2] - uu[2]) * c[1] / K[0]);
        invJ.SetValue(3, 3, -Au[0] * c[2] * c[1] / K[0]);
        invJ.SetValue(3, 4, 0.5 * Au[0] * Au[1] * c[2] / K[0]);
        invJ.SetValue(3, 5, 0.5 * Au[0] * c[1] * Au[2] / K[0]);
 
        invJ.SetValue(4, 0, Au[0] * Au[1] * (c[0] + uu[0]) * c[2] / K[1]);
        invJ.SetValue(4, 1, -Au[1] * (c[0] * Au[1] * c[2] + c[0] * uu[1] * Au[2]
                                                + c[2] * uu[1] * Au[0]) / K[1]);
        invJ.SetValue(4, 2, -Au[2] * Au[1] * (c[2] - uu[2]) * c[0] / K[1]);
        invJ.SetValue(4, 3, -c[0] * Au[1] * c[2] / K[1]);
        invJ.SetValue(4, 4, -0.5 * Au[1] * (c[0] * Au[2] + Au[0] * c[2]) / K[1]);
        invJ.SetValue(4, 5, 0.5 * Au[2] * Au[1] * c[0] / K[1]);
 
        invJ.SetValue(5, 0, Au[0] * Au[2] * (c[0] + uu[0]) * c[1] / K[2]);
        invJ.SetValue(5, 1, -Au[2] * Au[1] * (c[1] - uu[1]) * c[0] / K[2]);
        invJ.SetValue(5, 2, -Au[2] * (Au[2] * c[0] * c[1] + c[0] * uu[2] * Au[1]
                                                + c[1] * uu[2] * Au[0]) / K[2]);
        invJ.SetValue(5, 3, -Au[2] * c[0] * c[1] / K[2]);
        invJ.SetValue(5, 4, 0.5 * Au[2] * Au[1] * c[0] / K[2]);
        invJ.SetValue(5, 5, -0.5 * Au[2] * (Au[1] * c[0] + c[1] * Au[0]) / K[2]);
    }
    else if (type == "Merge")
    {
        NekMatrix<NekDouble> J(6, 6);
        Array<OneD, NekDouble> c(3);
        Array<OneD, NekDouble> K(3);
 
        for (int i = 0; i < 3; ++i)
        {
            GetC(c[i], beta[i], Au[i], A0[i], alpha[i]);
        }
 
        k = c[0] * Au[1] * c[2] + Au[0] * c[2] * c[1] + Au[2] * c[0] * c[1];
        K[0] = (c[0] - uu[0]) * k;
        K[1] = (c[1] + uu[1]) * k;
        K[2] = (c[2] + uu[2]) * k;
 
        invJ.SetValue(0, 0, (-c[1] * uu[0] * c[2] * Au[0] + Au[2] * c[1] * c[0]
                                   * c[0] + Au[1] * c[0] * c[0] * c[2]) / K[0]);
        invJ.SetValue(0, 1, Au[1] * (c[1] - uu[1]) * c[0] * c[2] / K[0]);
        invJ.SetValue(0, 2, Au[2] * (c[2] - uu[2]) * c[0] * c[1] / K[0]);
        invJ.SetValue(0, 3, c[0] * c[1] * c[2] / K[0]);
        invJ.SetValue(0, 4, -0.5 * c[0] * Au[1] * c[2] / K[0]);
        invJ.SetValue(0, 5, -0.5 * Au[2] * c[0] * c[1] / K[0]);
 
        invJ.SetValue(1, 0, Au[0] * (c[0] + uu[0]) * c[1] * c[2] / K[1]);
        invJ.SetValue(1, 1, (c[0] * uu[1] * c[2] * Au[1] + Au[2] * c[0] * c[1]
                                   * c[1] + c[2] * c[1] * c[1] * Au[0]) / K[1]);
        invJ.SetValue(1, 2, -Au[2] * (c[2] - uu[2]) * c[0] * c[1] / K[1]);
        invJ.SetValue(1, 3, -c[0] * c[1] * c[2] / K[1]);
        invJ.SetValue(1, 4, -0.5 * (c[0] * Au[2] + Au[0] * c[2]) * c[1] / K[1]);
        invJ.SetValue(1, 5, 0.5 * Au[2] * c[0] * c[1] / K[1]);
 
        invJ.SetValue(2, 0, Au[0] * (c[0] - uu[0]) * c[1] * c[2] / K[2]);
        invJ.SetValue(2, 1, -Au[1] * (c[1] + uu[1]) * c[0] * c[2] / K[2]);
        invJ.SetValue(2, 2, -(c[0] * uu[2] * c[1] * Au[2] - Au[1] * c[0] *
                              c[2] * c[2] - c[1] * c[2] * c[2] * Au[0]) / K[2]);
        invJ.SetValue(2, 3, -c[0] * c[1] * c[2] / K[2]);
        invJ.SetValue(2, 4, -0.5 * Au[1] * c[0] * c[2] / K[2]);
        invJ.SetValue(2, 5, 0.5 * (Au[1] * c[0] + Au[0] * c[1]) * c[2] / K[2]);
 
        invJ.SetValue(3, 0, -Au[0] * (Au[0] * c[2] * c[1] + uu[0] * c[2] *
                                         Au[1] + uu[0] * c[1] * Au[2]) / K[0]);
        invJ.SetValue(3, 1, Au[0] * Au[1] * (c[1] + uu[1]) * c[2] / K[0]);
 
        invJ.SetValue(3, 2, -Au[0] * Au[2] * (c[2] - uu[2]) * c[1] / K[0]);
        invJ.SetValue(3, 3, -Au[0] * c[2] * c[1] / K[0]);
        invJ.SetValue(3, 4, 0.5 * Au[0] * Au[1] * c[2] / K[0]);
        invJ.SetValue(3, 5, 0.5 * Au[0] * c[1] * Au[2] / K[0]);
 
        invJ.SetValue(4, 0, Au[0] * Au[1] * (c[0] + uu[0]) * c[2] / K[1]);
        invJ.SetValue(4, 1, -Au[1] * (c[0] * Au[1] * c[2] + c[0] * uu[1] * Au[2]
                                                + c[2] * uu[1] * Au[0]) / K[1]);
        invJ.SetValue(4, 2, -Au[2] * Au[1] * (c[2] - uu[2]) * c[0] / K[1]);
        invJ.SetValue(4, 3, -c[0] * Au[1] * c[2] / K[1]);
        invJ.SetValue(4, 4, -0.5 * Au[1] * (c[0] * Au[2] + Au[0] * c[2]) / K[1]);
        invJ.SetValue(4, 5, 0.5 * Au[2] * Au[1] * c[0] / K[1]);
 
        invJ.SetValue(5, 0, Au[0] * Au[2] * (c[0] + uu[0]) * c[1] / K[2]);
        invJ.SetValue(5, 1, -Au[2] * Au[1] * (c[1] - uu[1]) * c[0] / K[2]);
        invJ.SetValue(5, 2, -Au[2] * (Au[2] * c[0] * c[1] + c[0] * uu[2] * Au[1]
                                                + c[1] * uu[2] * Au[0]) / K[2]);
        invJ.SetValue(5, 3, -Au[2] * c[0] * c[1] / K[2]);
        invJ.SetValue(5, 4, 0.5 * Au[2] * Au[1] * c[0] / K[2]);
        invJ.SetValue(5, 5, -0.5 * Au[2] * (Au[1] * c[0] + c[1] * Au[0]) / K[2]);
    }
    else if (type == "Junction")
    {
        NekMatrix<NekDouble> J(4, 4);
        Array<OneD, NekDouble> c(2);
        Array<OneD, NekDouble> K(2);
 
        for (int i = 0; i < 2; ++i)
        {
            GetC(c[i], beta[i], Au[i], A0[i], alpha[i]);
        }
 
        k = (c[0] * Au[1] + Au[0] * c[1]);
        K[0] = (c[0] + uu[0]) * k;
        K[1] = (c[1] + uu[1]) * k;
 
        invJ.SetValue(0, 0, (Au[1] * c[0] * c[0] - c[1] * uu[0] * Au[0]) / K[0]);
        invJ.SetValue(0, 1, Au[1] * (c[1] - uu[1]) * c[0] / K[0]);
        invJ.SetValue(0, 2, c[0] * c[1] / K[0]);
        invJ.SetValue(0, 3, -0.5 * c[0] * Au[1] / K[0]);
 
        invJ.SetValue(1, 0, Au[0] * (c[0] + uu[0]) * c[1] / K[1]);
        invJ.SetValue(1, 1, (c[0] * uu[1] * Au[1] + c[1] * c[1] * Au[0]) / K[1]);
        invJ.SetValue(1, 2, -c[0] * c[1] / K[1]);
        invJ.SetValue(1, 3, -0.5 * Au[0] * c[1] / K[1]);
 
        invJ.SetValue(2, 0, Au[0] * (Au[0] * c[1] - uu[0] * Au[1]) / K[0]);
        invJ.SetValue(2, 1, -Au[0] * Au[1] * (c[1] - uu[1]) / K[0]);
        invJ.SetValue(2, 2, -Au[0] * c[1] / K[0]);
        invJ.SetValue(2, 3, 0.5 * Au[1] * Au[0] / K[0]);
 
        invJ.SetValue(3, 0, Au[0] * Au[1] * (c[0] + uu[0]) / K[1]);
        invJ.SetValue(3, 1, -Au[1] * (c[0] * Au[1] + uu[1] * Au[0]) / K[1]);
        invJ.SetValue(3, 2, -c[0] * Au[1] / K[1]);
        invJ.SetValue(3, 3, -0.5 * Au[1] * Au[0] / K[1]);
    }
}

References Nektar::PulseWavePressureArea::GetC().

v_GetPressure()

void Nektar::BetaPressureArea::v_GetPressure ( NekDouble &  P, const NekDouble &  beta, const NekDouble &  A, const NekDouble &  A0, const NekDouble &  dAUdx, const NekDouble &  gamma = 0, const NekDouble &  alpha = 0.5  )

protectedvirtual

Implements Nektar::PulseWavePressureArea.

Definition at line 59 of file BetaPressureArea.cpp.

{
    P = m_PExt + beta * (sqrt(A) - sqrt(A0)) - gamma * dAUdx / sqrt(A); // Viscoelasticity
}

References Nektar::PulseWavePressureArea::m_PExt, main::P, and tinysimd::sqrt().

v_GetUFromChars()

void Nektar::BetaPressureArea::v_GetUFromChars ( NekDouble &  u, const NekDouble &  W1, const NekDouble &  W2  )

protectedvirtual

Implements Nektar::PulseWavePressureArea.

Definition at line 99 of file BetaPressureArea.cpp.

{
    u = (W1 + W2) / 2; // Necessarily the case for all tube laws
}

v_GetW1()

void Nektar::BetaPressureArea::v_GetW1 ( NekDouble &  W1, const NekDouble &  u, const NekDouble &  beta, const NekDouble &  A, const NekDouble &  A0, const NekDouble &  alpha = 0.5  )

protectedvirtual

Implements Nektar::PulseWavePressureArea.

Definition at line 72 of file BetaPressureArea.cpp.

{
    NekDouble I = 0.0;
    GetCharIntegral(I, beta, A, A0);
 
    W1 = u + I; // Elastic and assumes u0 = 0
}

References Nektar::PulseWavePressureArea::GetCharIntegral().

v_GetW2()

void Nektar::BetaPressureArea::v_GetW2 ( NekDouble &  W2, const NekDouble &  u, const NekDouble &  beta, const NekDouble &  A, const NekDouble &  A0, const NekDouble &  alpha = 0.5  )

protectedvirtual

Implements Nektar::PulseWavePressureArea.

Definition at line 82 of file BetaPressureArea.cpp.

{
    NekDouble I = 0.0;
    GetCharIntegral(I, beta, A, A0);
 
    W2 = u - I; // Elastic and assumes u0 = 0
}

References Nektar::PulseWavePressureArea::GetCharIntegral().

Member Data Documentation

className

std::string Nektar::BetaPressureArea::className

static

Initial value:

=
    GetPressureAreaFactory().RegisterCreatorFunction(
        "Beta", BetaPressureArea::create,
        "Beta law pressure area relationship for the arterial system")

Definition at line 61 of file BetaPressureArea.h.