#include <FilterElectrogram.h>

Inheritance diagram for Nektar::FilterElectrogram:

Collaboration diagram for Nektar::FilterElectrogram:

Public Member Functions
	FilterElectrogram (const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
	Electrogram filter constructor.
virtual	~FilterElectrogram ()
	Electrogram filter destructor.
Public Member Functions inherited from Nektar::SolverUtils::Filter
SOLVER_UTILS_EXPORT	Filter (const LibUtilities::SessionReaderSharedPtr &pSession)
virtual SOLVER_UTILS_EXPORT	~Filter ()
SOLVER_UTILS_EXPORT void	Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
SOLVER_UTILS_EXPORT void	Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
SOLVER_UTILS_EXPORT void	Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
SOLVER_UTILS_EXPORT bool	IsTimeDependent ()

Static Public Member Functions
static SolverUtils::FilterSharedPtr	create (const LibUtilities::SessionReaderSharedPtr &pSession, const std::map< std::string, std::string > &pParams)
	Creates an instance of this class.

Static Public Attributes
static std::string	className
	Name of the class.

Protected Member Functions
virtual void	v_Initialise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
	Initialises the electrogram filter and open output file.
virtual void	v_Update (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
	Compute extracellular potential at egm points at current time.
virtual void	v_Finalise (const Array< OneD, const MultiRegions::ExpListSharedPtr > &pFields, const NekDouble &time)
	Finalise the electrogram filter and close output file.
virtual bool	v_IsTimeDependent ()
	Filter is time-dependent and should be called at each time-step.

Private Attributes
Array< OneD, Array< OneD, NekDouble > >	m_grad_R_x
	Gradient of the radius from each electrogram point in x-direction.
Array< OneD, Array< OneD, NekDouble > >	m_grad_R_y
	Gradient of the radius from each electrogram point in y-direction.
Array< OneD, Array< OneD, NekDouble > >	m_grad_R_z
	Gradient of the radius from each electrogram point in z-direction.
SpatialDomains::PointGeomVector	m_electrogramPoints
	List of electrogram points.
unsigned int	m_index
	Counts number of calls to update (number of timesteps)
unsigned int	m_outputFrequency
	Number of timesteps between outputs.
std::string	m_outputFile
	Filename to output electrogram data to.
std::ofstream	m_outputStream
	Output file stream for electrogram data.
std::stringstream	m_electrogramStream
	Point coordinate input string.

Friends
class	MemoryManager< FilterElectrogram >

Additional Inherited Members
Protected Attributes inherited from Nektar::SolverUtils::Filter
LibUtilities::SessionReaderSharedPtr	m_session

Detailed Description

Definition at line 44 of file FilterElectrogram.h.

Constructor & Destructor Documentation

Nektar::FilterElectrogram::FilterElectrogram	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::map< std::string, std::string > &	pParams
	)

Electrogram filter constructor.

Definition at line 50 of file FilterElectrogram.cpp.

References ASSERTL0, m_electrogramStream, m_index, m_outputFile, m_outputFrequency, and Nektar::SolverUtils::Filter::m_session.

        : Filter(pSession)
    {
        if (pParams.find("OutputFile") == pParams.end())
        {
            m_outputFile = m_session->GetSessionName();
        }
        else
        {
            ASSERTL0(!(pParams.find("OutputFile")->second.empty()),
                     "Missing parameter 'OutputFile'.");
            m_outputFile = pParams.find("OutputFile")->second;
        }
        if (!(m_outputFile.length() >= 4
              && m_outputFile.substr(m_outputFile.length() - 4) == ".ecg"))
        {
            m_outputFile += ".ecg";
        }
        if (pParams.find("OutputFrequency") == pParams.end())
        {
            m_outputFrequency = 1;
        }
        else
        {
            m_outputFrequency =
                    atoi(pParams.find("OutputFrequency")->second.c_str());
        }
        ASSERTL0(pParams.find("Points") != pParams.end(),
                 "Missing parameter 'Points'.");
        m_electrogramStream.str(pParams.find("Points")->second);
        m_index = 0;
    }

Nektar::FilterElectrogram::~FilterElectrogram ( )

virtual

Electrogram filter destructor.

Definition at line 91 of file FilterElectrogram.cpp.

    {
    }

Member Function Documentation

static SolverUtils::FilterSharedPtr Nektar::FilterElectrogram::create	(	const LibUtilities::SessionReaderSharedPtr &	pSession,
		const std::map< std::string, std::string > &	pParams
	)

inlinestatic

Creates an instance of this class.

Definition at line 50 of file FilterElectrogram.h.

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

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

protectedvirtual

Finalise the electrogram filter and close output file.

Implements Nektar::SolverUtils::Filter.

Definition at line 264 of file FilterElectrogram.cpp.

References m_outputStream.

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

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

protectedvirtual

Initialises the electrogram filter and open output file.

Implements Nektar::SolverUtils::Filter.

Definition at line 100 of file FilterElectrogram.cpp.

References ASSERTL0, m_electrogramPoints, m_electrogramStream, m_grad_R_x, m_grad_R_y, m_grad_R_z, m_index, m_outputFile, m_outputStream, Nektar::SolverUtils::Filter::m_session, npts, Vmath::Sadd(), Vmath::Sdiv(), v_Update(), Vmath::Vsqrt(), Vmath::Vvtvp(), and Vmath::Vvtvvtp().

    {
        ASSERTL0(!m_electrogramStream.fail(),
                 "No history points in stream.");
        m_index = 0;
        Array<OneD, NekDouble>  gloCoord(3,0.0);
        LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
        // Read electrogram points
        // Always use dim = 3 to allow electrode to be above surface
        const int dim = 3;
        int i = 0;
        while (!m_electrogramStream.fail())
        {
            m_electrogramStream >> gloCoord[0] >> gloCoord[1] >> gloCoord[2];
            if (!m_electrogramStream.fail())
            {
                SpatialDomains::PointGeomSharedPtr vert
                    = MemoryManager<SpatialDomains::PointGeom>
                    ::AllocateSharedPtr(dim, i, gloCoord[0],
                                        gloCoord[1], gloCoord[2]);
                m_electrogramPoints.push_back(vert);
                ++i;
            }
        }
        if (vComm->GetRank() == 0)
        {
            // Open output stream
            m_outputStream.open(m_outputFile.c_str());
            m_outputStream << "# Electrogram data for variables (:";
            for (i = 0; i < pFields.num_elements(); ++i)
            {
                m_outputStream << m_session->GetVariable(i) <<",";
            }
            m_outputStream << ") at points:" << endl;
            for (i = 0; i < m_electrogramPoints.size(); ++i)
            {
                m_electrogramPoints[i]->GetCoords(  gloCoord[0],
                                                gloCoord[1],
                                                gloCoord[2]);
                m_outputStream << "# \t" << i;
                m_outputStream.width(8);
                m_outputStream << gloCoord[0];
                m_outputStream.width(8);
                m_outputStream << gloCoord[1];
                m_outputStream.width(8);
                m_outputStream << gloCoord[2];
                m_outputStream << endl;
            }
        }
        // Compute the distance function for each electrogram point
        const unsigned int nq   = pFields[0]->GetNpoints();
        const unsigned int npts = m_electrogramPoints.size();
        NekDouble px, py, pz;
        m_grad_R_x = Array<OneD, Array<OneD, NekDouble> >(npts);
        m_grad_R_y = Array<OneD, Array<OneD, NekDouble> >(npts);
        m_grad_R_z = Array<OneD, Array<OneD, NekDouble> >(npts);
        Array<OneD, NekDouble> x(nq);
        Array<OneD, NekDouble> y(nq);
        Array<OneD, NekDouble> z(nq);
        Array<OneD, NekDouble> oneOverR(nq);
        for (unsigned int i = 0; i < npts; ++i)
        {
            m_grad_R_x[i] = Array<OneD, NekDouble>(nq);
            m_grad_R_y[i] = Array<OneD, NekDouble>(nq);
            m_grad_R_z[i] = Array<OneD, NekDouble>(nq);
            // Compute 1/R
            m_electrogramPoints[i]->GetCoords(px,py,pz);
            pFields[0]->GetCoords(x,y,z);
            Vmath::Sadd   (nq, -px, x, 1, x, 1);
            Vmath::Sadd   (nq, -py, y, 1, y, 1);
            Vmath::Sadd   (nq, -pz, z, 1, z, 1);
            Vmath::Vvtvvtp(nq, x, 1, x, 1, y, 1, y, 1, oneOverR, 1);
            Vmath::Vvtvp  (nq, z, 1, z, 1, oneOverR, 1, oneOverR, 1);
            Vmath::Vsqrt  (nq, oneOverR, 1, oneOverR, 1);
            Vmath::Sdiv   (nq, 1.0, oneOverR, 1, oneOverR, 1);
            // Compute grad 1/R
            pFields[0]->PhysDeriv(oneOverR, m_grad_R_x[i], m_grad_R_y[i],
                                            m_grad_R_z[i]);
        }
        // Compute electrogram point for initial condition
        v_Update(pFields, time);
    }

bool Nektar::FilterElectrogram::v_IsTimeDependent ( )

protectedvirtual

Filter is time-dependent and should be called at each time-step.

Implements Nektar::SolverUtils::Filter.

Definition at line 278 of file FilterElectrogram.cpp.

    {
        return true;
    }

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

protectedvirtual

Compute extracellular potential at egm points at current time.

Implements Nektar::SolverUtils::Filter.

Definition at line 206 of file FilterElectrogram.cpp.

References m_electrogramPoints, m_grad_R_x, m_grad_R_y, m_grad_R_z, m_index, m_outputFrequency, m_outputStream, Nektar::LibUtilities::ReduceSum, Vmath::Vvtvp(), and Vmath::Vvtvvtp().

Referenced by v_Initialise().

    {
        // Only output every m_outputFrequency.
        if ((m_index++) % m_outputFrequency)
        {
            return;
        }
        const unsigned int nq = pFields[0]->GetNpoints();
        const unsigned int npoints = m_electrogramPoints.size();
        LibUtilities::CommSharedPtr vComm = pFields[0]->GetComm();
        unsigned int i = 0;
        Array<OneD, NekDouble> e(npoints);
        // Compute grad V
        Array<OneD, NekDouble> grad_V_x(nq), grad_V_y(nq), grad_V_z(nq);
        pFields[0]->PhysDeriv(pFields[0]->GetPhys(),
                              grad_V_x, grad_V_y, grad_V_z);
        for (i = 0; i < npoints; ++i)
        {
            // Multiply together
            Array<OneD, NekDouble> output(nq);
            Vmath::Vvtvvtp(nq, m_grad_R_x[i], 1, grad_V_x, 1, m_grad_R_y[i], 1,
                               grad_V_y,      1, output,   1);
            Vmath::Vvtvp  (nq, m_grad_R_z[i], 1, grad_V_z, 1, output, 1,
                               output,        1);
            e[i] = pFields[0]->Integral(output);
        }
        // Exchange history data
        // This could be improved to reduce communication but works for now
        vComm->AllReduce(e, LibUtilities::ReduceSum);
        // Only the root process writes out electrogram data
        if (vComm->GetRank() == 0)
        {
            m_outputStream.width(8);
            m_outputStream << setprecision(6) << time;
            // Write data values point by point
            for (i = 0; i < m_electrogramPoints.size(); ++i)
            {
                m_outputStream.width(25);
                m_outputStream << setprecision(16) << e[i];
            }
            m_outputStream << endl;
        }
    }

Friends And Related Function Documentation

friend class MemoryManager< FilterElectrogram >

friend

Definition at line 47 of file FilterElectrogram.h.

Member Data Documentation

std::string Nektar::FilterElectrogram::className

static

Initial value:

            SolverUtils::GetFilterFactory().RegisterCreatorFunction(
                    "Electrogram",
                    FilterElectrogram::create)

Name of the class.

Definition at line 60 of file FilterElectrogram.h.

SpatialDomains::PointGeomVector Nektar::FilterElectrogram::m_electrogramPoints

private

List of electrogram points.

Definition at line 94 of file FilterElectrogram.h.

Referenced by v_Initialise(), and v_Update().

std::stringstream Nektar::FilterElectrogram::m_electrogramStream

private

Point coordinate input string.

Definition at line 104 of file FilterElectrogram.h.

Referenced by FilterElectrogram(), and v_Initialise().