37 #include "BidomainProblem.hpp"
38 #include "BidomainSolver.hpp"
39 #include "HeartConfig.hpp"
41 #include "DistributedVector.hpp"
42 #include "ReplicatableVector.hpp"
44 template <
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(...)?");
99 template<
unsigned DIM>
106 DistributedVector ic = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(init_cond);
115 voltage_stripe[index] = 0.0;
124 template<
unsigned DIM>
127 AnalyseMeshForBath();
129 return mpBidomainTissue;
132 template<
unsigned DIM>
147 this->mpBoundaryConditionsContainer.get());
152 mpSolver->SetRowForAverageOfPhiZeroed(mRowForAverageOfPhiZeroed);
163 template<
unsigned DIM>
167 mpBidomainTissue(NULL),
168 mRowForAverageOfPhiZeroed(INT_MAX),
174 template<
unsigned DIM>
177 mpBidomainTissue(NULL),
178 mRowForAverageOfPhiZeroed(INT_MAX)
183 template<
unsigned DIM>
186 mFixedExtracellularPotentialNodes.resize(nodes.size());
187 for (
unsigned i=0; i<nodes.size(); i++)
191 mFixedExtracellularPotentialNodes[i] = nodes[i];
195 template<
unsigned DIM>
198 mRowForAverageOfPhiZeroed = 2*node+1;
201 template<
unsigned DIM>
204 assert(mpBidomainTissue!=NULL);
205 return mpBidomainTissue;
208 template<
unsigned DIM>
213 std::cout <<
"Solved to time " << time <<
"\n" << std::flush;
216 double v_max, v_min, phi_max, phi_min;
218 VecStrideMax( this->mSolution, 0, PETSC_NULL, &v_max );
219 VecStrideMin( this->mSolution, 0, PETSC_NULL, &v_min );
221 VecStrideMax( this->mSolution, 1, PETSC_NULL, &phi_max );
222 VecStrideMin( this->mSolution, 1, PETSC_NULL, &phi_min );
226 std::cout <<
" V; phi_e = " <<
"[" <<v_min <<
", " << v_max <<
"]" <<
";\t"
227 <<
"[" <<phi_min <<
", " << phi_max <<
"]" <<
"\n"
232 template<
unsigned DIM>
238 mExtracelluarColumnId = this->mpWriter->GetVariableByName(
"Phi_e");
242 mExtracelluarColumnId = this->mpWriter->DefineVariable(
"Phi_e",
"mV");
247 template<
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);
258 template<
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");
278 template<
unsigned DIM>
287 assert(this->mpMesh!=NULL);
291 mpElectrodes = (boost::shared_ptr<Electrodes<DIM> >)
new Electrodes<DIM>(*(this->mpMesh));
295 template<
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();
336 template<
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;
372 template<
unsigned DIM>
377 rAdditionalStoppingTimes.push_back(mpElectrodes->GetSwitchOnTime());
378 rAdditionalStoppingTimes.push_back(mpElectrodes->GetSwitchOffTime());
382 template<
unsigned DIM>
void AnalyseMeshForBath()
ElementIterator GetElementIteratorBegin(bool skipDeletedElements=true)
static bool IsRegionBath(HeartRegionType regionId)
#define EXCEPTION(message)
static HeartRegionType GetValidBathId()
virtual Vec CreateInitialCondition()
BidomainTissue< DIM > * GetBidomainTissue()
void SetUpAdditionalStoppingTimes(std::vector< double > &rAdditionalStoppingTimes)
void SetFixedExtracellularPotentialNodes(std::vector< unsigned > nodes)
std::vector< unsigned > mFixedExtracellularPotentialNodes
virtual void PreSolveChecks()
void WriteInfo(double time)
void SetFixedExtracellularPotentialNodes(std::vector< unsigned > fixedExtracellularPotentialNodes)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
virtual void DefineWriterColumns(bool extending)
bool GetUseStateVariableInterpolation() const
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
Vec CreateInitialCondition()
unsigned GetNumNodes() const
void OnEndOfTimestep(double time)
virtual void WriteOneStep(double time, Vec voltageVec)
NodeIterator GetNodeIteratorBegin(bool skipDeletedNodes=true)
virtual void DefineWriterColumns(bool extending)
unsigned GetUnsignedAttribute()
const std::set< unsigned > & rGetBathIdentifiers()
void SetNodeForAverageOfPhiZeroed(unsigned node)
void AtBeginningOfTimestep(double time)
static HeartRegionType GetValidTissueId()
void WriteExtraVariablesOneStep()
static HeartConfig * Instance()
virtual AbstractCardiacTissue< DIM > * CreateCardiacTissue()
bool GetUseRelativeTolerance() const
virtual AbstractDynamicLinearPdeSolver< DIM, DIM, 2 > * CreateSolver()
void DefineExtraVariablesWriterColumns(bool extending)
static bool IsRegionTissue(HeartRegionType regionId)