Chaste Commit::ca8ccdedf819b6e02855bc0e8e6f50bdecbc5208
ExtendedBidomainTissue.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#include "ExtendedBidomainTissue.hpp"
37
38#include "DistributedVector.hpp"
39#include "OrthotropicConductivityTensors.hpp"
40#include "AxisymmetricConductivityTensors.hpp"
41#include "AbstractStimulusFunction.hpp"
42#include "ChastePoint.hpp"
43#include "AbstractChasteRegion.hpp"
44#include "HeartEventHandler.hpp"
45
46
55template <unsigned SPACE_DIM>
56void ExplicitCVectorCopy(c_vector<double, SPACE_DIM>& rFrom, c_vector<double, SPACE_DIM>& rTo)
57{
58 if constexpr (SPACE_DIM == 1u)
59 {
60 // In 1-D, explicitly copy the data
61 rTo[0] = rFrom[0];
62 }
63 else
64 {
65 // By default, use the standard copy-assignment overload
66 rTo = rFrom;
67 }
68}
69
70
71template <unsigned SPACE_DIM>
73 AbstractCardiacCellFactory<SPACE_DIM>* pCellFactorySecondCell,
74 AbstractStimulusFactory<SPACE_DIM>* pExtracellularStimulusFactory)
75 : AbstractCardiacTissue<SPACE_DIM>(pCellFactory),
76 mpIntracellularConductivityTensorsSecondCell(NULL),
77 mUserSuppliedExtracellularStimulus(false)
78{
79 //First, do the same that the abstract constructor does, but applied to the second cell
80
81 assert(pCellFactorySecondCell != NULL);
82 assert(pCellFactorySecondCell->GetMesh() != NULL);
83 assert(pCellFactorySecondCell->GetNumberOfCells() == pCellFactory->GetNumberOfCells() );
84 assert(pExtracellularStimulusFactory != NULL);
85 assert(pExtracellularStimulusFactory->GetMesh() != NULL);
86 assert(pExtracellularStimulusFactory->GetNumberOfCells() == pCellFactorySecondCell->GetNumberOfCells() );
87
88 unsigned num_local_nodes = this->mpDistributedVectorFactory->GetLocalOwnership();
89 unsigned ownership_range_low = this->mpDistributedVectorFactory->GetLow();
90 mCellsDistributedSecondCell.resize(num_local_nodes);
91 mGgapDistributed.resize(num_local_nodes);
92 mExtracellularStimuliDistributed.resize(num_local_nodes);
93
94 try
95 {
96 for (unsigned local_index = 0; local_index < num_local_nodes; local_index++)
97 {
98 unsigned global_index = local_index + ownership_range_low;
99 Node<SPACE_DIM>* p_node = this->mpMesh->GetNode(global_index);
100 mCellsDistributedSecondCell[local_index] = pCellFactorySecondCell->CreateCardiacCellForNode(p_node);
102 mGgapDistributed[local_index] = 0.0;//default. It will be changed by specific method later when user input will be obtained
103 mExtracellularStimuliDistributed[local_index] = pExtracellularStimulusFactory->CreateStimulusForNode(p_node);
104 }
105
106 pCellFactorySecondCell->FinaliseCellCreation(&mCellsDistributedSecondCell,
108 this->mpDistributedVectorFactory->GetHigh());
109 }
110 // LCOV_EXCL_START //don't really know how to cover this...
111 catch (const Exception& e)
112 {
113 // Errors thrown creating cells will often be process-specific
115 // Should really do this for other processes too, but this is all we need
116 // to get memory testing to pass, and leaking when we're about to die isn't
117 // that bad! Delete second cells
118 for (std::vector<AbstractCardiacCellInterface*>::iterator cell_iterator = mCellsDistributedSecondCell.begin();
119 cell_iterator != mCellsDistributedSecondCell.end();
120 ++cell_iterator)
121 {
122 delete (*cell_iterator);
123 }
124 throw e;
125 }
126 // LCOV_EXCL_STOP
128
129 HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
132 mGgapCacheReplicated.Resize(pCellFactorySecondCell->GetNumberOfCells());//this is a bit of a hack...
133 mExtracellularStimulusCacheReplicated.Resize(pExtracellularStimulusFactory->GetNumberOfCells());
134 HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
135
136 //Create the extracellular conductivity tensor
138}
139
140//archiving constructor
141template <unsigned SPACE_DIM>
142ExtendedBidomainTissue<SPACE_DIM>::ExtendedBidomainTissue(std::vector<AbstractCardiacCellInterface*> & rCellsDistributed,
143 std::vector<AbstractCardiacCellInterface*> & rSecondCellsDistributed,
144 std::vector<boost::shared_ptr<AbstractStimulusFunction> > & rExtraStimuliDistributed,
145 std::vector<double> & rGgapsDistributed,
147 c_vector<double, SPACE_DIM> intracellularConductivitiesSecondCell)
148 : AbstractCardiacTissue<SPACE_DIM>(pMesh),
149 mpIntracellularConductivityTensorsSecondCell(NULL),
150 mIntracellularConductivitiesSecondCell(intracellularConductivitiesSecondCell),
151 mCellsDistributedSecondCell(rSecondCellsDistributed),
152 mExtracellularStimuliDistributed(rExtraStimuliDistributed),
153 mGgapDistributed(rGgapsDistributed),
154 mUserSuppliedExtracellularStimulus(false)
155{
156 //segfault guards in case we failed to load anything from the archive
157 assert(mCellsDistributedSecondCell.size() > 0);
158 assert(mExtracellularStimuliDistributed.size() > 0);
159 assert(mGgapDistributed.size() > 0);
160 //allocate memory for the caches
165
168}
169
170
171template <unsigned SPACE_DIM>
172void ExtendedBidomainTissue<SPACE_DIM>::SetGgapHeterogeneities(std::vector<boost::shared_ptr<AbstractChasteRegion<SPACE_DIM> > >& rGgapHeterogeneityRegions,
173 std::vector<double> rGgapValues)
174{
175 assert( rGgapHeterogeneityRegions.size() == rGgapValues.size() );//problem class (which calls this method should have thrown otherwise)
176 mGgapHeterogeneityRegions = rGgapHeterogeneityRegions;
177 mGgapValues =rGgapValues;
178}
179
180template <unsigned SPACE_DIM>
182{
183 assert(mGgapHeterogeneityRegions.size() == mGgapValues.size());
184 assert(this->mpMesh != NULL);
185
186 unsigned ownership_range_low = this->mpDistributedVectorFactory->GetLow();
187 unsigned num_local_nodes = this->mpDistributedVectorFactory->GetLocalOwnership();
188 assert(mGgapDistributed.size() == num_local_nodes);//the constructor should have allocated memory.
189 try
190 {
191 for (unsigned local_index = 0; local_index < num_local_nodes; local_index++)
192 {
193 unsigned global_index = ownership_range_low + local_index;
194 Node<SPACE_DIM>* p_node = this->mpMesh->GetNode(global_index);
195 mGgapDistributed[local_index] = mGGap;//assign default uniform value everywhere first
196
197 // Then change where and if necessary
198 for (unsigned het_index = 0; het_index < mGgapHeterogeneityRegions.size(); het_index++)
199 {
200 if (mGgapHeterogeneityRegions[het_index]->DoesContain(p_node->GetPoint()))
202 mGgapDistributed[local_index] = mGgapValues[het_index];
203 }
204 }
205 }
206 }
207 // LCOV_EXCL_START
208 catch (const Exception& e)
209 {
211 throw e;
212 }
213 // LCOV_EXCL_STOP
214
217
218template <unsigned SPACE_DIM>
220{
221 HeartEventHandler::BeginEvent(HeartEventHandler::READ_MESH);
222 this->mpConfig = HeartConfig::Instance();
223
224 if (this->mpConfig->IsMeshProvided() && this->mpConfig->GetLoadMesh())
225 {
226 assert(this->mFibreFilePathNoExtension != "");
228 switch (this->mpConfig->GetConductivityMedia())
229 {
230 case cp::media_type::Orthotropic:
231 {
232 mpIntracellularConductivityTensorsSecondCell = new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
233 FileFinder ortho_file(this->mFibreFilePathNoExtension + ".ortho", RelativeTo::AbsoluteOrCwd);
234 assert(ortho_file.Exists());
235 mpIntracellularConductivityTensorsSecondCell->SetFibreOrientationFile(ortho_file);
236 break;
238
239 case cp::media_type::Axisymmetric:
240 {
241 mpIntracellularConductivityTensorsSecondCell = new AxisymmetricConductivityTensors<SPACE_DIM,SPACE_DIM>;
242 FileFinder axi_file(this->mFibreFilePathNoExtension + ".axi", RelativeTo::AbsoluteOrCwd);
243 assert(axi_file.Exists());
244 mpIntracellularConductivityTensorsSecondCell->SetFibreOrientationFile(axi_file);
245 break;
246 }
247
248 case cp::media_type::NoFibreOrientation:
249 mpIntracellularConductivityTensorsSecondCell = new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
250 break;
252 default:
254 }
255 }
256 else // Slab defined in config file or SetMesh() called; no fibre orientation assumed
257 {
258 mpIntracellularConductivityTensorsSecondCell = new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
260
261 // this definition must be here (and not inside the if statement) because SetNonConstantConductivities() will keep
262 // a pointer to it and we don't want it to go out of scope before Init() is called
263 unsigned num_elements = this->mpMesh->GetNumElements();
264 std::vector<c_vector<double, SPACE_DIM> > hetero_intra_conductivities;
265
266 c_vector<double, SPACE_DIM> intra_conductivities;
267 this->mpConfig->GetIntracellularConductivities(intra_conductivities);//this one is used just for resizing
268
269 if (this->mpConfig->GetConductivityHeterogeneitiesProvided())
270 {
271 try
272 {
273 assert(hetero_intra_conductivities.size()==0);
274 hetero_intra_conductivities.resize(num_elements);
275 for(auto& elem : hetero_intra_conductivities)
276 {
277 ExplicitCVectorCopy<SPACE_DIM>(intra_conductivities, elem);
278 }
279 }
280 // LCOV_EXCL_START
281 catch(std::bad_alloc &badAlloc)
282 {
283
284 std::cout << "Failed to allocate std::vector of size " << num_elements << std::endl;
286 throw badAlloc;
287 }
288 // LCOV_EXCL_STOP
289
291
292 std::vector<boost::shared_ptr<AbstractChasteRegion<SPACE_DIM> > > conductivities_heterogeneity_areas;
293 std::vector< c_vector<double,3> > intra_h_conductivities;
294 std::vector< c_vector<double,3> > extra_h_conductivities;
295 HeartConfig::Instance()->GetConductivityHeterogeneities(conductivities_heterogeneity_areas,
296 intra_h_conductivities,
297 extra_h_conductivities);
298 unsigned local_element_index = 0;
299 for (typename AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>::ElementIterator it = this->mpMesh->GetElementIteratorBegin();
300 it != this->mpMesh->GetElementIteratorEnd();
301 ++it)
302 {
303 //unsigned element_index = it->GetIndex();
304 // if element centroid is contained in the region
305 ChastePoint<SPACE_DIM> element_centroid(it->CalculateCentroid());
306 for (unsigned region_index=0; region_index< conductivities_heterogeneity_areas.size(); region_index++)
307 {
308 if (conductivities_heterogeneity_areas[region_index]->DoesContain(element_centroid))
309 {
310 // We don't use ublas vector assignment here, because we might be getting a subvector of a 3-vector
311 for (unsigned i=0; i<SPACE_DIM; i++)
313 hetero_intra_conductivities[local_element_index][i] = intra_h_conductivities[region_index][i];
314 }
315 }
316 }
317 local_element_index++;
319
320 mpIntracellularConductivityTensorsSecondCell->SetNonConstantConductivities(&hetero_intra_conductivities);
321 }
322 else
323 {
324 mpIntracellularConductivityTensorsSecondCell->SetConstantConductivities(mIntracellularConductivitiesSecondCell);
325 }
326
327 mpIntracellularConductivityTensorsSecondCell->Init(this->mpMesh);
328 HeartEventHandler::EndEvent(HeartEventHandler::READ_MESH);
329}
330
331template <unsigned SPACE_DIM>
333{
334 return mUserSuppliedExtracellularStimulus;
335}
336
337template <unsigned SPACE_DIM>
339{
340 mUserSuppliedExtracellularStimulus = flag;
341}
343template <unsigned SPACE_DIM>
344const std::vector<AbstractCardiacCellInterface*>& ExtendedBidomainTissue<SPACE_DIM>::rGetSecondCellsDistributed() const
345{
346 return mCellsDistributedSecondCell;
348
349template <unsigned SPACE_DIM>
351{
352 return mGgapDistributed;
353}
354
355template <unsigned SPACE_DIM>
356const std::vector<boost::shared_ptr<AbstractStimulusFunction> >& ExtendedBidomainTissue<SPACE_DIM>::rGetExtracellularStimulusDistributed() const
358 return mExtracellularStimuliDistributed;
359}
360
361template <unsigned SPACE_DIM>
363{
364 if (this->mpConfig->IsMeshProvided() && this->mpConfig->GetLoadMesh())
365 {
366 assert(this->mFibreFilePathNoExtension != "");
367 switch (this->mpConfig->GetConductivityMedia())
368 {
369 case cp::media_type::Orthotropic:
370 {
371 mpExtracellularConductivityTensors = new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
372 FileFinder ortho_file(this->mFibreFilePathNoExtension + ".ortho", RelativeTo::AbsoluteOrCwd);
373 assert(ortho_file.Exists());
374 mpExtracellularConductivityTensors->SetFibreOrientationFile(ortho_file);
375 break;
376 }
378 case cp::media_type::Axisymmetric:
379 {
380 mpExtracellularConductivityTensors = new AxisymmetricConductivityTensors<SPACE_DIM,SPACE_DIM>;
381 FileFinder axi_file(this->mFibreFilePathNoExtension + ".axi", RelativeTo::AbsoluteOrCwd);
382 assert(axi_file.Exists());
383 mpExtracellularConductivityTensors->SetFibreOrientationFile(axi_file);
384 break;
385 }
386
387 case cp::media_type::NoFibreOrientation:
388 mpExtracellularConductivityTensors = new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
389 break;
390
391 default:
393 }
394 }
395 else // no fibre orientation assumed
396 {
397 mpExtracellularConductivityTensors = new OrthotropicConductivityTensors<SPACE_DIM,SPACE_DIM>;
398 }
399
400 c_vector<double, SPACE_DIM> extra_conductivities;
401 this->mpConfig->GetExtracellularConductivities(extra_conductivities);
402
403 // this definition must be here (and not inside the if statement) because SetNonConstantConductivities() will keep
404 // a pointer to it and we don't want it to go out of scope before Init() is called
405 unsigned num_elements = this->mpMesh->GetNumElements();
406 std::vector<c_vector<double, SPACE_DIM> > hetero_extra_conductivities;
407
408 if (this->mpConfig->GetConductivityHeterogeneitiesProvided())
409 {
410 try
411 {
412 assert(hetero_extra_conductivities.size()==0);
413 //initialise with the values of teh default conductivity tensor
414 hetero_extra_conductivities.resize(num_elements);
415 for(auto& elem : hetero_extra_conductivities)
416 {
417 ExplicitCVectorCopy<SPACE_DIM>(extra_conductivities, elem);
418 }
419 }
420 // LCOV_EXCL_START
421 catch(std::bad_alloc &badAlloc)
422 {
423 std::cout << "Failed to allocate std::vector of size " << num_elements << std::endl;
425 throw badAlloc;
426 }
427 // LCOV_EXCL_STOP
428
430
431 std::vector<boost::shared_ptr<AbstractChasteRegion<SPACE_DIM> > > conductivities_heterogeneity_areas;
432 std::vector< c_vector<double,3> > intra_h_conductivities;
433 std::vector< c_vector<double,3> > extra_h_conductivities;
434 HeartConfig::Instance()->GetConductivityHeterogeneities(conductivities_heterogeneity_areas,
435 intra_h_conductivities,
436 extra_h_conductivities);
437 unsigned local_element_index = 0;
438 for (typename AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>::ElementIterator iter = (this->mpMesh)->GetElementIteratorBegin();
439 iter != (this->mpMesh)->GetElementIteratorEnd();
440 ++iter)
441 {
442 //unsigned element_index = iter->GetIndex();
443 // if element centroid is contained in the region
444 ChastePoint<SPACE_DIM> element_centroid(iter->CalculateCentroid());
445 for (unsigned region_index=0; region_index< conductivities_heterogeneity_areas.size(); region_index++)
446 {
447 // If element centroid is contained in the region
448 if (conductivities_heterogeneity_areas[region_index]->DoesContain(element_centroid))
449 {
450 // We don't use ublas vector assignment here, because we might be getting a subvector of a 3-vector
451 for (unsigned i=0; i<SPACE_DIM; i++)
452 {
453 hetero_extra_conductivities[local_element_index][i] = extra_h_conductivities[region_index][i];
454 }
455 }
456 }
457 local_element_index++;
458 }
459
460 mpExtracellularConductivityTensors->SetNonConstantConductivities(&hetero_extra_conductivities);
461 }
462 else
463 {
464 mpExtracellularConductivityTensors->SetConstantConductivities(extra_conductivities);
465 }
466 mpExtracellularConductivityTensors->Init(this->mpMesh);
467}
468
469template <unsigned SPACE_DIM>
471{
472 // Delete (second) cells
473 for (std::vector<AbstractCardiacCellInterface*>::iterator cell_iterator = mCellsDistributedSecondCell.begin();
474 cell_iterator != mCellsDistributedSecondCell.end();
475 ++cell_iterator)
476 {
477 delete (*cell_iterator);
478 }
479
480 if (mpExtracellularConductivityTensors)
481 {
482 delete mpExtracellularConductivityTensors;
483 }
484
485 if (mpIntracellularConductivityTensorsSecondCell)
486 {
487 delete mpIntracellularConductivityTensorsSecondCell;
488 }
489}
490
491template <unsigned SPACE_DIM>
493{
494 for (unsigned i = 0; i < SPACE_DIM; i++)
495 {
496 mIntracellularConductivitiesSecondCell[i] = conductivities[i];
497 }
498}
499
500template <unsigned SPACE_DIM>
502{
503 return mIntracellularConductivitiesSecondCell;
504}
505
506template <unsigned SPACE_DIM>
507const c_matrix<double, SPACE_DIM, SPACE_DIM>& ExtendedBidomainTissue<SPACE_DIM>::rGetExtracellularConductivityTensor(unsigned elementIndex)
508{
509 assert(mpExtracellularConductivityTensors);
510 if (this->mpConductivityModifier==NULL)
511 {
512 return (*mpExtracellularConductivityTensors)[elementIndex];
513 }
514 else
515 {
516 return this->mpConductivityModifier->rGetModifiedConductivityTensor(elementIndex, (*mpExtracellularConductivityTensors)[elementIndex], 1u);
517 }
518}
519
520template <unsigned SPACE_DIM>
521const c_matrix<double, SPACE_DIM, SPACE_DIM>& ExtendedBidomainTissue<SPACE_DIM>::rGetIntracellularConductivityTensorSecondCell(unsigned elementIndex)
522{
523 assert(mpIntracellularConductivityTensorsSecondCell);
524 if (this->mpConductivityModifier==NULL)
525 {
526 return (*mpIntracellularConductivityTensorsSecondCell)[elementIndex];
527 }
528 else
529 {
530 return this->mpConductivityModifier->rGetModifiedConductivityTensor(elementIndex, (*mpIntracellularConductivityTensorsSecondCell)[elementIndex], 2u);
531 }
532}
533
534template <unsigned SPACE_DIM>
536{
537 return mCellsDistributedSecondCell[globalIndex - this->mpDistributedVectorFactory->GetLow()];
538}
539
540template <unsigned SPACE_DIM>
541boost::shared_ptr<AbstractStimulusFunction> ExtendedBidomainTissue<SPACE_DIM>::GetExtracellularStimulus( unsigned globalIndex )
542{
543 return mExtracellularStimuliDistributed[globalIndex - this->mpDistributedVectorFactory->GetLow()];
544}
545
546template <unsigned SPACE_DIM>
547void ExtendedBidomainTissue<SPACE_DIM>::SolveCellSystems(Vec existingSolution, double time, double nextTime, bool updateVoltage)
548{
549 HeartEventHandler::BeginEvent(HeartEventHandler::SOLVE_ODES);
550
551 DistributedVector dist_solution = this->mpDistributedVectorFactory->CreateDistributedVector(existingSolution);
552 DistributedVector::Stripe V_first_cell(dist_solution, 0);
553 DistributedVector::Stripe V_second_cell(dist_solution, 1);
554 DistributedVector::Stripe phi_e(dist_solution, 2);
555
556 for (DistributedVector::Iterator index = dist_solution.Begin();
557 index != dist_solution.End();
558 ++index)
559 {
560 // overwrite the voltage with the input value
561 this->mCellsDistributed[index.Local]->SetVoltage( V_first_cell[index] );
562 mCellsDistributedSecondCell[index.Local]->SetVoltage( V_second_cell[index] );
563 try
564 {
565 // solve
566 // Note: Voltage should not be updated. GetIIonic will be called later
567 // and needs the old voltage. The voltage will be updated from the pde.
568 this->mCellsDistributed[index.Local]->ComputeExceptVoltage(time, nextTime);
569 mCellsDistributedSecondCell[index.Local]->ComputeExceptVoltage(time, nextTime);
570 }
571 // LCOV_EXCL_START
572 catch (Exception &e)
573 {
575 throw e;
576 }
577 // LCOV_EXCL_STOP
578
579 // update the Iionic and stimulus caches
580 this->UpdateCaches(index.Global, index.Local, nextTime);//in parent class
581 UpdateAdditionalCaches(index.Global, index.Local, nextTime);//extended bidomain specific caches
582 }
584 HeartEventHandler::EndEvent(HeartEventHandler::SOLVE_ODES);
585
586 HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
587 if (this->mDoCacheReplication)
588 {
589 this->ReplicateCaches();
590 ReplicateAdditionalCaches();//extended bidomain specific caches
591 }
592 HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
593}
594
595template <unsigned SPACE_DIM>
596void ExtendedBidomainTissue<SPACE_DIM>::UpdateAdditionalCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
597{
598 mIionicCacheReplicatedSecondCell[globalIndex] = mCellsDistributedSecondCell[localIndex]->GetIIonic();
599 mIntracellularStimulusCacheReplicatedSecondCell[globalIndex] = mCellsDistributedSecondCell[localIndex]->GetIntracellularStimulus(nextTime);
600 mExtracellularStimulusCacheReplicated[globalIndex] = mExtracellularStimuliDistributed[localIndex]->GetStimulus(nextTime);
601 mGgapCacheReplicated[globalIndex] = mGgapDistributed[localIndex];
602}
603
604template <unsigned SPACE_DIM>
606{
607 mIionicCacheReplicatedSecondCell.Replicate(this->mpDistributedVectorFactory->GetLow(), this->mpDistributedVectorFactory->GetHigh());
608 mIntracellularStimulusCacheReplicatedSecondCell.Replicate(this->mpDistributedVectorFactory->GetLow(), this->mpDistributedVectorFactory->GetHigh());
609 mExtracellularStimulusCacheReplicated.Replicate(this->mpDistributedVectorFactory->GetLow(), this->mpDistributedVectorFactory->GetHigh());
610 mGgapCacheReplicated.Replicate(this->mpDistributedVectorFactory->GetLow(), this->mpDistributedVectorFactory->GetHigh());
611}
612
613template <unsigned SPACE_DIM>
618
619template <unsigned SPACE_DIM>
621{
622 return mIntracellularStimulusCacheReplicatedSecondCell;
623}
624
625template <unsigned SPACE_DIM>
627{
628 return mExtracellularStimulusCacheReplicated;
629}
630
631template <unsigned SPACE_DIM>
636
637template <unsigned SPACE_DIM>
639{
640 return mAmFirstCell;
641}
642
643template <unsigned SPACE_DIM>
645{
646 return mAmSecondCell;
647}
648
649template <unsigned SPACE_DIM>
651{
652 return mAmGap;
653}
654
655template <unsigned SPACE_DIM>
657{
658 return mCmFirstCell;
659}
660
661template <unsigned SPACE_DIM>
663{
664 return mCmSecondCell;
665}
666
667template <unsigned SPACE_DIM>
669{
670 return mGGap;
671}
672
673template <unsigned SPACE_DIM>
675{
676 mAmFirstCell = value;
677}
678
679template <unsigned SPACE_DIM>
681{
682 mAmSecondCell = value;
683}
684
685template <unsigned SPACE_DIM>
687{
688 mAmGap = value;
689}
690
691template <unsigned SPACE_DIM>
693{
694 mGGap = value;
695}
696
697template <unsigned SPACE_DIM>
699{
700 mCmFirstCell = value;
701}
702
703template <unsigned SPACE_DIM>
705{
706 mCmSecondCell = value;
707}
708
709// Explicit instantiation
710template class ExtendedBidomainTissue<1>;
711template class ExtendedBidomainTissue<2>;
712template class ExtendedBidomainTissue<3>;
713
714// Serialization for Boost >= 1.36
#define NEVER_REACHED
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
virtual AbstractCardiacCellInterface * CreateCardiacCellForNode(Node< SPACE_DIM > *pNode)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * GetMesh()
virtual void FinaliseCellCreation(std::vector< AbstractCardiacCellInterface * > *pCellsDistributed, unsigned lo, unsigned hi)
DistributedVectorFactory * mpDistributedVectorFactory
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
virtual boost::shared_ptr< AbstractStimulusFunction > CreateStimulusForNode(Node< SPACE_DIM > *pNode)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * GetMesh()
ReplicatableVector mExtracellularStimulusCacheReplicated
void UpdateAdditionalCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
ReplicatableVector & rGetExtracellularStimulusCacheReplicated()
const std::vector< boost::shared_ptr< AbstractStimulusFunction > > & rGetExtracellularStimulusDistributed() const
virtual void SolveCellSystems(Vec existingSolution, double time, double nextTime, bool updateVoltage=false)
void SetGgapHeterogeneities(std::vector< boost::shared_ptr< AbstractChasteRegion< SPACE_DIM > > > &rGgapHeterogeneityRegions, std::vector< double > rGgapValues)
ExtendedBidomainTissue(AbstractCardiacCellFactory< SPACE_DIM > *pCellFactory, AbstractCardiacCellFactory< SPACE_DIM > *pCellFactorySecondCell, AbstractStimulusFactory< SPACE_DIM > *pExtracellularStimulusFactory)
ReplicatableVector mIionicCacheReplicatedSecondCell
ReplicatableVector & rGetIionicCacheReplicatedSecondCell()
ReplicatableVector & rGetGgapCacheReplicated()
const std::vector< AbstractCardiacCellInterface * > & rGetSecondCellsDistributed() const
std::vector< AbstractCardiacCellInterface * > mCellsDistributedSecondCell
ReplicatableVector mIntracellularStimulusCacheReplicatedSecondCell
boost::shared_ptr< AbstractStimulusFunction > GetExtracellularStimulus(unsigned globalIndex)
std::vector< boost::shared_ptr< AbstractStimulusFunction > > mExtracellularStimuliDistributed
ReplicatableVector & rGetIntracellularStimulusCacheReplicatedSecondCell()
c_vector< double, SPACE_DIM > GetIntracellularConductivitiesSecondCell() const
const c_matrix< double, SPACE_DIM, SPACE_DIM > & rGetExtracellularConductivityTensor(unsigned elementIndex)
std::vector< double > mGgapDistributed
AbstractCardiacCellInterface * GetCardiacSecondCell(unsigned globalIndex)
const std::vector< double > & rGetGapsDistributed() const
void SetUserSuppliedExtracellularStimulus(bool flag)
void SetIntracellularConductivitiesSecondCell(c_vector< double, SPACE_DIM > conductivities)
const c_matrix< double, SPACE_DIM, SPACE_DIM > & rGetIntracellularConductivityTensorSecondCell(unsigned elementIndex)
ReplicatableVector mGgapCacheReplicated
bool Exists() const
void GetConductivityHeterogeneities(std::vector< boost::shared_ptr< AbstractChasteRegion< DIM > > > &conductivitiesHeterogeneityAreas, std::vector< c_vector< double, 3 > > &intraConductivities, std::vector< c_vector< double, 3 > > &extraConductivities) const
static HeartConfig * Instance()
Definition Node.hpp:59
ChastePoint< SPACE_DIM > GetPoint() const
Definition Node.cpp:133
static void ReplicateException(bool flag)
void Resize(unsigned size)