Skip to content

GitLab

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

Merge branch 'feature/cardiac-stim' of localhost:nektar

parents a69acfcc c9115957
Expand all Hide whitespace changes
Inline Side-by-side
solvers/CardiacEPSolver/CMakeLists.txt
View file @ e2f24cc5
......@@ -13,6 +13,13 @@ SET(CardiacEPSolverSource
./CellModels/TenTusscher06M.cpp
./CellModels/TenTusscher06Endo.cpp
./Filters/FilterCheckpointCellModel.cpp
./Stimuli/Stimulus.cpp
./Stimuli/StimulusCircle.cpp
./Stimuli/StimulusRect.cpp
./Stimuli/Protocol.cpp
./Stimuli/ProtocolSingle.cpp
./Stimuli/ProtocolS1.cpp
./Stimuli/ProtocolS1S2.cpp
)
ADD_SOLVER_EXECUTABLE(CardiacEPSolver solvers-extra
......
solvers/CardiacEPSolver/EquationSystems/Monodomain.cpp
View file @ e2f24cc5
......@@ -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.
......@@ -171,17 +186,6 @@ namespace Nektar
}
}
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.
int k = 0;
......@@ -198,6 +202,9 @@ namespace Nektar
}
}
// Load stimuli
m_stimulus = Stimulus::LoadStimuli(m_session, m_fields[0]);
if (!m_explicitDiffusion)
{
m_ode.DefineImplicitSolve (&Monodomain::DoImplicitSolve, this);
......@@ -260,6 +267,9 @@ namespace Nektar
}
/**
*
*/
void Monodomain::DoOdeRhs(
const Array<OneD, const Array<OneD, NekDouble> >&inarray,
Array<OneD, Array<OneD, NekDouble> >&outarray,
......@@ -270,30 +280,23 @@ namespace Nektar
m_cell->TimeIntegrate(inarray, outarray, time);
if (m_stimDuration > 0 && time < m_stimDuration)
{
Array<OneD,NekDouble> x0(nq);
Array<OneD,NekDouble> x1(nq);
Array<OneD,NekDouble> x2(nq);
Array<OneD,NekDouble> result(nq);
// get the coordinates
m_fields[0]->GetCoords(x0,x1,x2);
LibUtilities::EquationSharedPtr ifunc
= m_session->GetFunction("Stimulus", "u");
ifunc->Evaluate(x0,x1,x2,time, result);
Vmath::Vadd(nq, outarray[0], 1, result, 1, outarray[0], 1);
for (unsigned int i = 0; i < m_stimulus.size(); ++i)
{
m_stimulus[i]->Update(outarray, time);
}
Vmath::Smul(nq, 1.0/m_capMembrane, outarray[0], 1, outarray[0], 1);
}
/**
*
*/
void Monodomain::v_SetInitialConditions(NekDouble initialtime,
bool dumpInitialConditions)
{
EquationSystem::v_SetInitialConditions(initialtime, dumpInitialConditions);
EquationSystem::v_SetInitialConditions(initialtime,
dumpInitialConditions);
m_cell->Initialise();
}
......@@ -305,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()
......@@ -327,5 +333,4 @@ namespace Nektar
}
m_cell->PrintSummary(out);
}
}
solvers/CardiacEPSolver/EquationSystems/Monodomain.h
View file @ e2f24cc5
......@@ -38,6 +38,7 @@
#include <SolverUtils/UnsteadySystem.h>
#include <CardiacEPSolver/CellModels/CellModel.h>
#include <CardiacEPSolver/Stimuli/Stimulus.h>
using namespace Nektar::SolverUtils;
......@@ -97,6 +98,8 @@ namespace Nektar
/// Cell model.
CellModelSharedPtr m_cell;
std::vector<StimulusSharedPtr> m_stimulus;
/// Variable diffusivity
StdRegions::VarCoeffMap m_vardiff;
......@@ -106,6 +109,7 @@ namespace Nektar
/// Stimulus current
NekDouble m_stimDuration;
void LoadStimuli();
};
}
......
solvers/CardiacEPSolver/Examples/crn.xml
View file @ e2f24cc5
<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:noNamespaceSchemaLocation="http://www.nektar.info/schema/nektar.xsd">
<GEOMETRY DIM="2" SPACE="2">
<VERTEX>
<V ID="0">2.285e+00 3.636e-01 0.000e+00</V>
<V ID="1">2.083e+00 0.000e+00 0.000e+00</V>
......@@ -505,23 +505,15 @@
<EXPANSIONS>
<E COMPOSITE="C[0]" NUMMODES="7" FIELDS="u,v" TYPE="MODIFIED" />
</EXPANSIONS>
<CONDITIONS>
<PARAMETERS>
<P> TimeStep = 0.02 </P>
<P> FinTime = 500 </P>
<P> FinTime = 200 </P>
<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. -->
......@@ -550,24 +542,52 @@
<N VAR="u" VALUE="0.0" />
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="InitialConditions">
<E VAR="u" VALUE="-81.19" />
</FUNCTION>
<FUNCTION NAME="IsotropicConductivity">
<E VAR="intensity" VALUE="0.13341" />
</FUNCTION>
<EXPRESSIONS>
<E NAME="P" VALUE="1.0/(1.0 + exp( 32*(x-1)))" />
</EXPRESSIONS>
<FUNCTION NAME="Stimulus">
<E VAR="u" VALUE="C*P" />
</FUNCTION>
</CONDITIONS>
<STIMULI>
<STIMULUS ID="0" TYPE = "StimulusRect">
<p_x1> 0.3</p_x1>
<p_y1> 0.3</p_y1>
<p_z1> 0.0</p_z1>
<p_x2> 0.7</p_x2>
<p_y2> 0.7</p_y2>
<p_z2> 0.0</p_z2>
<p_is> 1.0</p_is>
<p_strength> 3.0 </p_strength>
<PROTOCOL TYPE = "ProtocolS1S2">
<START> 2.0 </START>
<DURATION> 2.0 </DURATION>
<S1CYCLELENGTH> 300.0 </S1CYCLELENGTH>
<NUM_S1> 2.0 </NUM_S1>
<S2CYCLELENGTH>100.0 </S2CYCLELENGTH>
</PROTOCOL>
</STIMULUS>
<STIMULUS ID="0" TYPE = "StimulusCirc">
<p_x1> 9.8</p_x1>
<p_y1> 0.5</p_y1>
<p_z1> 0.0</p_z1>
<p_r1> 0.1</p_r1>
<p_is> 1.0</p_is>
<p_strength> 5.0 </p_strength>
<PROTOCOL TYPE = "ProtocolSingle">
<START> 280.0 </START>
<DURATION> 2.0 </DURATION>
</PROTOCOL>
</STIMULUS>
</STIMULI>
<FILTERS>
<FILTER TYPE="HistoryPoints">
<PARAM NAME="OutputFile">crn.his</PARAM>
......
solvers/CardiacEPSolver/Examples/strip.xml 0 → 100644
View file @ e2f24cc5
This diff is collapsed.
solvers/CardiacEPSolver/Stimuli/Protocol.cpp 0 → 100644
View file @ e2f24cc5
///////////////////////////////////////////////////////////////////////////////
//
// File Protocol.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Protocol base class.
//
///////////////////////////////////////////////////////////////////////////////
#include <CardiacEPSolver/Stimuli/Protocol.h>
namespace Nektar
{
ProtocolFactory& GetProtocolFactory()
{
typedef Loki::SingletonHolder<ProtocolFactory,
Loki::CreateUsingNew,
Loki::NoDestroy > Type;
return Type::Instance();
}
/**
* @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.
*/
/**
* Protocol base class constructor.
*/
Protocol::Protocol(const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
{
m_session = pSession;
}
/**
* Initialise the protocol. Allocate workspace and variable storage.
*/
void Protocol::Initialise()
{
}
}
solvers/CardiacEPSolver/Stimuli/Protocol.h 0 → 100644
View file @ e2f24cc5
///////////////////////////////////////////////////////////////////////////////
//
// File Protocol.h
//
// For more information, please see: http://www.nektar.info
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Stimulus protocol base class.
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERS_CARDIACEPSOLVER_STIMULI_PROTOCOL
#define NEKTAR_SOLVERS_CARDIACEPSOLVER_STIMULI_PROTOCOL
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
namespace Nektar
{
// Forward declaration
class Protocol;
/// 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();
/// Protocol base class.
class Protocol
{
public:
Protocol(const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml);
virtual ~Protocol() {}
/// Initialise the protocol storage and set initial conditions
void Initialise();
/// Returns amplitude of stimulus (1 or 0) at given time
NekDouble GetAmplitude(const NekDouble time)
{
return v_GetAmplitude(time);
}
/// Print a summary of the cell model
void PrintSummary(std::ostream &out)
{
v_PrintSummary(out);
}
protected:
/// Session
LibUtilities::SessionReaderSharedPtr m_session;
virtual NekDouble v_GetAmplitude(const NekDouble time) = 0;
virtual void v_PrintSummary(std::ostream &out) = 0;
};
}
#endif /*PROTOCOL_H_ */
solvers/CardiacEPSolver/Stimuli/ProtocolS1.cpp 0 → 100644
View file @ e2f24cc5
///////////////////////////////////////////////////////////////////////////////
//
// File ProtocolS1.cpp
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: S1 impulse protocol.
//
///////////////////////////////////////////////////////////////////////////////
#include <tinyxml/tinyxml.h>
#include <CardiacEPSolver/Stimuli/ProtocolS1.h>
namespace Nektar
{
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.
*/
/**
* Protocol base class constructor.
*/
ProtocolS1::ProtocolS1(
const LibUtilities::SessionReaderSharedPtr& pSession,
const TiXmlElement* pXml)
: Protocol(pSession, pXml)
{
m_session = pSession;
if (!pXml)
{
return;
}
// 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());
pXmlparameter = pXml->FirstChildElement("DURATION");
m_dur = atof(pXmlparameter->GetText());
pXmlparameter = pXml->FirstChildElement("S1CYCLELENGTH");
m_s1cyclelength = atof(pXmlparameter->GetText());
pXmlparameter = pXml->FirstChildElement("NUM_S1");
m_num_s1 = atof(pXmlparameter->GetText());
}
/**
* Initialise the protocol. Allocate workspace and variable storage.
*/
void ProtocolS1::Initialise()
{
}
/**
*
*/
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) )
{
return 1.0;
}
return 0.0;
}
/**
*
*/
void ProtocolS1::v_PrintSummary(std::ostream &out)
{
}
/**
*
*/
void ProtocolS1::v_SetInitialConditions()