Chaste  Release::3.4
AbstractCardiacTissue.hpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 #ifndef ABSTRACTCARDIACTISSUE_HPP_
36 #define ABSTRACTCARDIACTISSUE_HPP_
37 
38 #include <set>
39 #include <vector>
40 #include <boost/shared_ptr.hpp>
41 
42 #include "UblasMatrixInclude.hpp"
43 
44 #include "ChasteSerialization.hpp"
45 #include "ClassIsAbstract.hpp"
47 #include <boost/serialization/base_object.hpp>
48 #include <boost/serialization/shared_ptr.hpp>
49 #include <boost/serialization/vector.hpp>
50 #include <boost/serialization/string.hpp>
51 #include <boost/serialization/split_member.hpp>
52 
53 #include "AbstractCardiacCellInterface.hpp"
54 #include "FakeBathCell.hpp"
55 #include "AbstractCardiacCellFactory.hpp"
56 #include "AbstractConductivityTensors.hpp"
57 #include "AbstractPurkinjeCellFactory.hpp"
58 #include "ReplicatableVector.hpp"
59 #include "HeartConfig.hpp"
60 #include "ArchiveLocationInfo.hpp"
61 #include "AbstractDynamicallyLoadableEntity.hpp"
62 #include "DynamicModelLoaderRegistry.hpp"
63 #include "AbstractConductivityModifier.hpp"
64 
77 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM = ELEMENT_DIM>
78 class AbstractCardiacTissue : private boost::noncopyable
79 {
80 private:
81 
83  friend class boost::serialization::access;
84  friend class TestMonodomainTissue;
85 
92  template<class Archive>
93  void save(Archive & archive, const unsigned int version) const
94  {
95  if (version >= 3)
96  {
97  archive & mHasPurkinje;
98  }
99  if (version >= 2)
100  {
101  archive & mExchangeHalos;
102  }
103  // Don't use the std::vector serialization for cardiac cells, so that we can load them
104  // more cleverly when migrating checkpoints.
106 
107  // archive & mpMesh; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
108  // archive & mpIntracellularConductivityTensors; Loaded from HeartConfig every time constructor is called
109  if (HeartConfig::Instance()->IsMeshProvided() && HeartConfig::Instance()->GetLoadMesh())
110  {
111  switch (HeartConfig::Instance()->GetConductivityMedia())
112  {
113  case cp::media_type::Orthotropic:
114  {
116  assert(source_file.Exists());
118 
119  TRY_IF_MASTER(source_file.CopyTo(dest_file));
120  break;
121  }
122 
123  case cp::media_type::Axisymmetric:
124  {
126  assert(source_file.Exists());
129 
130  TRY_IF_MASTER(source_file.CopyTo(dest_file));
131  break;
132  }
133 
134  case cp::media_type::NoFibreOrientation:
135  break;
136 
137  default :
139  }
140  }
141 
142  // archive & mIionicCacheReplicated; // will be regenerated
143  // archive & mIntracellularStimulusCacheReplicated; // will be regenerated
144  archive & mDoCacheReplication;
145  // archive & mMeshUnarchived; Not archived since set to true when archiving constructor is called.
146 
148 
149  // Paranoia: check we agree with the mesh on who owns what
150  assert(mpDistributedVectorFactory == mpMesh->GetDistributedVectorFactory());
151  assert(mpDistributedVectorFactory->GetLow()==mpMesh->GetDistributedVectorFactory()->GetLow());
152  assert(mpDistributedVectorFactory->GetLocalOwnership()==mpMesh->GetDistributedVectorFactory()->GetLocalOwnership());
153  }
154 
161  template<class Archive>
162  void load(Archive & archive, const unsigned int version)
163  {
164  // archive & mpMesh; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
165  // archive & mpIntracellularConductivityTensors; Loaded from HeartConfig every time constructor is called
166 
167  if (version >= 3)
168  {
169  archive & mHasPurkinje;
170  if (mHasPurkinje)
171  {
173  }
174  }
175  if (version >= 2)
176  {
177  archive & mExchangeHalos;
178  if (mExchangeHalos)
179  {
182  unsigned num_halo_nodes = mHaloNodes.size();
183  mHaloCellsDistributed.resize( num_halo_nodes );
184  for (unsigned local_index = 0; local_index < num_halo_nodes; local_index++)
185  {
186  unsigned global_index = mHaloNodes[local_index];
187  mHaloGlobalToLocalIndexMap[global_index] = local_index;
188  }
189  }
190  }
191 
192  // mCellsDistributed & mHaloCellsDistributed:
194 
195  // archive & mIionicCacheReplicated; // will be regenerated
196  // archive & mIntracellularStimulusCacheReplicated; // will be regenerated
197  archive & mDoCacheReplication;
198 
199  // we no longer have a bool mDoOneCacheReplication, but to maintain backwards compatibility
200  // we archive something if version==0
201  if (version==0)
202  {
203  bool do_one_cache_replication = true;
204  archive & do_one_cache_replication;
205  }
206 
208 
209  // Paranoia: check we agree with the mesh on who owns what
210  assert(mpDistributedVectorFactory == mpMesh->GetDistributedVectorFactory());
211  assert(mpDistributedVectorFactory->GetLow()==mpMesh->GetDistributedVectorFactory()->GetLow());
212  assert(mpDistributedVectorFactory->GetLocalOwnership()==mpMesh->GetDistributedVectorFactory()->GetLocalOwnership());
213  // archive & mMeshUnarchived; Not archived since set to true when archiving constructor is called.
214 
215  // not archiving mpConductivityModifier for the time being (mechanics simulations are only use-case at the moment, and they
216  // do not get archived...). mpConductivityModifier has to be reset to NULL upon load.
217  mpConductivityModifier = NULL;
218  }
219  BOOST_SERIALIZATION_SPLIT_MEMBER()
220 
221 
225 
226 protected:
227 
229  AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* mpMesh;
230 
234 
237 
241 
247 
253 
259 
265 
268 
278 
283 
290 
293 
304 
309 
315 
317  std::vector<unsigned> mHaloNodes;
318 
321 
323  std::map<unsigned, unsigned> mHaloGlobalToLocalIndexMap;
324 
330  std::vector<std::vector<unsigned> > mNodesToSendPerProcess;
331 
336  std::vector<std::vector<unsigned> > mNodesToReceivePerProcess;
337 
344 
351  void SetUpHaloCells(AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory);
352 
353 public:
364  AbstractCardiacTissue(AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory, bool exchangeHalos=false);
365 
371  AbstractCardiacTissue(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);
372 
374  virtual ~AbstractCardiacTissue();
375 
377  bool HasPurkinje();
378 
385  void SetCacheReplication(bool doCacheReplication);
386 
392  bool GetDoCacheReplication();
393 
397  const c_matrix<double, SPACE_DIM, SPACE_DIM>& rGetIntracellularConductivityTensor(unsigned elementIndex);
398 
406  virtual const c_matrix<double, SPACE_DIM, SPACE_DIM>& rGetExtracellularConductivityTensor(unsigned elementIndex);
407 
416  AbstractCardiacCellInterface* GetCardiacCell( unsigned globalIndex );
417 
426  AbstractCardiacCellInterface* GetPurkinjeCell( unsigned globalIndex );
427 
437 
447  virtual void SolveCellSystems(Vec existingSolution, double time, double nextTime, bool updateVoltage=false);
448 
451 
454 
457 
460 
461 
469  void UpdateCaches(unsigned globalIndex, unsigned localIndex, double nextTime);
470 
478  void UpdatePurkinjeCaches(unsigned globalIndex, unsigned localIndex, double nextTime);
479 
483  void ReplicateCaches();
484 
488  const std::vector<AbstractCardiacCellInterface*>& rGetCellsDistributed() const;
489 
493  const std::vector<AbstractCardiacCellInterface*>& rGetPurkinjeCellsDistributed() const;
494 
500  const AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pGetMesh() const;
501 
507  void SetConductivityModifier(AbstractConductivityModifier<ELEMENT_DIM,SPACE_DIM>* pModifier);
508 
520  template<class Archive>
521  void SaveCardiacCells(Archive & archive, const unsigned int version) const
522  {
523  const std::vector<AbstractCardiacCellInterface*> & r_cells_distributed = rGetCellsDistributed();
524  assert(mpDistributedVectorFactory == this->mpMesh->GetDistributedVectorFactory());
525  archive & mpDistributedVectorFactory; // Needed when loading
526  const unsigned num_cells = r_cells_distributed.size();
527  archive & num_cells;
528  for (unsigned i=0; i<num_cells; i++)
529  {
530  AbstractDynamicallyLoadableEntity* p_entity = dynamic_cast<AbstractDynamicallyLoadableEntity*>(r_cells_distributed[i]);
531  bool is_dynamic = (p_entity != NULL);
532  archive & is_dynamic;
533  if (is_dynamic)
534  {
535 #ifdef CHASTE_CAN_CHECKPOINT_DLLS
536  archive & p_entity->GetLoader()->GetLoadableModulePath();
537 #else
538  // We should have thrown an exception before this point
540 #endif // CHASTE_CAN_CHECKPOINT_DLLS
541  }
542  archive & r_cells_distributed[i];
543  if (mHasPurkinje)
544  {
545  archive & rGetPurkinjeCellsDistributed()[i];
546  }
547  }
548  }
549 
562  template<class Archive>
563  void LoadCardiacCells(Archive & archive, const unsigned int version)
564  {
565  // Note that p_factory loaded from this archive might not be the same as our mesh's factory,
566  // since we're loading a process-specific archive that could have been written by any process.
567  // We therefore need to use p_mesh_factory to determine the partitioning in use for the resumed
568  // simulation, and p_factory to determine what the original partitioning was when the simulation
569  // was saved.
570  DistributedVectorFactory* p_factory;
571  DistributedVectorFactory* p_mesh_factory = this->mpMesh->GetDistributedVectorFactory();
572  archive & p_factory;
573  unsigned num_cells;
574  archive & num_cells;
575  if (mCellsDistributed.empty())
576  {
577  mCellsDistributed.resize(p_mesh_factory->GetLocalOwnership());
578 #ifndef NDEBUG
579  // Paranoia
580  for (unsigned i=0; i<mCellsDistributed.size(); i++)
581  {
582  assert(mCellsDistributed[i] == NULL);
583  }
584 #endif
585  }
586  else
587  {
588  assert(mCellsDistributed.size() == p_mesh_factory->GetLocalOwnership());
589  }
590  if (mHasPurkinje)
591  {
592  if (mPurkinjeCellsDistributed.empty())
593  {
594  mPurkinjeCellsDistributed.resize(p_mesh_factory->GetLocalOwnership());
595  }
596  else
597  {
598  assert(mPurkinjeCellsDistributed.size() == p_mesh_factory->GetLocalOwnership());
599  }
600  }
601 
602  // We don't store a cell index in the archive, so need to work out what global index this tissue starts at.
603  // If we have an original factory we use the original low index; otherwise we use the current low index.
604  unsigned index_low = p_factory->GetOriginalFactory() ? p_factory->GetOriginalFactory()->GetLow() : p_mesh_factory->GetLow();
605 
606  // Track fake cells (which might have multiple pointers to the same object) to make sure we only delete non-local ones
607  std::set<FakeBathCell*> fake_cells_non_local, fake_cells_local;
608 
609  /*
610  * Historical note:
611  *
612  * We always do a dumb partition when we unarchive.
613  *
614  * When unarchive was first implemented in parallel (migration #1199) it was thought that we might want to repartition the mesh. This would be feasible and would give
615  * better partitions when we move to different numbers of processes. However it would require us to apply a new permutation to entire data structure.
616  *
617  * In the case where the original mesh was permuted and *copied* into the archive, we need to apply the stored permutation to the mesh but not to the archive (cells). That
618  * is, any permutation stored with the mesh can be safely ignored. (If we also had to repartition/permute the archive, then we would be applying a double permutation to the
619  * mesh and a single permutation to archive.)
620  *
621  */
622  for (unsigned local_index=0; local_index<num_cells; local_index++)
623  {
624  // Figure out where this cell goes
625  unsigned global_index = index_low + local_index;
626  bool local = p_mesh_factory->IsGlobalIndexLocal(global_index);
627 
628  // Check if this will be a halo cell
629  std::map<unsigned, unsigned>::const_iterator halo_position;
630  bool halo = ((halo_position=mHaloGlobalToLocalIndexMap.find(global_index)) != mHaloGlobalToLocalIndexMap.end());
631  // halo_position->second is local halo index
632 
633  bool is_dynamic;
634  archive & is_dynamic;
635  if (is_dynamic)
636  {
637 #ifdef CHASTE_CAN_CHECKPOINT_DLLS
638  // Ensure the shared object file for this cell model is loaded.
639  // We need to do this here, rather than in the class' serialization code,
640  // because that code won't be available until this is done...
641  std::string shared_object_path;
642  archive & shared_object_path;
643  DynamicModelLoaderRegistry::Instance()->GetLoader(shared_object_path);
644 #else
645  // Since checkpoints with dynamically loadable cells can only be
646  // created on Boost>=1.37, trying to load such a checkpoint on an
647  // earlier Boost would give an error when first opening the archive.
649 #endif // CHASTE_CAN_CHECKPOINT_DLLS
650  }
652  archive & p_cell;
653  AbstractCardiacCellInterface* p_purkinje_cell = NULL;
654  if (mHasPurkinje)
655  {
656  archive & p_purkinje_cell;
657  }
658  // Check if it's a fake cell
659  FakeBathCell* p_fake = dynamic_cast<FakeBathCell*>(p_cell);
660  if (p_fake)
661  {
662  if (halo || local)
663  {
664  fake_cells_local.insert(p_fake);
665  }
666  else
667  {
668  fake_cells_non_local.insert(p_fake);
669  }
670  }
671  FakeBathCell* p_fake_purkinje = dynamic_cast<FakeBathCell*>(p_purkinje_cell);
672  if (p_fake_purkinje)
673  {
674  if (halo || local)
675  {
676  fake_cells_local.insert(p_fake_purkinje);
677  }
678  else
679  {
680  fake_cells_non_local.insert(p_fake_purkinje);
681  }
682  }
683  // Add real cells to the local or halo vectors
684  if (local)
685  {
686  // Note that the original local_index was relative to the archived partition (distributed vector)
687  // The new_local_index is local relative to the new partition in memory
688  unsigned new_local_index = global_index - p_mesh_factory->GetLow();
689  assert(mCellsDistributed[new_local_index] == NULL);
690  mCellsDistributed[new_local_index] = p_cell;
691  if (mHasPurkinje)
692  {
693  assert(mPurkinjeCellsDistributed[new_local_index] == NULL);
694  mPurkinjeCellsDistributed[new_local_index] = p_purkinje_cell;
695  }
696  }
697  else if (halo)
698  {
699  assert(mHaloCellsDistributed[halo_position->second] == NULL);
700  mHaloCellsDistributed[halo_position->second] = p_cell;
701  }
702  else
703  {
704  if (!p_fake)
705  {
706  // Non-local real cell, so free the memory.
707  delete p_cell;
708  }
709  if (!p_fake_purkinje)
710  {
711  // This will be NULL if there's no Purkinje, so a delete is OK.
712  delete p_purkinje_cell;
713  }
714  }
715  }
716 
717  // Delete any unused fake cells
718  for (std::set<FakeBathCell*>::iterator it = fake_cells_non_local.begin();
719  it != fake_cells_non_local.end();
720  ++it)
721  {
722  if (fake_cells_local.find(*it) == fake_cells_local.end())
723  {
724  delete (*it);
725  }
726  }
727  }
728 };
729 
731 
732 namespace boost {
733 namespace serialization {
740 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
741 struct version<AbstractCardiacTissue<ELEMENT_DIM, SPACE_DIM> >
742 {
745 };
746 } // namespace serialization
747 } // namespace boost
748 
749 #endif /*ABSTRACTCARDIACTISSUE_HPP_*/
750 
ReplicatableVector mIntracellularStimulusCacheReplicated
#define TEMPLATED_CLASS_IS_ABSTRACT_2_UNSIGNED(T)
static DynamicModelLoaderRegistry * Instance()
AbstractConductivityTensors< ELEMENT_DIM, SPACE_DIM > * mpIntracellularConductivityTensors
void SetConductivityModifier(AbstractConductivityModifier< ELEMENT_DIM, SPACE_DIM > *pModifier)
DistributedVectorFactory * mpDistributedVectorFactory
const DynamicCellModelLoaderPtr GetLoader() const
ReplicatableVector & rGetPurkinjeIionicCacheReplicated()
std::vector< std::vector< unsigned > > mNodesToSendPerProcess
virtual void SolveCellSystems(Vec existingSolution, double time, double nextTime, bool updateVoltage=false)
static std::string GetMeshFilename()
FileFinder CopyTo(const FileFinder &rDest) const
Definition: FileFinder.cpp:308
void UpdateCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
void LoadCardiacCells(Archive &archive, const unsigned int version)
ReplicatableVector mPurkinjeIionicCacheReplicated
std::vector< AbstractCardiacCellInterface * > mPurkinjeCellsDistributed
#define NEVER_REACHED
Definition: Exception.hpp:206
ReplicatableVector mPurkinjeIntracellularStimulusCacheReplicated
void load(Archive &archive, const unsigned int version)
std::vector< AbstractCardiacCellInterface * > mHaloCellsDistributed
std::vector< AbstractCardiacCellInterface * > mCellsDistributed
bool IsGlobalIndexLocal(unsigned globalIndex)
static Archive * Get(void)
AbstractCardiacCellInterface * GetCardiacCellOrHaloCell(unsigned globalIndex)
ReplicatableVector & rGetIntracellularStimulusCacheReplicated()
std::vector< std::vector< unsigned > > mNodesToReceivePerProcess
AbstractCardiacCellInterface * GetPurkinjeCell(unsigned globalIndex)
void SetUpHaloCells(AbstractCardiacCellFactory< ELEMENT_DIM, SPACE_DIM > *pCellFactory)
void SaveCardiacCells(Archive &archive, const unsigned int version) const
#define CHASTE_VERSION_CONTENT(N)
DynamicCellModelLoaderPtr GetLoader(const std::string &rPath)
ReplicatableVector mIionicCacheReplicated
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
const std::vector< AbstractCardiacCellInterface * > & rGetCellsDistributed() const
virtual const c_matrix< double, SPACE_DIM, SPACE_DIM > & rGetExtracellularConductivityTensor(unsigned elementIndex)
AbstractConductivityModifier< ELEMENT_DIM, SPACE_DIM > * mpConductivityModifier
std::map< unsigned, unsigned > mHaloGlobalToLocalIndexMap
bool Exists() const
Definition: FileFinder.cpp:180
void Resize(unsigned size)
const std::vector< AbstractCardiacCellInterface * > & rGetPurkinjeCellsDistributed() const
void SetCacheReplication(bool doCacheReplication)
Forward declaration which is going to be used for friendship.
ReplicatableVector & rGetPurkinjeIntracellularStimulusCacheReplicated()
const c_matrix< double, SPACE_DIM, SPACE_DIM > & rGetIntracellularConductivityTensor(unsigned elementIndex)
static std::string GetArchiveRelativePath()
std::vector< unsigned > mHaloNodes
AbstractCardiacCellInterface * GetCardiacCell(unsigned globalIndex)
void UpdatePurkinjeCaches(unsigned globalIndex, unsigned localIndex, double nextTime)
#define TRY_IF_MASTER(method)
Definition: PetscTools.hpp:90
static HeartConfig * Instance()
const AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * pGetMesh() const
void save(Archive &archive, const unsigned int version) const
ReplicatableVector & rGetIionicCacheReplicated()