37#include "BidomainProblem.hpp"
38#include "BidomainSolver.hpp"
39#include "HeartConfig.hpp"
41#include "DistributedVector.hpp"
42#include "ReplicatableVector.hpp"
44template <
unsigned DIM>
52 iter != this->mpMesh->GetNodeIteratorEnd();
58 bool any_bath_element_found =
false;
63 it != this->mpMesh->GetElementIteratorEnd();
78 any_bath_element_found =
true;
92 if (!(bath_identifiers.size()==1 && bath_identifiers.find(1)==bath_identifiers.begin()))
94 EXCEPTION(
"User has set bath identifiers, but the BidomainProblem isn't expecting a bath. Did you mean to use BidomainProblem(..., true)? Or alternatively, BidomainWithBathProblem(...)?");
106 DistributedVector ic = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(init_cond);
115 voltage_stripe[index] = 0.0;
124template<
unsigned DIM>
127 AnalyseMeshForBath();
129 return mpBidomainTissue;
132template<
unsigned DIM>
147 this->mpBoundaryConditionsContainer.get());
151 mpSolver->SetFixedExtracellularPotentialNodes(mFixedExtracellularPotentialNodes);
152 mpSolver->SetRowForAverageOfPhiZeroed(mRowForAverageOfPhiZeroed);
163template<
unsigned DIM>
167 mpBidomainTissue(NULL),
168 mRowForAverageOfPhiZeroed(INT_MAX),
174template<
unsigned DIM>
177 mpBidomainTissue(NULL),
178 mRowForAverageOfPhiZeroed(INT_MAX)
183template<
unsigned DIM>
186 mFixedExtracellularPotentialNodes.resize(nodes.size());
187 for (
unsigned i=0; i<nodes.size(); i++)
191 mFixedExtracellularPotentialNodes[i] = nodes[i];
195template<
unsigned DIM>
198 mRowForAverageOfPhiZeroed = 2*node+1;
201template<
unsigned DIM>
204 assert(mpBidomainTissue!=NULL);
205 return mpBidomainTissue;
208template<
unsigned DIM>
213 std::cout <<
"Solved to time " << time <<
"\n" << std::flush;
216 double v_max, v_min, phi_max, phi_min;
226 std::cout <<
" V; phi_e = " <<
"[" <<v_min <<
", " << v_max <<
"]" <<
";\t"
227 <<
"[" <<phi_min <<
", " << phi_max <<
"]" <<
"\n"
232template<
unsigned DIM>
238 mExtracelluarColumnId = this->mpWriter->GetVariableByName(
"Phi_e");
242 mExtracelluarColumnId = this->mpWriter->DefineVariable(
"Phi_e",
"mV");
247template<
unsigned DIM>
250 this->mpWriter->PutUnlimitedVariable(time);
251 std::vector<int> variable_ids;
252 variable_ids.push_back(this->mVoltageColumnId);
253 variable_ids.push_back(mExtracelluarColumnId);
254 this->mpWriter->PutStripedVector(variable_ids, voltageVec);
258template<
unsigned DIM>
262 if (mFixedExtracellularPotentialNodes.empty())
265 if (mRowForAverageOfPhiZeroed==INT_MAX)
272 EXCEPTION(
"Bidomain external voltage is not bounded in this simulation - use KSP *absolute* tolerance");
278template<
unsigned DIM>
287 assert(this->mpMesh!=NULL);
291 mpElectrodes = (boost::shared_ptr<Electrodes<DIM> >)
new Electrodes<DIM>(*(this->mpMesh));
295template<
unsigned DIM>
298 if (mpElectrodes && mpElectrodes->SwitchOn(time))
301 assert(this->mpBoundaryConditionsContainer);
306 mpSolver->ResetBoundaryConditionsContainer(mpElectrodes->GetBoundaryConditionsContainer().get());
310 this->mpBoundaryConditionsContainer = mpElectrodes->GetBoundaryConditionsContainer();
313 this->mpDefaultBoundaryConditionsContainer = this->mpBoundaryConditionsContainer;
317 if (mpSolver->GetLinearSystem() != NULL)
321 if (mpElectrodes->HasGroundedElectrode())
323 this->mpBoundaryConditionsContainer->ApplyDirichletToLinearProblem( *(mpSolver->GetLinearSystem()),
328 if (mpElectrodes->HasGroundedElectrode())
330 mpSolver->GetLinearSystem()->RemoveNullSpace();
336template<
unsigned DIM>
339 if (mpElectrodes && mpElectrodes->SwitchOff(time))
343 assert(this->mpBoundaryConditionsContainer);
349 for (
unsigned problem_index=0; problem_index<2; problem_index++)
351 this->mpDefaultBoundaryConditionsContainer->DefineZeroNeumannOnMeshBoundary(this->mpMesh, problem_index);
356 if (mpElectrodes->HasGroundedElectrode())
365 mpSolver->ResetBoundaryConditionsContainer(this->mpDefaultBoundaryConditionsContainer.get());
368 this->mpBoundaryConditionsContainer = this->mpDefaultBoundaryConditionsContainer;
372template<
unsigned DIM>
377 rAdditionalStoppingTimes.push_back(mpElectrodes->GetSwitchOnTime());
378 rAdditionalStoppingTimes.push_back(mpElectrodes->GetSwitchOffTime());
382template<
unsigned DIM>
#define EXCEPTION(message)
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
virtual Vec CreateInitialCondition()
void DefineExtraVariablesWriterColumns(bool extending)
virtual void PreSolveChecks()
void WriteExtraVariablesOneStep()
virtual void DefineWriterColumns(bool extending)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
unsigned GetUnsignedAttribute()
unsigned GetNumNodes() const
void SetNodeForAverageOfPhiZeroed(unsigned node)
void AnalyseMeshForBath()
void SetUpAdditionalStoppingTimes(std::vector< double > &rAdditionalStoppingTimes)
void AtBeginningOfTimestep(double time)
virtual AbstractDynamicLinearPdeSolver< DIM, DIM, 2 > * CreateSolver()
Vec CreateInitialCondition()
virtual void DefineWriterColumns(bool extending)
std::vector< unsigned > mFixedExtracellularPotentialNodes
virtual AbstractCardiacTissue< DIM > * CreateCardiacTissue()
void OnEndOfTimestep(double time)
void SetFixedExtracellularPotentialNodes(std::vector< unsigned > nodes)
void WriteInfo(double time)
BidomainTissue< DIM > * GetBidomainTissue()
virtual void WriteOneStep(double time, Vec voltageVec)
const std::set< unsigned > & rGetBathIdentifiers()
static HeartConfig * Instance()
static bool IsRegionTissue(HeartRegionType regionId)
static HeartRegionType GetValidTissueId()
static bool IsRegionBath(HeartRegionType regionId)
static HeartRegionType GetValidBathId()