Skip to content

GitLab

Commit c9115957 authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Tidy up stimuli code.

parent cb4d39c6
Hide whitespace changes
Inline Side-by-side
solvers/CardiacEPSolver/EquationSystems/Monodomain.cpp
View file @ c9115957
......@@ -73,6 +73,10 @@ namespace Nektar
{
}
/**
*
*/
void Monodomain::v_InitObject()
{
UnsteadySystem::v_InitObject();
......@@ -85,7 +89,8 @@ namespace Nektar
ASSERTL0(vCellModel != "", "Cell Model not specified.");
m_cell = GetCellModelFactory().CreateInstance(vCellModel, m_session, m_fields[0]);
m_cell = GetCellModelFactory().CreateInstance(
vCellModel, m_session, m_fields[0]);
m_intVariables.push_back(0);
......@@ -139,7 +144,9 @@ namespace Nektar
NekDouble f_range = f_max - f_min;
NekDouble o_min = m_session->GetParameter("o_min");
NekDouble o_max = m_session->GetParameter("o_max");
Vmath::Sadd(nq, -f_min, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
Vmath::Sadd(nq, -f_min,
m_vardiff[varCoeffEnum[i]], 1,
m_vardiff[varCoeffEnum[i]], 1);
for (int j = 0; j < nq; ++j)
{
if (m_vardiff[varCoeffEnum[i]][j] < 0)
......@@ -151,10 +158,18 @@ namespace Nektar
m_vardiff[varCoeffEnum[i]][j] = f_range;
}
}
Vmath::Smul(nq, -1.0/f_range, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
Vmath::Sadd(nq, 1.0, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
Vmath::Smul(nq, o_max-o_min, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
Vmath::Sadd(nq, o_min, m_vardiff[varCoeffEnum[i]], 1, m_vardiff[varCoeffEnum[i]], 1);
Vmath::Smul(nq, -1.0/f_range,
m_vardiff[varCoeffEnum[i]], 1,
m_vardiff[varCoeffEnum[i]], 1);
Vmath::Sadd(nq, 1.0,
m_vardiff[varCoeffEnum[i]], 1,
m_vardiff[varCoeffEnum[i]], 1);
Vmath::Smul(nq, o_max-o_min,
m_vardiff[varCoeffEnum[i]], 1,
m_vardiff[varCoeffEnum[i]], 1);
Vmath::Sadd(nq, o_min,
m_vardiff[varCoeffEnum[i]], 1,
m_vardiff[varCoeffEnum[i]], 1);
}
// Transform variable coefficient and write out to file.
......@@ -170,17 +185,6 @@ namespace Nektar
WriteFld(filename.str());
}
}
//
// if (m_session->DefinesParameter("StimulusDuration"))
// {
// ASSERTL0(m_session->DefinesFunction("Stimulus", "u"),
// "Stimulus function not defined.");
// m_session->LoadParameter("StimulusDuration", m_stimDuration);
// }
// else
// {
// m_stimDuration = 0;
// }
// Search through the loaded filters and pass the cell model to any
// CheckpointCellModel filters loaded.
......@@ -263,6 +267,9 @@ namespace Nektar
}
/**
*
*/
void Monodomain::DoOdeRhs(
const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
......@@ -282,10 +289,14 @@ namespace Nektar
}
/**
*
*/
void Monodomain::v_SetInitialConditions(NekDouble initialtime,
bool dumpInitialConditions)
{
EquationSystem::v_SetInitialConditions(initialtime, dumpInitialConditions);
EquationSystem::v_SetInitialConditions(initialtime,
dumpInitialConditions);
m_cell->Initialise();
}
......@@ -297,21 +308,24 @@ namespace Nektar
{
UnsteadySystem::v_PrintSummary(out);
if (m_session->DefinesFunction("d00") &&
m_session->GetFunctionType("d00", "intensity") == LibUtilities::eFunctionTypeExpression)
m_session->GetFunctionType("d00", "intensity")
== LibUtilities::eFunctionTypeExpression)
{
out << "\tDiffusivity-x : "
<< m_session->GetFunction("d00", "intensity")->GetExpression()
<< endl;
}
if (m_session->DefinesFunction("d11") &&
m_session->GetFunctionType("d11", "intensity") == LibUtilities::eFunctionTypeExpression)
m_session->GetFunctionType("d11", "intensity")
== LibUtilities::eFunctionTypeExpression)
{
out << "\tDiffusivity-x : "
<< m_session->GetFunction("d11", "intensity")->GetExpression()
<< endl;
}
if (m_session->DefinesFunction("d22") &&
m_session->GetFunctionType("d22", "intensity") == LibUtilities::eFunctionTypeExpression)
m_session->GetFunctionType("d22", "intensity")
== LibUtilities::eFunctionTypeExpression)
{
out << "\tDiffusivity-x : "
<< m_session->GetFunction("d22", "intensity")->GetExpression()
......@@ -319,6 +333,4 @@ namespace Nektar
}
m_cell->PrintSummary(out);
}
}
solvers/CardiacEPSolver/Examples/crn.xml
View file @ c9115957
......@@ -514,15 +514,6 @@
<P> NumSteps = FinTime/TimeStep </P>
<P> IO_CheckSteps = 0.1/TimeStep </P>
<P> IO_InfoSteps = 1 </P>
<!-- <P> StimulusDuration = 2.0 </P> <!-- ms --
<P> StimulusStrength = 3 </P>
<P> ix = 1.0 </P>
<P> iy = 0.5 </P>
<P> iz = 0.0 </P>
<P> ir = 0.05 </P>
<P> is = 1 </P>
<P> C = StimulusStrength </P>
<P> D = StimulusDuration/4.0 </P>-->
<P> IterativeSolverTolerance = 1e-05 </P>
<P> Chi = 28 </P> <!-- larger: wavefront moves slower -->
<P> Cm = 0.125 </P> <!-- smaller: higher peak mag. of act. -->
......@@ -571,7 +562,6 @@
<p_x2> 0.7</p_x2>
<p_y2> 0.7</p_y2>
<p_z2> 0.0</p_z2>
<!--<p_r1> 0.1</p_r1>-->
<p_is> 1.0</p_is>
<p_strength> 3.0 </p_strength>
......@@ -581,7 +571,7 @@
<S1CYCLELENGTH> 300.0 </S1CYCLELENGTH>
<NUM_S1> 2.0 </NUM_S1>
<S2CYCLELENGTH>100.0 </S2CYCLELENGTH>
</PROTOCOL>
</PROTOCOL>
</STIMULUS>
<STIMULUS ID="0" TYPE = "StimulusCirc">
<p_x1> 9.8</p_x1>
......
solvers/CardiacEPSolver/Examples/strip.xml
View file @ c9115957
......@@ -734,15 +734,6 @@
<P> NumSteps = FinTime/TimeStep </P>
<P> IO_CheckSteps = 0.1/TimeStep </P>
<P> IO_InfoSteps = 1 </P>
<P> StimulusDuration = 2.0 </P> <!-- ms -->
<P> StimulusStrength = 3 </P>
<P> ix = 1.0 </P>
<P> iy = 0.5 </P>
<P> iz = 0.0 </P>
<P> ir = 0.05 </P>
<P> is = 1 </P>
<P> C = StimulusStrength </P>
<P> D = StimulusDuration/4.0 </P>
<P> IterativeSolverTolerance = 1e-05 </P>
<P> Chi = 28 </P> <!-- larger: wavefront moves slower -->
<P> Cm = 0.125 </P> <!-- smaller: higher peak mag. of act. -->
......@@ -755,7 +746,7 @@
<I PROPERTY="Projection" VALUE="Continuous"/>
<I PROPERTY="DiffusionAdvancement" VALUE="Implicit"/>
<I PROPERTY="TimeIntegrationMethod" VALUE="IMEXdirk_3_4_3"/>
<I PROPERTY="GlobalSysSoln" VALUE="DirectStaticCond" /><!--IterativeStaticCond-->
<I PROPERTY="GlobalSysSoln" VALUE="DirectStaticCond" />
</SOLVERINFO>
<VARIABLES>
......@@ -780,7 +771,6 @@
<E VAR="intensity" VALUE="0.13341" />
</FUNCTION>
</CONDITIONS>
<STIMULI>
......@@ -791,7 +781,6 @@
<p_x2> 0.8</p_x2>
<p_y2> 0.8</p_y2>
<p_z2> 0.0</p_z2>
<!--<p_r1> 0.1</p_r1>> -->
<p_is> 1.0</p_is>
<p_strength> 3.0 </p_strength>
......
solvers/CardiacEPSolver/Stimuli/Protocol.cpp
View file @ c9115957
......@@ -33,13 +33,8 @@
//
///////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <CardiacEPSolver/Stimuli/Protocol.h>
#include <StdRegions/StdNodalTriExp.h>
//#include <LibUtilities/LinearAlgebra/Blas.hpp>
namespace Nektar
{
ProtocolFactory& GetProtocolFactory()
......@@ -55,9 +50,9 @@ namespace Nektar
*
* The Stimuli class and derived classes implement a range of stimuli.
* The stimulus contains input stimuli that can be applied throughout the
* domain, on specified regions determined by the derived classes of Stimulus,
* at specified frequencies determined by the derived classes of Protocol.
*
* domain, on specified regions determined by the derived classes of
* Stimulus, at specified frequencies determined by the derived classes of
* Protocol.
*/
/**
......
solvers/CardiacEPSolver/Stimuli/Protocol.h
View file @ c9115957
......@@ -37,10 +37,6 @@
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
//#include <SpatialDomains/SpatialData.h>
#include <MultiRegions/ExpList.h>
#include <StdRegions/StdNodalTriExp.h>
#include <StdRegions/StdNodalTetExp.h>
namespace Nektar
{
......@@ -49,11 +45,13 @@ namespace Nektar
/// A shared pointer to an EquationSystem object
typedef boost::shared_ptr<Protocol> ProtocolSharedPtr;
/// Datatype of the NekFactory used to instantiate classes derived from
/// the EquationSystem class.
typedef LibUtilities::NekFactory< std::string, Protocol,
const LibUtilities::SessionReaderSharedPtr&,
const TiXmlElement*> ProtocolFactory;
ProtocolFactory& GetProtocolFactory();
......@@ -70,8 +68,7 @@ namespace Nektar
void Initialise();
/// Returns amplitude of stimulus (1 or 0) at given time
NekDouble GetAmplitude(
const NekDouble time)
NekDouble GetAmplitude(const NekDouble time)
{
return v_GetAmplitude(time);
}
......@@ -85,10 +82,8 @@ namespace Nektar
protected:
/// Session
LibUtilities::SessionReaderSharedPtr m_session;
/// Transmembrane potential field from PDE system
virtual NekDouble v_GetAmplitude(
const NekDouble time) = 0;
virtual NekDouble v_GetAmplitude(const NekDouble time) = 0;
virtual void v_PrintSummary(std::ostream &out) = 0;
......
solvers/CardiacEPSolver/Stimuli/ProtocolS1.cpp
View file @ c9115957
......@@ -37,25 +37,30 @@
#include <CardiacEPSolver/Stimuli/ProtocolS1.h>
namespace Nektar
{ std::string ProtocolS1::className
= GetProtocolFactory().RegisterCreatorFunction(
{
std::string ProtocolS1::className
= GetProtocolFactory().RegisterCreatorFunction(
"ProtocolS1",
ProtocolS1::create,
"S1 stimulus protocol.");
/**
* @class ProtocolS1
*
* The Stimuli class and derived classes implement a range of stimuli.
* The stimulus contains input stimuli that can be applied throughout the
* domain, on specified regions determined by the derived classes of Stimulus,
* at specified frequencies determined by the derived classes of Protocol.
*
* domain, on specified regions determined by the derived classes of
* Stimulus, at specified frequencies determined by the derived classes of
* Protocol.
*/
/**
* Protocol base class constructor.
*/
ProtocolS1::ProtocolS1(const LibUtilities::SessionReaderSharedPtr& pSession,const TiXmlElement* pXml)
: Protocol(pSession, pXml)
ProtocolS1::ProtocolS1(
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
: Protocol(pSession, pXml)
{
m_session = pSession;
......@@ -64,13 +69,12 @@ namespace Nektar
return;
}
const TiXmlElement *pXmlparameter; //Declaring variable called pxml...
// See if we have parameters defined. They are optional so we go on if not.
//member variables m_* defined in ProtocolS1.h
// Declare temporary XML element pointer
const TiXmlElement *pXmlparameter;
// Read each variable, extract text and convert to floating-point
pXmlparameter = pXml->FirstChildElement("START");
m_start = atof(pXmlparameter->GetText()); //text value within px1, convert to a floating pt and save in m_px1
m_start = atof(pXmlparameter->GetText());
pXmlparameter = pXml->FirstChildElement("DURATION");
m_dur = atof(pXmlparameter->GetText());
......@@ -80,12 +84,9 @@ namespace Nektar
pXmlparameter = pXml->FirstChildElement("NUM_S1");
m_num_s1 = atof(pXmlparameter->GetText());
}
/**
* Initialise the protocol. Allocate workspace and variable storage.
*/
......@@ -94,25 +95,38 @@ namespace Nektar
}
NekDouble ProtocolS1::v_GetAmplitude(
const NekDouble time)
/**
*
*/
NekDouble ProtocolS1::v_GetAmplitude(const NekDouble time)
{
time1 = time - floor((time-m_start)/m_s1cyclelength)*m_s1cyclelength - m_start;
if( time1 > 0 && (m_s1cyclelength * (m_num_s1) + m_start) && time1 < (m_dur))
time1 = time - floor((time-m_start)/m_s1cyclelength)*m_s1cyclelength
- m_start;
if( (time1 > 0) &&
(m_s1cyclelength * (m_num_s1) + m_start) &&
(time1 < m_dur) )
{
return 1.0;
}
return 0.0;
return 0.0;
}
/**
*
*/
void ProtocolS1::v_PrintSummary(std::ostream &out)
{
}
/**
*
*/
void ProtocolS1::v_SetInitialConditions()
{
......
solvers/CardiacEPSolver/Stimuli/ProtocolS1.h
View file @ c9115957
......@@ -35,6 +35,7 @@
#ifndef NEKTAR_SOLVERS_CARDIACEPSOLVER_STIMULI_PROTOCOLS1
#define NEKTAR_SOLVERS_CARDIACEPSOLVER_STIMULI_PROTOCOLS1
#include <CardiacEPSolver/Stimuli/Protocol.h>
namespace Nektar
......@@ -48,8 +49,8 @@ namespace Nektar
public:
/// Creates an instance of this class
static ProtocolSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
{
return MemoryManager<ProtocolS1>
::AllocateSharedPtr(pSession, pXml);
......@@ -59,7 +60,7 @@ namespace Nektar
static std::string className;
ProtocolS1(const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml);
const TiXmlElement* pXml);
virtual ~ProtocolS1() {}
......@@ -73,8 +74,7 @@ namespace Nektar
NekDouble m_s1cyclelength;
NekDouble time1;
virtual NekDouble v_GetAmplitude(
const NekDouble time);
virtual NekDouble v_GetAmplitude(const NekDouble time);
virtual void v_PrintSummary(std::ostream &out);
......
solvers/CardiacEPSolver/Stimuli/ProtocolS1S2.h
View file @ c9115957
......@@ -48,18 +48,18 @@ namespace Nektar
public:
/// Creates an instance of this class
static ProtocolSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
{
return MemoryManager<ProtocolS1S2>
::AllocateSharedPtr(pSession, pXml);
::AllocateSharedPtr(pSession, pXml);
}
/// Name of class
static std::string className;
ProtocolS1S2(const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml);
const TiXmlElement* pXml);
virtual ~ProtocolS1S2() {}
......@@ -75,8 +75,7 @@ namespace Nektar
NekDouble m_s2start;
NekDouble time1;
virtual NekDouble v_GetAmplitude(
const NekDouble time);
virtual NekDouble v_GetAmplitude(const NekDouble time);
virtual void v_PrintSummary(std::ostream &out);
......@@ -85,4 +84,4 @@ namespace Nektar
}
#endif /* ProtocolS1S2_H_ */
#endif
solvers/CardiacEPSolver/Stimuli/ProtocolSingle.cpp
View file @ c9115957
......@@ -37,25 +37,31 @@
#include <CardiacEPSolver/Stimuli/ProtocolSingle.h>
namespace Nektar
{ std::string ProtocolSingle::className
= GetProtocolFactory().RegisterCreatorFunction(
{
std::string ProtocolSingle::className
= GetProtocolFactory().RegisterCreatorFunction(
"ProtocolSingle",
ProtocolSingle::create,
"Single stimulus protocol.");
/**
* @class Protocol
*
* The Stimuli class and derived classes implement a range of stimuli.
* The stimulus contains input stimuli that can be applied throughout the
* domain, on specified regions determined by the derived classes of Stimulus,
* at specified frequencies determined by the derived classes of Protocol.
* domain, on specified regions determined by the derived classes of
* Stimulus, at specified frequencies determined by the derived classes of
* Protocol.
*
*/
/**
* Protocol base class constructor.
*/
ProtocolSingle::ProtocolSingle(const LibUtilities::SessionReaderSharedPtr& pSession,const TiXmlElement* pXml)
: Protocol(pSession, pXml)
ProtocolSingle::ProtocolSingle(
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
: Protocol(pSession, pXml)
{
m_session = pSession;
......@@ -65,22 +71,16 @@ namespace Nektar
}
const TiXmlElement *pXmlparameter; //Declaring variable called pxml...
// See if we have parameters defined. They are optional so we go on if not.
//member variables m_* defined in ProtocolSingle.h
const TiXmlElement *pXmlparameter;
pXmlparameter = pXml->FirstChildElement("START");
m_start = atof(pXmlparameter->GetText()); //text value within px1, convert to a floating pt and save in m_px1
m_start = atof(pXmlparameter->GetText());
pXmlparameter = pXml->FirstChildElement("DURATION");
m_dur = atof(pXmlparameter->GetText());
}
/**
* Initialise the protocol. Allocate workspace and variable storage.
*/
......@@ -90,6 +90,9 @@ namespace Nektar
}
/**
*
*/
NekDouble ProtocolSingle::v_GetAmplitude(
const NekDouble time)
{
......@@ -104,11 +107,19 @@ namespace Nektar
}
/**
*
*/
void ProtocolSingle::v_PrintSummary(std::ostream &out)
{
}
/**
*
*/
void ProtocolSingle::v_SetInitialConditions()
{
......
solvers/CardiacEPSolver/Stimuli/ProtocolSingle.h
View file @ c9115957
......@@ -35,6 +35,7 @@
#ifndef NEKTAR_SOLVERS_CARDIACEPSOLVER_STIMULI_PROTOCOLSINGLE
#define NEKTAR_SOLVERS_CARDIACEPSOLVER_STIMULI_PROTOCOLSINGLE
#include <CardiacEPSolver/Stimuli/Protocol.h>
namespace Nektar
......@@ -48,18 +49,18 @@ namespace Nektar
public:
/// Creates an instance of this class
static ProtocolSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
{
return MemoryManager<ProtocolSingle>
::AllocateSharedPtr(pSession, pXml);
::AllocateSharedPtr(pSession, pXml);
}
/// Name of class
static std::string className;
ProtocolSingle(const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml);
const TiXmlElement* pXml);
virtual ~ProtocolSingle() {}
......@@ -70,8 +71,7 @@ namespace Nektar
NekDouble m_start;
NekDouble m_dur;
virtual NekDouble v_GetAmplitude(
const NekDouble time);
virtual NekDouble v_GetAmplitude(const NekDouble time);
virtual void v_PrintSummary(std::ostream &out);
......@@ -80,4 +80,4 @@ namespace Nektar
}
#endif /* PROTOCOLSINGLE_H_ */
#endif
solvers/CardiacEPSolver/Stimuli/Stimulus.cpp
View file @ c9115957
......@@ -38,9 +38,6 @@
#include <CardiacEPSolver/Stimuli/Stimulus.h>
#include <StdRegions/StdNodalTriExp.h>
//#include <LibUtilities/LinearAlgebra/Blas.hpp>
namespace Nektar
{