Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ExtendedBidomainProblem.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 "ExtendedBidomainProblem.hpp"
38#include "ExtendedBidomainSolver.hpp"
39#include "AbstractDynamicLinearPdeSolver.hpp"
40#include "HeartConfig.hpp"
41#include "Exception.hpp"
42#include "DistributedVector.hpp"
43#include "ReplicatableVector.hpp"
44
45template<unsigned DIM>
47 AbstractCardiacCellFactory<DIM>* pCellFactory, AbstractCardiacCellFactory<DIM>* pSecondCellFactory, bool hasBath)
48 : AbstractCardiacProblem<DIM,DIM, 3>(pCellFactory),
49 mpSecondCellFactory(pSecondCellFactory),
50 mpExtendedBidomainTissue(NULL),
51 mUserSpecifiedSecondCellConductivities(false),
52 mUserHasSetBidomainValuesExplicitly(false),
53 mpExtracellularStimulusFactory(NULL),
54 mRowForAverageOfPhiZeroed(INT_MAX),
55 mApplyAveragePhieZeroConstraintAfterSolving(false),
56 mUserSuppliedExtracellularStimulus(false),
57 mHasBath(hasBath)
58{
60}
61
62template<unsigned DIM>
64 : AbstractCardiacProblem<DIM,DIM, 3>(),
65 mpSecondCellFactory(NULL),
66 mpExtendedBidomainTissue(NULL),
67 mUserSpecifiedSecondCellConductivities(false),
68 mUserHasSetBidomainValuesExplicitly(false),
69 mpExtracellularStimulusFactory(NULL),
70 mRowForAverageOfPhiZeroed(INT_MAX),
71 mApplyAveragePhieZeroConstraintAfterSolving(false),
72 mUserSuppliedExtracellularStimulus(false)
73{
75}
76
77template<unsigned DIM>
79{
80
81 DistributedVectorFactory* p_factory = this->mpMesh->GetDistributedVectorFactory();
82 Vec initial_condition = p_factory->CreateVec(3);
83 DistributedVector ic = p_factory->CreateDistributedVector(initial_condition);
84 std::vector<DistributedVector::Stripe> stripe;
85 stripe.reserve(3);
86
87 stripe.push_back(DistributedVector::Stripe(ic, 0));
88 stripe.push_back(DistributedVector::Stripe(ic, 1));
89 stripe.push_back(DistributedVector::Stripe(ic, 2));
90
91 for (DistributedVector::Iterator index = ic.Begin();
92 index != ic.End();
93 ++index)
94 {
95 stripe[0][index] = mpExtendedBidomainTissue->GetCardiacCell(index.Global)->GetVoltage();//phi_i of frrst cell = Vm first cell at the start
96 stripe[1][index] = mpExtendedBidomainTissue->GetCardiacSecondCell(index.Global)->GetVoltage();//phi_i of second cell = Vm second cell at the start
97 stripe[2][index] = 0.0;//
98 }
99 ic.Restore();
100
101 return initial_condition;
102}
103
104template<unsigned DIM>
106{
107 if (mpExtracellularStimulusFactory == NULL) // user has not passed in any extracellular stimulus in any form
108 {
109 mpExtracellularStimulusFactory = new AbstractStimulusFactory<DIM>();
110 // Create one (with default implementation to zero stimulus everywhere)
111 }
112
113 assert(mpExtracellularStimulusFactory); // should be created by now, either above or by the user...
114 mpExtracellularStimulusFactory->SetMesh(this->mpMesh);//so, set the mesh into it.
115 mpExtracellularStimulusFactory->SetCompatibleExtracellularStimulus();//make sure compatibility condition will be valid
116
117 std::vector<AbstractChasteRegion<DIM>* > grounded_regions = mpExtracellularStimulusFactory->GetRegionsToBeGrounded();
118
119 if ((mUserSuppliedExtracellularStimulus) && grounded_regions.size() > 0) //we check for grunded nodes here
120 {
121 std::vector<unsigned> grounded_indices;
122 for (unsigned global_node_index = 0; global_node_index < this->mpMesh->GetNumNodes(); global_node_index++)
123 {
124 if (this->mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(global_node_index))
125 {
126 for (unsigned region_index = 0; region_index <grounded_regions.size(); region_index++)
127 {
128 if (grounded_regions[region_index]->DoesContain(this->mpMesh->GetNode(global_node_index)->GetPoint()))
129 {
130 grounded_indices.push_back( this->mpMesh->GetNode(global_node_index)->GetIndex() );
131 }
132 }
133 }
134 }
136 SetFixedExtracellularPotentialNodes(grounded_indices);
137 }
138}
139
140template<unsigned DIM>
142{
143 //set the mesh into the second cell factory as well.
144 mpSecondCellFactory->SetMesh(this->mpMesh);
145
146 //deal with extracellular stimulus, if any
147 ProcessExtracellularStimulus();
148
149 //Now create the tissue object
150 mpExtendedBidomainTissue = new ExtendedBidomainTissue<DIM>(this->mpCellFactory, mpSecondCellFactory,mpExtracellularStimulusFactory);
151
152 //Let the Tissue know if the user wants an extracellular stimulus (or if we had to create a default zero one).
153 mpExtendedBidomainTissue->SetUserSuppliedExtracellularStimulus(mUserSuppliedExtracellularStimulus);
154
155 //if the user remembered to set a different value for the sigma of the second cell...
156 if (mUserSpecifiedSecondCellConductivities)
157 {
158 mpExtendedBidomainTissue->SetIntracellularConductivitiesSecondCell(mIntracellularConductivitiesSecondCell);
159 }
160 else //..otherwise it gets the same as the first cell (according to heartconfig...)
161 {
162 c_vector<double, DIM> intra_conductivities;
164 mpExtendedBidomainTissue->SetIntracellularConductivitiesSecondCell(intra_conductivities);
165 }
166
167 //the conductivities for the first cell are created within the tissue constructor in the abstract class
168 //here we create the ones for the second cell
169 mpExtendedBidomainTissue->CreateIntracellularConductivityTensorSecondCell();
170
171 if (mUserHasSetBidomainValuesExplicitly)
172 {
173 mpExtendedBidomainTissue->SetAmFirstCell(mAmFirstCell);
174 mpExtendedBidomainTissue->SetAmSecondCell(mAmSecondCell);
175 mpExtendedBidomainTissue->SetAmGap(mAmGap);
176 mpExtendedBidomainTissue->SetGGap(mGGap);
177 mpExtendedBidomainTissue->SetCmFirstCell(mCmFirstCell);
178 mpExtendedBidomainTissue->SetCmSecondCell(mCmSecondCell);
179 }
180 else//we set all the Am and Cm to the values set by the heartconfig (only one value for all Am and one value for all Cms)
181 {
182 mpExtendedBidomainTissue->SetAmFirstCell(HeartConfig::Instance()->GetSurfaceAreaToVolumeRatio());
183 mpExtendedBidomainTissue->SetAmSecondCell(HeartConfig::Instance()->GetSurfaceAreaToVolumeRatio());
184 mpExtendedBidomainTissue->SetAmGap(HeartConfig::Instance()->GetSurfaceAreaToVolumeRatio());
185 mpExtendedBidomainTissue->SetGGap(0.0);
186 mpExtendedBidomainTissue->SetCmFirstCell(HeartConfig::Instance()->GetCapacitance());
187 mpExtendedBidomainTissue->SetCmSecondCell(HeartConfig::Instance()->GetCapacitance());
188 }
189
190 mpExtendedBidomainTissue->SetGgapHeterogeneities(mGgapHeterogeneityRegions, mGgapHeterogenousValues);//set user input into the tissue class
191 mpExtendedBidomainTissue->CreateGGapConductivities();//if vectors are empty, mGgap will be put everywhere by this method.
192
193 return mpExtendedBidomainTissue;
194}
195
196template<unsigned DIM>
197void ExtendedBidomainProblem<DIM>::SetExtendedBidomainParameters(double Am1, double Am2, double AmGap, double Cm1, double Cm2, double Ggap)
198{
199 mAmFirstCell = Am1;
200 mAmSecondCell = Am2;
201 mAmGap = AmGap;
202 mCmFirstCell = Cm1;
203 mCmSecondCell = Cm2;
204 mGGap = Ggap;
205
206 mUserHasSetBidomainValuesExplicitly = true;
207}
208
209template <unsigned DIM>
210void ExtendedBidomainProblem<DIM>::SetGgapHeterogeneities ( std::vector<boost::shared_ptr<AbstractChasteRegion<DIM> > >& rGgapHeterogeneityRegions, std::vector<double>& rGgapValues)
211{
212 if (rGgapHeterogeneityRegions.size() != rGgapValues.size() )
213 {
214 EXCEPTION ("Gap junction heterogeneity areas must be of the same number as the heterogeneity values");
215 }
216 mGgapHeterogeneityRegions = rGgapHeterogeneityRegions;
217 mGgapHeterogenousValues =rGgapValues;
218}
219
220template <unsigned DIM>
222{
223 mpExtracellularStimulusFactory = pFactory;
224 mUserSuppliedExtracellularStimulus = true;
225}
226
227template<unsigned DIM>
229{
230 /*
231 * NOTE: The this->mpBoundaryConditionsContainer.get() lines below convert a
232 * boost::shared_ptr to a normal pointer, as this is what the assemblers are
233 * expecting. We have to be a bit careful though as boost could decide to delete
234 * them whenever it feels like as it won't count the assemblers as using them.
235 *
236 * As long as they are kept as member variables here for as long as they are
237 * required in the assemblers it should all work OK.
238 */
239
240 mpSolver = new ExtendedBidomainSolver<DIM,DIM>( mHasBath,
241 this->mpMesh,
242 mpExtendedBidomainTissue,
243 this->mpBoundaryConditionsContainer.get());
244
245 try
246 {
247 mpSolver->SetFixedExtracellularPotentialNodes(mFixedExtracellularPotentialNodes);
248 mpSolver->SetRowForAverageOfPhiZeroed(mRowForAverageOfPhiZeroed);
249 }
250 catch (const Exception& e)
251 {
252 delete mpSolver;
253 throw e;
254 }
255
256 return mpSolver;
257}
258
259template<unsigned DIM>
261{
262 if (!mUserSuppliedExtracellularStimulus)
263 {
264 delete mpExtracellularStimulusFactory;
265 }
266}
267
268template<unsigned DIM>
270{
271 for (unsigned i = 0; i < DIM; i++)
272 {
273 mIntracellularConductivitiesSecondCell[i] = conductivities[i];
274 }
275 mUserSpecifiedSecondCellConductivities = true;
276}
277
278template<unsigned DIM>
280{
281 assert(mFixedExtracellularPotentialNodes.size() == 0);
282 mFixedExtracellularPotentialNodes.resize(nodes.size());
283 for (unsigned i=0; i<nodes.size(); i++)
284 {
285 // the assembler checks that the nodes[i] is less than
286 // the number of nodes in the mesh so this is not done here
287 mFixedExtracellularPotentialNodes[i] = nodes[i];
288 }
289}
290
291template<unsigned DIM>
293{
294 if (node==0)
295 {
296 mRowForAverageOfPhiZeroed = 2;
297 }
298 else
299 {
300 //Phie is every three lines, starting from zero...
301 mRowForAverageOfPhiZeroed = 3*node - 1;
302 }
303}
304
305template<unsigned DIM>
307{
308 assert(mpExtendedBidomainTissue!=NULL);
309 return mpExtendedBidomainTissue;
310}
311
312template<unsigned DIM>
314{
316 {
317 std::cout << "Solved to time " << time << "\n" << std::flush;
318 }
319
320 double V_max_first_cell, V_min_first_cell, V_max_second_cell, V_min_second_cell, phi_e_min, phi_e_max;
321
322 VecStrideMax( this->mSolution, 0, CHASTE_PETSC_NULLPTR, &V_max_first_cell );
323 VecStrideMin( this->mSolution, 0, CHASTE_PETSC_NULLPTR, &V_min_first_cell );
324
325 VecStrideMax( this->mSolution, 1, CHASTE_PETSC_NULLPTR, &V_max_second_cell );
326 VecStrideMin( this->mSolution, 1, CHASTE_PETSC_NULLPTR, &V_min_second_cell );
327
328 VecStrideMax( this->mSolution, 2, CHASTE_PETSC_NULLPTR, &phi_e_max );
329 VecStrideMin( this->mSolution, 2, CHASTE_PETSC_NULLPTR, &phi_e_min );
330
332 {
333 std::cout << " V first cell = " << "[" <<V_min_first_cell << ", " << V_max_first_cell << "]" << ";\n"
334 << " V second cell = " << "[" <<V_min_second_cell << ", " << V_max_second_cell << "]" << ";\n"
335 << " Phi_e = " << "[" <<phi_e_min << ", " << phi_e_max << "]" << ";\n"
336 << std::flush;
337 }
338}
339
340template<unsigned DIM>
342{
343 if (!extending)
344 {
345 if (this->mNodesToOutput.empty())
346 {
347 // Set writer to output all nodes
348 this->mpWriter->DefineFixedDimension(this->mpMesh->GetNumNodes());
349 }
350// else
351// {
352// // Output only the nodes indicted
353// this->mpWriter->DefineFixedDimension( this->mNodesToOutput, this->mpMesh->GetNumNodes() );
354// }
355 //mNodeColumnId = mpWriter->DefineVariable("Node", "dimensionless");
356 mVoltageColumnId_Vm1 = this->mpWriter->DefineVariable("V","mV");
357 mVoltageColumnId_Vm2 = this->mpWriter->DefineVariable("V_2","mV");
358 mVoltageColumnId_Phie = this->mpWriter->DefineVariable("Phi_e","mV");
359 mVariablesIDs.push_back(mVoltageColumnId_Vm1);
360 mVariablesIDs.push_back(mVoltageColumnId_Vm2);
361 mVariablesIDs.push_back(mVoltageColumnId_Phie);
362
363 // Only used to get an estimate of the # of timesteps below (copied from Abstract class)
365 HeartConfig::Instance()->GetSimulationDuration(),
366 HeartConfig::Instance()->GetPrintingTimeStep());
367 this->mpWriter->DefineUnlimitedDimension("Time", "msecs", stepper.EstimateTimeSteps()+1); // +1 for start and end
368 }
369 else
370 {
371 mVoltageColumnId_Vm1 = this->mpWriter->GetVariableByName("V");
372 mVoltageColumnId_Vm2 = this->mpWriter->GetVariableByName("V_2");
373 mVoltageColumnId_Phie = this->mpWriter->GetVariableByName("Phi_e");
374 }
375 //define any extra variable. NOTE: it must be in the first cell (not the second)
377
378}
379
380template<unsigned DIM>
382{
383 this->mpWriter->PutUnlimitedVariable(time);
384
385 // Create a striped vector
386 Vec ordered_voltages = this->mpMesh->GetDistributedVectorFactory()->CreateVec(3);
387 DistributedVector wrapped_ordered_solution = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(ordered_voltages);
388 DistributedVector wrapped_solution = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(voltageVec);
389
390 DistributedVector::Stripe V_first_cell_stripe(wrapped_solution,0);
391 DistributedVector::Stripe V_second_cell_stripe(wrapped_solution,1);
392 DistributedVector::Stripe phi_e_stripe(wrapped_solution,2);
393
394 DistributedVector::Stripe wrapped_ordered_solution_first_stripe(wrapped_ordered_solution,0);
395 DistributedVector::Stripe wrapped_ordered_solution_second_stripe(wrapped_ordered_solution,1);
396 DistributedVector::Stripe wrapped_ordered_solution_third_stripe(wrapped_ordered_solution,2);
397
398 for (DistributedVector::Iterator index = wrapped_solution.Begin();
399 index != wrapped_solution.End();
400 ++index)
401 {
402 wrapped_ordered_solution_first_stripe[index] = V_first_cell_stripe[index];
403 wrapped_ordered_solution_second_stripe[index] = V_second_cell_stripe[index];
404 wrapped_ordered_solution_third_stripe[index] = phi_e_stripe[index];
405 }
406 wrapped_solution.Restore();
407 wrapped_ordered_solution.Restore();
408
409 this->mpWriter->PutStripedVector(mVariablesIDs, ordered_voltages);
410 PetscTools::Destroy(ordered_voltages);
411 //write any extra variable. Note that this method in the parent class will
412 //take the extra variable only from the first cell.
415}
416
417template<unsigned DIM>
419{
421 if (mFixedExtracellularPotentialNodes.empty())
422 {
423 // We're not pinning any nodes.
424 if (mRowForAverageOfPhiZeroed==INT_MAX)
425 {
426 // We're not using the constrain Average phi_e to 0 method, hence use a null space
427 // Check that the KSP solver isn't going to do anything stupid:
428 // phi_e is not bounded, so it'd be wrong to use a relative tolerance
429 if (HeartConfig::Instance()->GetUseRelativeTolerance())
430 {
431 EXCEPTION("Bidomain external voltage is not bounded in this simulation - use KSP *absolute* tolerance");
432 }
433 }
434 }
435}
436
437template<unsigned DIM>
439{
440 return mHasBath;
441}
442
443template<unsigned DIM>
445{
446 mHasBath = hasBath;
447}
448
449// Explicit instantiation
450template class ExtendedBidomainProblem<1>;
451template class ExtendedBidomainProblem<2>;
452template class ExtendedBidomainProblem<3>;
453
454// 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)
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
void SetFixedExtracellularPotentialNodes(std::vector< unsigned > nodes)
virtual void DefineWriterColumns(bool extending)
void SetGgapHeterogeneities(std::vector< boost::shared_ptr< AbstractChasteRegion< DIM > > > &rGgapHeterogeneityRegions, std::vector< double > &rGgapValues)
void SetExtendedBidomainParameters(double Am1, double Am2, double AmGap, double Cm1, double Cm2, double Ggap)
void SetNodeForAverageOfPhiZeroed(unsigned node)
void SetIntracellularConductivitiesForSecondCell(c_vector< double, DIM > conductivities)
virtual AbstractCardiacTissue< DIM > * CreateCardiacTissue()
virtual void WriteOneStep(double time, Vec voltageVec)
void SetExtracellularStimulusFactory(AbstractStimulusFactory< DIM > *pFactory)
ExtendedBidomainTissue< DIM > * GetExtendedBidomainTissue()
virtual AbstractDynamicLinearPdeSolver< DIM, DIM, 3 > * CreateSolver()
std::vector< unsigned > mFixedExtracellularPotentialNodes
void GetIntracellularConductivities(c_vector< double, 3 > &rIntraConductivities) const
static HeartConfig * Instance()
static void Destroy(Vec &rVec)
static bool AmMaster()
static void Barrier(const std::string callerId="")
unsigned EstimateTimeSteps() const