36#include "ExtendedBidomainTissue.hpp"
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"
55template <
unsigned SPACE_DIM>
56void ExplicitCVectorCopy(c_vector<double, SPACE_DIM>& rFrom, c_vector<double, SPACE_DIM>& rTo)
58 if constexpr (SPACE_DIM == 1u)
71template <
unsigned SPACE_DIM>
76 mpIntracellularConductivityTensorsSecondCell(NULL),
77 mUserSuppliedExtracellularStimulus(false)
81 assert(pCellFactorySecondCell != NULL);
82 assert(pCellFactorySecondCell->
GetMesh() != NULL);
84 assert(pExtracellularStimulusFactory != NULL);
85 assert(pExtracellularStimulusFactory->
GetMesh() != NULL);
96 for (
unsigned local_index = 0; local_index < num_local_nodes; local_index++)
98 unsigned global_index = local_index + ownership_range_low;
108 this->mpDistributedVectorFactory->GetHigh());
122 delete (*cell_iterator);
141template <
unsigned SPACE_DIM>
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)
149 mpIntracellularConductivityTensorsSecondCell(NULL),
150 mIntracellularConductivitiesSecondCell(intracellularConductivitiesSecondCell),
151 mCellsDistributedSecondCell(rSecondCellsDistributed),
152 mExtracellularStimuliDistributed(rExtraStimuliDistributed),
153 mGgapDistributed(rGgapsDistributed),
154 mUserSuppliedExtracellularStimulus(false)
171template <
unsigned SPACE_DIM>
173 std::vector<double> rGgapValues)
175 assert( rGgapHeterogeneityRegions.size() == rGgapValues.size() );
176 mGgapHeterogeneityRegions = rGgapHeterogeneityRegions;
177 mGgapValues =rGgapValues;
180template <
unsigned SPACE_DIM>
183 assert(mGgapHeterogeneityRegions.size() == mGgapValues.size());
184 assert(this->mpMesh != NULL);
186 unsigned ownership_range_low = this->mpDistributedVectorFactory->GetLow();
187 unsigned num_local_nodes = this->mpDistributedVectorFactory->GetLocalOwnership();
188 assert(mGgapDistributed.size() == num_local_nodes);
191 for (
unsigned local_index = 0; local_index < num_local_nodes; local_index++)
193 unsigned global_index = ownership_range_low + local_index;
195 mGgapDistributed[local_index] = mGGap;
198 for (
unsigned het_index = 0; het_index < mGgapHeterogeneityRegions.size(); het_index++)
200 if (mGgapHeterogeneityRegions[het_index]->DoesContain(p_node->
GetPoint()))
202 mGgapDistributed[local_index] = mGgapValues[het_index];
218template <
unsigned SPACE_DIM>
224 if (this->mpConfig->IsMeshProvided() && this->mpConfig->GetLoadMesh())
226 assert(this->mFibreFilePathNoExtension !=
"");
228 switch (this->mpConfig->GetConductivityMedia())
230 case cp::media_type::Orthotropic:
234 assert(ortho_file.
Exists());
235 mpIntracellularConductivityTensorsSecondCell->SetFibreOrientationFile(ortho_file);
239 case cp::media_type::Axisymmetric:
243 assert(axi_file.
Exists());
244 mpIntracellularConductivityTensorsSecondCell->SetFibreOrientationFile(axi_file);
248 case cp::media_type::NoFibreOrientation:
263 unsigned num_elements = this->mpMesh->GetNumElements();
264 std::vector<c_vector<double, SPACE_DIM> > hetero_intra_conductivities;
266 c_vector<double, SPACE_DIM> intra_conductivities;
267 this->mpConfig->GetIntracellularConductivities(intra_conductivities);
269 if (this->mpConfig->GetConductivityHeterogeneitiesProvided())
273 assert(hetero_intra_conductivities.size()==0);
274 hetero_intra_conductivities.resize(num_elements);
275 for(
auto& elem : hetero_intra_conductivities)
277 ExplicitCVectorCopy<SPACE_DIM>(intra_conductivities, elem);
281 catch(std::bad_alloc &badAlloc)
284 std::cout <<
"Failed to allocate std::vector of size " << num_elements << std::endl;
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;
296 intra_h_conductivities,
297 extra_h_conductivities);
298 unsigned local_element_index = 0;
300 it != this->mpMesh->GetElementIteratorEnd();
306 for (
unsigned region_index=0; region_index< conductivities_heterogeneity_areas.size(); region_index++)
308 if (conductivities_heterogeneity_areas[region_index]->DoesContain(element_centroid))
311 for (
unsigned i=0; i<SPACE_DIM; i++)
313 hetero_intra_conductivities[local_element_index][i] = intra_h_conductivities[region_index][i];
317 local_element_index++;
320 mpIntracellularConductivityTensorsSecondCell->SetNonConstantConductivities(&hetero_intra_conductivities);
324 mpIntracellularConductivityTensorsSecondCell->SetConstantConductivities(mIntracellularConductivitiesSecondCell);
327 mpIntracellularConductivityTensorsSecondCell->Init(this->mpMesh);
331template <
unsigned SPACE_DIM>
334 return mUserSuppliedExtracellularStimulus;
337template <
unsigned SPACE_DIM>
340 mUserSuppliedExtracellularStimulus = flag;
343template <
unsigned SPACE_DIM>
346 return mCellsDistributedSecondCell;
349template <
unsigned SPACE_DIM>
352 return mGgapDistributed;
355template <
unsigned SPACE_DIM>
358 return mExtracellularStimuliDistributed;
361template <
unsigned SPACE_DIM>
364 if (this->mpConfig->IsMeshProvided() && this->mpConfig->GetLoadMesh())
366 assert(this->mFibreFilePathNoExtension !=
"");
367 switch (this->mpConfig->GetConductivityMedia())
369 case cp::media_type::Orthotropic:
373 assert(ortho_file.Exists());
374 mpExtracellularConductivityTensors->SetFibreOrientationFile(ortho_file);
378 case cp::media_type::Axisymmetric:
383 mpExtracellularConductivityTensors->SetFibreOrientationFile(axi_file);
387 case cp::media_type::NoFibreOrientation:
400 c_vector<double, SPACE_DIM> extra_conductivities;
401 this->mpConfig->GetExtracellularConductivities(extra_conductivities);
405 unsigned num_elements = this->mpMesh->GetNumElements();
406 std::vector<c_vector<double, SPACE_DIM> > hetero_extra_conductivities;
408 if (this->mpConfig->GetConductivityHeterogeneitiesProvided())
412 assert(hetero_extra_conductivities.size()==0);
414 hetero_extra_conductivities.resize(num_elements);
415 for(
auto& elem : hetero_extra_conductivities)
417 ExplicitCVectorCopy<SPACE_DIM>(extra_conductivities, elem);
421 catch(std::bad_alloc &badAlloc)
423 std::cout <<
"Failed to allocate std::vector of size " << num_elements << std::endl;
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;
435 intra_h_conductivities,
436 extra_h_conductivities);
437 unsigned local_element_index = 0;
439 iter != (this->mpMesh)->GetElementIteratorEnd();
445 for (
unsigned region_index=0; region_index< conductivities_heterogeneity_areas.size(); region_index++)
448 if (conductivities_heterogeneity_areas[region_index]->DoesContain(element_centroid))
451 for (
unsigned i=0; i<SPACE_DIM; i++)
453 hetero_extra_conductivities[local_element_index][i] = extra_h_conductivities[region_index][i];
457 local_element_index++;
460 mpExtracellularConductivityTensors->SetNonConstantConductivities(&hetero_extra_conductivities);
464 mpExtracellularConductivityTensors->SetConstantConductivities(extra_conductivities);
466 mpExtracellularConductivityTensors->Init(this->mpMesh);
469template <
unsigned SPACE_DIM>
473 for (std::vector<AbstractCardiacCellInterface*>::iterator cell_iterator = mCellsDistributedSecondCell.begin();
474 cell_iterator != mCellsDistributedSecondCell.end();
477 delete (*cell_iterator);
480 if (mpExtracellularConductivityTensors)
482 delete mpExtracellularConductivityTensors;
485 if (mpIntracellularConductivityTensorsSecondCell)
487 delete mpIntracellularConductivityTensorsSecondCell;
491template <
unsigned SPACE_DIM>
494 for (
unsigned i = 0; i < SPACE_DIM; i++)
496 mIntracellularConductivitiesSecondCell[i] = conductivities[i];
500template <
unsigned SPACE_DIM>
503 return mIntracellularConductivitiesSecondCell;
506template <
unsigned SPACE_DIM>
509 assert(mpExtracellularConductivityTensors);
510 if (this->mpConductivityModifier==NULL)
512 return (*mpExtracellularConductivityTensors)[elementIndex];
516 return this->mpConductivityModifier->rGetModifiedConductivityTensor(elementIndex, (*mpExtracellularConductivityTensors)[elementIndex], 1u);
520template <
unsigned SPACE_DIM>
523 assert(mpIntracellularConductivityTensorsSecondCell);
524 if (this->mpConductivityModifier==NULL)
526 return (*mpIntracellularConductivityTensorsSecondCell)[elementIndex];
530 return this->mpConductivityModifier->rGetModifiedConductivityTensor(elementIndex, (*mpIntracellularConductivityTensorsSecondCell)[elementIndex], 2u);
534template <
unsigned SPACE_DIM>
537 return mCellsDistributedSecondCell[globalIndex - this->mpDistributedVectorFactory->GetLow()];
540template <
unsigned SPACE_DIM>
543 return mExtracellularStimuliDistributed[globalIndex - this->mpDistributedVectorFactory->GetLow()];
546template <
unsigned SPACE_DIM>
551 DistributedVector dist_solution = this->mpDistributedVectorFactory->CreateDistributedVector(existingSolution);
557 index != dist_solution.
End();
561 this->mCellsDistributed[index.Local]->SetVoltage( V_first_cell[index] );
562 mCellsDistributedSecondCell[index.Local]->SetVoltage( V_second_cell[index] );
568 this->mCellsDistributed[index.Local]->ComputeExceptVoltage(time, nextTime);
569 mCellsDistributedSecondCell[index.Local]->ComputeExceptVoltage(time, nextTime);
580 this->UpdateCaches(index.Global, index.Local, nextTime);
581 UpdateAdditionalCaches(index.Global, index.Local, nextTime);
587 if (this->mDoCacheReplication)
589 this->ReplicateCaches();
590 ReplicateAdditionalCaches();
595template <
unsigned SPACE_DIM>
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];
604template <
unsigned SPACE_DIM>
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());
613template <
unsigned SPACE_DIM>
616 return mIionicCacheReplicatedSecondCell;
619template <
unsigned SPACE_DIM>
622 return mIntracellularStimulusCacheReplicatedSecondCell;
625template <
unsigned SPACE_DIM>
628 return mExtracellularStimulusCacheReplicated;
631template <
unsigned SPACE_DIM>
634 return mGgapCacheReplicated;
637template <
unsigned SPACE_DIM>
643template <
unsigned SPACE_DIM>
646 return mAmSecondCell;
649template <
unsigned SPACE_DIM>
655template <
unsigned SPACE_DIM>
661template <
unsigned SPACE_DIM>
664 return mCmSecondCell;
667template <
unsigned SPACE_DIM>
673template <
unsigned SPACE_DIM>
676 mAmFirstCell = value;
679template <
unsigned SPACE_DIM>
682 mAmSecondCell = value;
685template <
unsigned SPACE_DIM>
691template <
unsigned SPACE_DIM>
697template <
unsigned SPACE_DIM>
700 mCmFirstCell = value;
703template <
unsigned SPACE_DIM>
706 mCmSecondCell = value;
#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)
virtual unsigned GetNumberOfCells()
void SetUsedInTissueSimulation(bool tissue=true)
DistributedVectorFactory * mpDistributedVectorFactory
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
virtual unsigned GetNumberOfCells()
virtual boost::shared_ptr< AbstractStimulusFunction > CreateStimulusForNode(Node< SPACE_DIM > *pNode)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * GetMesh()
unsigned GetProblemSize() const
unsigned GetLocalOwnership() const
ReplicatableVector mExtracellularStimulusCacheReplicated
void UpdateAdditionalCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
void ReplicateAdditionalCaches()
void SetGGap(double value)
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)
void CreateExtracellularConductivityTensors()
ExtendedBidomainTissue(AbstractCardiacCellFactory< SPACE_DIM > *pCellFactory, AbstractCardiacCellFactory< SPACE_DIM > *pCellFactorySecondCell, AbstractStimulusFactory< SPACE_DIM > *pExtracellularStimulusFactory)
ReplicatableVector mIionicCacheReplicatedSecondCell
ReplicatableVector & rGetIionicCacheReplicatedSecondCell()
void SetCmFirstCell(double value)
ReplicatableVector & rGetGgapCacheReplicated()
const std::vector< AbstractCardiacCellInterface * > & rGetSecondCellsDistributed() const
std::vector< AbstractCardiacCellInterface * > mCellsDistributedSecondCell
ReplicatableVector mIntracellularStimulusCacheReplicatedSecondCell
void SetAmFirstCell(double value)
boost::shared_ptr< AbstractStimulusFunction > GetExtracellularStimulus(unsigned globalIndex)
bool HasTheUserSuppliedExtracellularStimulus()
std::vector< boost::shared_ptr< AbstractStimulusFunction > > mExtracellularStimuliDistributed
void SetAmSecondCell(double value)
ReplicatableVector & rGetIntracellularStimulusCacheReplicatedSecondCell()
c_vector< double, SPACE_DIM > GetIntracellularConductivitiesSecondCell() const
void CreateIntracellularConductivityTensorSecondCell()
const c_matrix< double, SPACE_DIM, SPACE_DIM > & rGetExtracellularConductivityTensor(unsigned elementIndex)
std::vector< double > mGgapDistributed
virtual ~ExtendedBidomainTissue()
AbstractCardiacCellInterface * GetCardiacSecondCell(unsigned globalIndex)
const std::vector< double > & rGetGapsDistributed() const
void SetUserSuppliedExtracellularStimulus(bool flag)
void SetAmGap(double value)
void SetIntracellularConductivitiesSecondCell(c_vector< double, SPACE_DIM > conductivities)
const c_matrix< double, SPACE_DIM, SPACE_DIM > & rGetIntracellularConductivityTensorSecondCell(unsigned elementIndex)
void CreateGGapConductivities()
ReplicatableVector mGgapCacheReplicated
void SetCmSecondCell(double value)
static void BeginEvent(unsigned event)
static void EndEvent(unsigned event)
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()
ChastePoint< SPACE_DIM > GetPoint() const
void Resize(unsigned size)