Commit ce8a6cba authored by Pavel Burovskiy's avatar Pavel Burovskiy
Browse files

Merge branch 'master' into fix/cg-svtvp

parents 39b19188 64ffa599
......@@ -1434,6 +1434,11 @@ namespace Nektar
{
m_geometricInfo.clear();
if (!geometry)
{
return;
}
TiXmlElement *geometricInfoElement =
geometry->FirstChildElement("GEOMINFO");
......
......@@ -21,6 +21,7 @@ IF( NEKTAR_SOLVER_CARDIAC_EP )
./Stimuli/Stimulus.cpp
./Stimuli/StimulusCircle.cpp
./Stimuli/StimulusRect.cpp
./Stimuli/StimulusPoint.cpp
./Stimuli/Protocol.cpp
./Stimuli/ProtocolSingle.cpp
./Stimuli/ProtocolS1.cpp
......
......@@ -125,7 +125,7 @@ namespace Nektar
void CellModel::Initialise()
{
ASSERTL1(m_nvar > 0, "Cell model must have at least 1 variable.");
m_cellSol = Array<OneD, Array<OneD, NekDouble> >(m_nvar);
m_wsp = Array<OneD, Array<OneD, NekDouble> >(m_nvar);
for (unsigned int i = 0; i < m_nvar; ++i)
......
///////////////////////////////////////////////////////////////////////////////
//
// File StimulusPoint.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: Rectangular stimulus class
//
///////////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <tinyxml/tinyxml.h>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <CardiacEPSolver/Stimuli/StimulusPoint.h>
namespace Nektar
{
std::string StimulusPoint::className
= GetStimulusFactory().RegisterCreatorFunction(
"StimulusPoint",
StimulusPoint::create,
"Point stimulus.");
/**
* @class StimulusPoint
*
* The Stimulus 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.
*/
/**
* Stimulus base class constructor.
*/
StimulusPoint::StimulusPoint(
const LibUtilities::SessionReaderSharedPtr& pSession,
const MultiRegions::ExpListSharedPtr& pField,
const TiXmlElement* pXml)
: Stimulus(pSession, pField, pXml)
{
m_session = pSession;
m_field = pField;
m_nq = pField->GetTotPoints();
if (!pXml)
{
return;
}
const TiXmlElement *pXmlparameter;
pXmlparameter = pXml->FirstChildElement("p_strength");
m_strength = atof(pXmlparameter->GetText());
}
/**
* Initialise the stimulus. Allocate workspace and variable storage.
*/
void StimulusPoint::Initialise()
{
}
/**
*
*/
void StimulusPoint::v_Update(
Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time)
{
// Get the protocol amplitude
NekDouble v_amp = m_Protocol->GetAmplitude(time) * m_strength;
outarray[0][0] += v_amp;
}
/**
*
*/
void StimulusPoint::v_PrintSummary(std::ostream &out)
{
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File StimulusPoint.h
//
// 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: Rectangular stimulus header file
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERS_CARDIACEPSOLVER_STIMULI_STIMULUSRECT
#define NEKTAR_SOLVERS_CARDIACEPSOLVER_STIMULI_STIMULUSRECT
#include <LibUtilities/BasicUtils/NekFactory.hpp>
#include <LibUtilities/BasicUtils/SessionReader.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <MultiRegions/ExpList.h>
#include <CardiacEPSolver/Stimuli/Stimulus.h>
namespace Nektar
{
/// Protocol base class.
class StimulusPoint: public Stimulus
{
public:
/// Creates an instance of this class
static StimulusSharedPtr create(
const LibUtilities::SessionReaderSharedPtr& pSession,
const MultiRegions::ExpListSharedPtr& pField,
const TiXmlElement* pXml)
{
return MemoryManager<StimulusPoint>
::AllocateSharedPtr(pSession, pField, pXml);
}
/// Name of class
static std::string className;
StimulusPoint(const LibUtilities::SessionReaderSharedPtr& pSession,
const MultiRegions::ExpListSharedPtr& pField,
const TiXmlElement* pXml);
virtual ~StimulusPoint() {}
/// Initialise the stimulus storage and set initial conditions
void Initialise();
protected:
NekDouble m_strength;
virtual void v_Update(Array<OneD, Array<OneD, NekDouble> >&outarray,
const NekDouble time);
virtual void v_PrintSummary(std::ostream &out);
};
}
#endif
ADD_SUBDIRECTORY(PrePacing)
ADD_SUBDIRECTORY(FibreRegister)
ADD_SUBDIRECTORY(NavXToVtk)
ADD_SUBDIRECTORY(FibreToNektar)
SET(LinkLibraries
optimized MultiRegions debug MultiRegions-g
optimized LocalRegions debug LocalRegions-g
optimized SpatialDomains debug SpatialDomains-g
optimized StdRegions debug StdRegions-g
optimized LibUtilities debug LibUtilities-g
optimized ${Boost_THREAD_LIBRARY_RELEASE} debug ${Boost_THREAD_LIBRARY_DEBUG}
optimized ${Boost_IOSTREAMS_LIBRARY_RELEASE} debug ${Boost_IOSTREAMS_LIBRARY_DEBUG}
optimized ${Boost_ZLIB_LIBRARY_RELEASE} debug ${Boost_ZLIB_LIBRARY_DEBUG}
optimized ${TINYXML_LIB} debug ${TINYXML_LIB}
optimized sbtk debug sbtk
)
SET(PP_SOURCES ./Prepacing.cpp
../../CellModels/CellModel.cpp
../../CellModels/CourtemancheRamirezNattel98.cpp
../../Stimuli/Stimulus.cpp
../../Stimuli/StimulusPoint.cpp
../../Stimuli/Protocol.cpp
../../Stimuli/ProtocolS1S2.cpp)
ADD_SOLVER_EXECUTABLE(PrePacing solvers-extra ${PP_SOURCES})
TARGET_LINK_LIBRARIES(PrePacing ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(PrePacing)
#include <SpatialDomains/MeshComponents.h>
#include <MultiRegions/ExpList0D.h>
#include <CardiacEPSolver/CellModels/CellModel.h>
#include <CardiacEPSolver/Stimuli/Stimulus.h>
using namespace Nektar;
int main(int argc, char *argv[])
{
SpatialDomains::VertexComponentSharedPtr vPoint;
MultiRegions::ExpListSharedPtr vExp;
LibUtilities::SessionReaderSharedPtr vSession;
std::string vCellModel;
CellModelSharedPtr vCell;
std::vector<StimulusSharedPtr> vStimulus;
Array<OneD, Array<OneD, NekDouble> > vWsp(1);
Array<OneD, Array<OneD, NekDouble> > vSol(1);
NekDouble vDeltaT;
NekDouble vTime;
unsigned int nSteps;
// Create a session reader to read pacing parameters
vSession = LibUtilities::SessionReader::CreateInstance(argc, argv);
try
{
// Construct a field consisting of a single vertex
vPoint = MemoryManager<SpatialDomains::VertexComponent>
::AllocateSharedPtr(3, 0, 0.0, 0.0, 0.0);
vExp = MemoryManager<MultiRegions::ExpList0D>
::AllocateSharedPtr(vPoint);
// Get cell model name and create it
vSession->LoadSolverInfo("CELLMODEL", vCellModel, "");
ASSERTL0(vCellModel != "", "Cell Model not specified.");
vCell = GetCellModelFactory().CreateInstance(
vCellModel, vSession, vExp);
vCell->Initialise();
// Load the stimuli
vStimulus = Stimulus::LoadStimuli(vSession, vExp);
// Set up solution arrays, workspace and read in parameters
vSol[0] = Array<OneD, NekDouble>(1, 0.0);
vWsp[0] = Array<OneD, NekDouble>(1, 0.0);
vDeltaT = vSession->GetParameter("TimeStep");
vTime = 0.0;
nSteps = vSession->GetParameter("NumSteps");
vSol[0][0] = -81.0;
// Time integrate cell model
for (unsigned int i = 0; i < nSteps; ++i)
{
// Compute J_ion
vCell->TimeIntegrate(vSol, vWsp, vTime);
// Add stimuli J_stim
for (unsigned int i = 0; i < vStimulus.size(); ++i)
{
vStimulus[i]->Update(vWsp, vTime);
}
// Time-step with forward Euler
Vmath::Svtvp(1, vDeltaT, vWsp[0], 1, vSol[0], 1, vSol[0], 1);
// Increment time
vTime += vDeltaT;
// Output current solution to stdout
cout << vTime << " " << vSol[0][0] << endl;
}
}
catch (...)
{
cerr << "An error occured" << endl;
}
return 0;
}
<NEKTAR>
<CONDITIONS>
<PARAMETERS>
<P> TimeStep = 0.02 </P>
<P> FinTime = 15000 </P>
<P> NumSteps = FinTime/TimeStep </P>
<P> SubSteps = 1 </P>
</PARAMETERS>
<SOLVERINFO>
<I PROPERTY="CellModel" VALUE="CourtemancheRamirezNattel98" />
</SOLVERINFO>
</CONDITIONS>
<STIMULI>
<STIMULUS ID="0" TYPE="StimulusPoint">
<p_strength> 20.0 </p_strength>
<PROTOCOL TYPE = "ProtocolS1S2">
<START> 2.0 </START>
<DURATION> 2.0 </DURATION>
<S1CYCLELENGTH> 700.0 </S1CYCLELENGTH>
<NUM_S1> 50 </NUM_S1>
<S2CYCLELENGTH>0.0 </S2CYCLELENGTH>
</PROTOCOL>
</STIMULUS>
</STIMULI>
</NEKTAR>
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment