Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
BidomainProblem.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36
37#include "BidomainProblem.hpp"
38#include "BidomainSolver.hpp"
39#include "HeartConfig.hpp"
40#include "Exception.hpp"
41#include "DistributedVector.hpp"
42#include "ReplicatableVector.hpp"
43
44template <unsigned DIM>
46{
47 // Annotate bath notes with correct region code
48 if (mHasBath)
49 {
50 // Initialize all nodes to be bath nodes
51 for (typename AbstractTetrahedralMesh<DIM,DIM>::NodeIterator iter=this->mpMesh->GetNodeIteratorBegin();
52 iter != this->mpMesh->GetNodeIteratorEnd();
53 ++iter)
54 {
55 (*iter).SetRegion(HeartRegionCode::GetValidBathId());
56 }
57
58 bool any_bath_element_found = false;
59
60 // Set nodes that are part of a heart element to be heart nodes
61 //for (unsigned i=0; i<this->mpMesh->GetNumElements(); i++)
62 for (typename AbstractTetrahedralMesh<DIM,DIM>::ElementIterator it = this->mpMesh->GetElementIteratorBegin();
63 it != this->mpMesh->GetElementIteratorEnd();
64 ++it)
65 {
66 Element<DIM, DIM>& r_element = *it;
67
69 {
70 for (unsigned j=0; j<r_element.GetNumNodes(); j++)
71 {
72 r_element.GetNode(j)->SetRegion(HeartRegionCode::GetValidTissueId());
73 }
74 }
75 else
76 {
78 any_bath_element_found = true;
79 }
80 }
81
82 if (!PetscTools::ReplicateBool(any_bath_element_found))
83 {
84 EXCEPTION("No bath element found");
85 }
86 }
87 else
88 {
89 // The problem hasn't been set up with a bath, so check that the user hasn't set any options
90 // which would suggest they're expecting there to be a bath!
91 std::set<unsigned> bath_identifiers = HeartConfig::Instance()->rGetBathIdentifiers();
92 if (!(bath_identifiers.size()==1 && bath_identifiers.find(1)==bath_identifiers.begin())) // IF NOT only 1 in the bath identifiers set
93 {
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(...)?");
95 }
96 }
97}
98
99template<unsigned DIM>
101{
103 if (mHasBath)
104 {
105 // get the voltage stripe
106 DistributedVector ic = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(init_cond);
108
109 for (DistributedVector::Iterator index = ic.Begin();
110 index!= ic.End();
111 ++index)
112 {
113 if (HeartRegionCode::IsRegionBath( this->mpMesh->GetNode( index.Global )->GetRegion() ))
114 {
115 voltage_stripe[index] = 0.0;
116 }
117 }
118 ic.Restore();
119 }
120
121 return init_cond;
122}
123
124template<unsigned DIM>
126{
127 AnalyseMeshForBath();
128 mpBidomainTissue = new BidomainTissue<DIM>(this->mpCellFactory, HeartConfig::Instance()->GetUseStateVariableInterpolation());
129 return mpBidomainTissue;
130}
131
132template<unsigned DIM>
134{
135 /*
136 * NOTE: The this->mpBoundaryConditionsContainer.get() lines below convert a
137 * boost::shared_ptr to a normal pointer, as this is what the solvers are
138 * expecting. We have to be a bit careful though as boost could decide to delete
139 * them whenever it feels like as it won't count the assembers as using them.
140 *
141 * As long as they are kept as member variables here for as long as they are
142 * required in the solvers it should all work OK.
143 */
144 mpSolver = new BidomainSolver<DIM,DIM>(mHasBath,
145 this->mpMesh,
146 mpBidomainTissue,
147 this->mpBoundaryConditionsContainer.get());
148
149 try
150 {
151 mpSolver->SetFixedExtracellularPotentialNodes(mFixedExtracellularPotentialNodes);
152 mpSolver->SetRowForAverageOfPhiZeroed(mRowForAverageOfPhiZeroed);
153 }
154 catch (const Exception& e)
155 {
156 delete mpSolver;
157 throw e;
158 }
159
160 return mpSolver;
161}
162
163template<unsigned DIM>
165 AbstractCardiacCellFactory<DIM>* pCellFactory, bool hasBath)
166 : AbstractCardiacProblem<DIM,DIM, 2>(pCellFactory),
167 mpBidomainTissue(NULL),
168 mRowForAverageOfPhiZeroed(INT_MAX),
169 mHasBath(hasBath)
170{
172}
173
174template<unsigned DIM>
176 : AbstractCardiacProblem<DIM, DIM, 2>(),
177 mpBidomainTissue(NULL),
178 mRowForAverageOfPhiZeroed(INT_MAX)
179{
181}
182
183template<unsigned DIM>
185{
186 mFixedExtracellularPotentialNodes.resize(nodes.size());
187 for (unsigned i=0; i<nodes.size(); i++)
188 {
189 // the solver checks that the nodes[i] is less than
190 // the number of nodes in the mesh so this is not done here
191 mFixedExtracellularPotentialNodes[i] = nodes[i];
192 }
193}
194
195template<unsigned DIM>
197{
198 mRowForAverageOfPhiZeroed = 2*node+1;
199}
200
201template<unsigned DIM>
203{
204 assert(mpBidomainTissue!=NULL);
205 return mpBidomainTissue;
206}
207
208template<unsigned DIM>
210{
212 {
213 std::cout << "Solved to time " << time << "\n" << std::flush;
214 }
215
216 double v_max, v_min, phi_max, phi_min;
217
218 VecStrideMax( this->mSolution, 0, CHASTE_PETSC_NULLPTR, &v_max );
219 VecStrideMin( this->mSolution, 0, CHASTE_PETSC_NULLPTR, &v_min );
220
221 VecStrideMax( this->mSolution, 1, CHASTE_PETSC_NULLPTR, &phi_max );
222 VecStrideMin( this->mSolution, 1, CHASTE_PETSC_NULLPTR, &phi_min );
223
225 {
226 std::cout << " V; phi_e = " << "[" <<v_min << ", " << v_max << "]" << ";\t"
227 << "[" <<phi_min << ", " << phi_max << "]" << "\n"
228 << std::flush;
229 }
230}
231
232template<unsigned DIM>
234{
236 if (extending)
237 {
238 mExtracelluarColumnId = this->mpWriter->GetVariableByName("Phi_e");
239 }
240 else
241 {
242 mExtracelluarColumnId = this->mpWriter->DefineVariable("Phi_e","mV");
243 }
245}
246
247template<unsigned DIM>
248void BidomainProblem<DIM>::WriteOneStep(double time, Vec voltageVec)
249{
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);
256}
257
258template<unsigned DIM>
260{
262 if (mFixedExtracellularPotentialNodes.empty())
263 {
264 // We're not pinning any nodes.
265 if (mRowForAverageOfPhiZeroed==INT_MAX)
266 {
267 // We're not using the constrain Average phi_e to 0 method, hence use a null space
268 // Check that the KSP solver isn't going to do anything stupid:
269 // phi_e is not bounded, so it'd be wrong to use a relative tolerance
270 if (HeartConfig::Instance()->GetUseRelativeTolerance())
271 {
272 EXCEPTION("Bidomain external voltage is not bounded in this simulation - use KSP *absolute* tolerance");
273 }
274 }
275 }
276}
277
278template<unsigned DIM>
280{
281 if (!mHasBath)
282 {
283 //Cannot set electrodes when problem has been defined to not have a bath
284 return;
285 }
286
287 assert(this->mpMesh!=NULL);
288
289 if (HeartConfig::Instance()->IsElectrodesPresent())
290 {
291 mpElectrodes = (boost::shared_ptr<Electrodes<DIM> >) new Electrodes<DIM>(*(this->mpMesh));
292 }
293}
294
295template<unsigned DIM>
297{
298 if (mpElectrodes && mpElectrodes->SwitchOn(time))
299 {
300 // At the moment mpBcc and mpDefaultBcc point to a set default BC
301 assert(this->mpBoundaryConditionsContainer);
302 //assert(this->mpDefaultBoundaryConditionsContainer);
303
304 // Note, no point calling this->SetBoundaryConditionsContainer() as the
305 // solver has already been created..
306 mpSolver->ResetBoundaryConditionsContainer(mpElectrodes->GetBoundaryConditionsContainer().get());
307
308 // ..but we set mpBcc anyway, so the local mpBcc is
309 // the same as the one being used in the solver...
310 this->mpBoundaryConditionsContainer = mpElectrodes->GetBoundaryConditionsContainer();
311
313 this->mpDefaultBoundaryConditionsContainer = this->mpBoundaryConditionsContainer;
314
315 // At t==0 or after checkpointing we won't have a system assembled at this stage: BCs will be applied once the matrix
316 // is assembled. Dirichlet BCs will be present at the time of assembly and no null space will be created either.
317 if (mpSolver->GetLinearSystem() != NULL)
318 {
319 // System matrix is assembled once at the beginning of the simulation. After that, nobody will take care
320 // of applying new BC to the system matrix. Must be triggered explicitly.
321 if (mpElectrodes->HasGroundedElectrode())
322 {
323 this->mpBoundaryConditionsContainer->ApplyDirichletToLinearProblem( *(mpSolver->GetLinearSystem()),
324 true, false);
325 }
326
327 // If a grounded electrode is switched on, the linear system is not singular anymore. Remove the null space.
328 if (mpElectrodes->HasGroundedElectrode())
329 {
330 mpSolver->GetLinearSystem()->RemoveNullSpace();
331 }
332 }
333 }
334}
335
336template<unsigned DIM>
338{
339 if (mpElectrodes && mpElectrodes->SwitchOff(time))
340 {
341 // At the moment mpBcc should exist and therefore
342 // mpDefaultBcc should be empty (not if electrodes switched on after 0ms)
343 assert(this->mpBoundaryConditionsContainer);
344 //assert(! this->mpDefaultBoundaryConditionsContainer);
345
346 // Set up default boundary conditions container - no Neumann fluxes
347 // or Dirichlet fixed nodes
348 this->mpDefaultBoundaryConditionsContainer.reset(new BoundaryConditionsContainer<DIM,DIM,2>);
349 for (unsigned problem_index=0; problem_index<2; problem_index++)
350 {
351 this->mpDefaultBoundaryConditionsContainer->DefineZeroNeumannOnMeshBoundary(this->mpMesh, problem_index);
352 }
353
354 // If there's a grounded electrode, we must remove BC from linear system. At the moment, we don't
355 // have a sensible way of doing this, therefore we reassemble the system.
356 if (mpElectrodes->HasGroundedElectrode())
357 {
358 delete mpSolver;
360 mpSolver->SetTimeStep(HeartConfig::Instance()->GetPdeTimeStep());
361 }
362
363 // Note, no point calling this->SetBoundaryConditionsContainer() as the
364 // solver has already been created..
365 mpSolver->ResetBoundaryConditionsContainer(this->mpDefaultBoundaryConditionsContainer.get());
366 // ..but we set mpBcc to be mpDefaultBcc anyway, so the local mpBcc is
367 // the same as the one being used in the solver...
368 this->mpBoundaryConditionsContainer = this->mpDefaultBoundaryConditionsContainer;
369 }
370}
371
372template<unsigned DIM>
373void BidomainProblem<DIM>::SetUpAdditionalStoppingTimes(std::vector<double>& rAdditionalStoppingTimes)
374{
375 if (mpElectrodes)
376 {
377 rAdditionalStoppingTimes.push_back(mpElectrodes->GetSwitchOnTime());
378 rAdditionalStoppingTimes.push_back(mpElectrodes->GetSwitchOffTime());
379 }
380}
381
382template<unsigned DIM>
384{
385 return mHasBath;
386}
387
388// Explicit instantiation
389template class BidomainProblem<1>;
390template class BidomainProblem<2>;
391template class BidomainProblem<3>;
392
393
394// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define CHASTE_PETSC_NULLPTR
A macro to define PETSc null pointer based on the PETSc version.
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
void DefineExtraVariablesWriterColumns(bool extending)
virtual void DefineWriterColumns(bool extending)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
unsigned GetUnsignedAttribute()
unsigned GetNumNodes() const
void SetNodeForAverageOfPhiZeroed(unsigned node)
void SetUpAdditionalStoppingTimes(std::vector< double > &rAdditionalStoppingTimes)
void AtBeginningOfTimestep(double time)
virtual AbstractDynamicLinearPdeSolver< DIM, DIM, 2 > * CreateSolver()
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()
static bool ReplicateBool(bool flag)
static bool AmMaster()