64 "Bidomain model of cardiac electrophysiology with 3D diffusion.");
78 UnsteadySystem::v_InitObject();
82 std::string vCellModel;
83 m_session->LoadSolverInfo(
"CELLMODEL", vCellModel,
"");
85 ASSERTL0(vCellModel !=
"",
"Cell Model not specified.");
97 std::string varName[3] = {
98 "AnisotropicConductivityX",
99 "AnisotropicConductivityY",
100 "AnisotropicConductivityZ"
104 if (
m_session->DefinesFunction(
"IntracellularConductivity") &&
m_session->DefinesFunction(
"ExtracellularConductivity"))
123 =
m_session->GetFunction(
"IntracellularConductivity", varName[i]);
125 =
m_session->GetFunction(
"ExtracellularConductivity", varName[i]);
126 for(
int j = 0; j < nq; j++)
128 tmp1[i][j] = ifunc1->Evaluate(x0[j],x1[j],x2[j],0.0);
129 tmp2[i][j] = ifunc2->Evaluate(x0[j],x1[j],x2[j],0.0);
131 Vmath::Vadd(nq, tmp1[i], 1, tmp2[i], 1, tmp3[i], 1);
138 if (
m_session->DefinesParameter(
"StimulusDuration"))
141 "Stimulus function not defined.");
154 LibUtilities::FilterMap::const_iterator x;
155 for (x = f.begin(); x != f.end(); ++x, ++k)
157 if (x->first ==
"CheckpointCellModel")
159 boost::shared_ptr<FilterCheckpointCellModel> c
195 int nvariables = inarray.num_elements();
205 for (
int i = 0; i < nvariables; ++i)
209 Vmath::Vcopy(nq, &inarray[i][0], 1, &outarray[i][0], 1);
217 m_fields[i]->PhysDeriv(inarray[1],ggrad0);
219 m_fields[i]->PhysDeriv(0,ggrad0,ggrad0);
221 if (
m_session->DefinesFunction(
"IntracellularConductivity") &&
m_session->DefinesFunction(
"ExtracellularConductivity"))
237 outarray[i] =
m_fields[i]->GetPhys();
242 m_fields[i]->PhysDeriv(inarray[1],ggrad0,ggrad1);
244 m_fields[i]->PhysDeriv(0,ggrad0,ggrad0);
245 m_fields[i]->PhysDeriv(1,ggrad1,ggrad1);
247 if (
m_session->DefinesFunction(
"IntracellularConductivity") &&
m_session->DefinesFunction(
"ExtracellularConductivity"))
264 outarray[i] =
m_fields[i]->GetPhys();
269 m_fields[i]->PhysDeriv(inarray[1],ggrad0,ggrad1,ggrad2);
271 m_fields[i]->PhysDeriv(0,ggrad0,ggrad0);
272 m_fields[i]->PhysDeriv(1,ggrad1,ggrad1);
273 m_fields[i]->PhysDeriv(2,ggrad2,ggrad2);
275 if (
m_session->DefinesFunction(
"IntracellularConductivity") &&
m_session->DefinesFunction(
"ExtracellularConductivity"))
294 outarray[i] =
m_fields[i]->GetPhys();
305 m_fields[i]->PhysDeriv(0,grad0,grad0);
307 if (
m_session->DefinesFunction(
"IntracellularConductivity") &&
m_session->DefinesFunction(
"ExtracellularConductivity"))
320 outarray[i] =
m_fields[i]->GetPhys();
327 m_fields[i]->PhysDeriv(0,grad0,grad0);
328 m_fields[i]->PhysDeriv(1,grad1,grad1);
330 if (
m_session->DefinesFunction(
"IntracellularConductivity") &&
m_session->DefinesFunction(
"ExtracellularConductivity"))
344 outarray[i] =
m_fields[i]->GetPhys();
351 m_fields[i]->PhysDeriv(0,grad0,grad0);
352 m_fields[i]->PhysDeriv(1,grad1,grad1);
353 m_fields[i]->PhysDeriv(2,grad2,grad2);
355 if (
m_session->DefinesFunction(
"IntracellularConductivity") &&
m_session->DefinesFunction(
"ExtracellularConductivity"))
371 outarray[i] =
m_fields[i]->GetPhys();
385 m_cell->TimeIntegrate(inarray, outarray, time);
397 =
m_session->GetFunction(
"Stimulus",
"u");
398 ifunc->Evaluate(x0,x1,x2,time, result);
400 Vmath::Vadd(nq, outarray[0], 1, result, 1, outarray[0], 1);
407 bool dumpInitialConditions,
410 EquationSystem::v_SetInitialConditions(initialtime, dumpInitialConditions, domain);
419 UnsteadySystem::v_GenerateSummary(s);
422 ASSERTL0(
false,
"Update the generate summary");
468 m_cell->GenerateSummary(s);
#define ASSERTL0(condition, msg)
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
bool m_explicitDiffusion
Indicates if explicit or implicit treatment of diffusion is used.
void DefineImplicitSolve(FuncPointerT func, ObjectPointerT obj)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Computes the reaction terms and .
void SetCellModel(CellModelSharedPtr &pCellModel)
std::map< ConstFactorType, NekDouble > ConstFactorMap
std::vector< std::pair< std::string, std::string > > SummaryList
virtual ~Bidomain()
Desctructor.
static std::string className
Name of class.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Bidomain(const LibUtilities::SessionReaderSharedPtr &pSession)
Constructor.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
StdRegions::VarCoeffMap m_vardiffie
NekDouble m_stimDuration
Stimulus current.
virtual void v_GenerateSummary(SummaryList &s)
Prints a summary of the model parameters.
Base class for unsteady solvers.
std::vector< std::pair< std::string, FilterParams > > FilterMap
int m_spacedim
Spatial dimension (>= expansion dim).
virtual void v_SetInitialConditions(NekDouble initialtime, bool dumpInitialConditions, const int domain)
Sets a custom initial condition.
boost::shared_ptr< Equation > EquationSharedPtr
StdRegions::VarCoeffMap m_vardiffi
EquationSystemFactory & GetEquationSystemFactory()
CellModelFactory & GetCellModelFactory()
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
CellModelSharedPtr m_cell
Cell model.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
std::vector< FilterSharedPtr > m_filters
virtual void v_InitObject()
Init object for UnsteadySystem class.
Array< OneD, Array< OneD, NekDouble > > tmp1
Array< OneD, Array< OneD, NekDouble > > tmp2
Array< OneD, Array< OneD, NekDouble > > tmp3
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void DoImplicitSolve(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, NekDouble time, NekDouble lambda)
Solve for the diffusion term.
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
static FlagList NullFlagList
An empty flag list.
std::vector< int > m_intVariables
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.