AbstractCardiacTissue.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 #ifndef ABSTRACTCARDIACTISSUE_HPP_
00029 #define ABSTRACTCARDIACTISSUE_HPP_
00030 
00031 #include <set>
00032 #include <vector>
00033 #include <boost/shared_ptr.hpp>
00034 
00035 #include "UblasMatrixInclude.hpp"
00036 
00037 #include "ChasteSerialization.hpp"
00038 #include "ClassIsAbstract.hpp"
00039 #include "ChasteSerializationVersion.hpp"
00040 #include <boost/serialization/base_object.hpp>
00041 #include <boost/serialization/shared_ptr.hpp>
00042 #include <boost/serialization/vector.hpp>
00043 #include <boost/serialization/string.hpp>
00044 #include <boost/serialization/split_member.hpp>
00045 
00046 #include "AbstractCardiacCell.hpp"
00047 #include "FakeBathCell.hpp"
00048 #include "AbstractCardiacCellFactory.hpp"
00049 #include "AbstractConductivityTensors.hpp"
00050 
00051 #include "ReplicatableVector.hpp"
00052 #include "HeartConfig.hpp"
00053 #include "ArchiveLocationInfo.hpp"
00054 #include "AbstractDynamicallyLoadableEntity.hpp"
00055 #include "DynamicModelLoaderRegistry.hpp"
00056 #include "AbstractConductivityModifier.hpp"
00057 
00070 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM = ELEMENT_DIM>
00071 class AbstractCardiacTissue : boost::noncopyable
00072 {
00073 private:
00074 
00076     friend class boost::serialization::access;
00083     template<class Archive>
00084     void save(Archive & archive, const unsigned int version) const
00085     {
00086         if (version >= 2)
00087         {
00088             archive & mExchangeHalos;
00089         }
00090         // Don't use the std::vector serialization for cardiac cells, so that we can load them
00091         // more cleverly when migrating checkpoints.
00092         SaveCardiacCells(*ProcessSpecificArchive<Archive>::Get(), version);
00093 
00094         // archive & mpMesh; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
00095         // archive & mpIntracellularConductivityTensors; Loaded from HeartConfig every time constructor is called
00096         if (HeartConfig::Instance()->IsMeshProvided() && HeartConfig::Instance()->GetLoadMesh())
00097         {
00098             switch (HeartConfig::Instance()->GetConductivityMedia())
00099             {
00100                 case cp::media_type::Orthotropic:
00101                 {
00102                     FileFinder source_file(mFibreFilePathNoExtension + ".ortho", RelativeTo::AbsoluteOrCwd);
00103                     assert(source_file.Exists());
00104                     FileFinder dest_file(ArchiveLocationInfo::GetArchiveRelativePath() + ArchiveLocationInfo::GetMeshFilename() + ".ortho", RelativeTo::ChasteTestOutput);
00105 
00106                     if (PetscTools::AmMaster())
00107                     {
00108                         MPIABORTIFNON0(system,"cp " + source_file.GetAbsolutePath() + " " + dest_file.GetAbsolutePath());
00109                     }
00110                     PetscTools::Barrier();
00111                     break;
00112                 }
00113 
00114                 case cp::media_type::Axisymmetric:
00115                 {
00116                     FileFinder source_file(mFibreFilePathNoExtension + ".axi", RelativeTo::AbsoluteOrCwd);
00117                     assert(source_file.Exists());
00118                     FileFinder dest_file(ArchiveLocationInfo::GetArchiveRelativePath() + ArchiveLocationInfo::GetMeshFilename() + ".axi", RelativeTo::ChasteTestOutput);
00119 
00120                     if (PetscTools::AmMaster())
00121                     {
00122                         MPIABORTIFNON0(system,"cp " + source_file.GetAbsolutePath() + " " + dest_file.GetAbsolutePath());
00123                     }
00124                     PetscTools::Barrier();
00125 
00126                     break;
00127                 }
00128 
00129                 case cp::media_type::NoFibreOrientation:
00130                     break;
00131 
00132                 default :
00133                     NEVER_REACHED;
00134             }
00135         }
00136 
00137         // archive & mCellsDistributed; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
00138         // archive & mIionicCacheReplicated; // will be regenerated
00139         // archive & mIntracellularStimulusCacheReplicated; // will be regenerated
00140         archive & mDoCacheReplication;
00141 
00142         (*ProcessSpecificArchive<Archive>::Get()) & mpDistributedVectorFactory;
00143 
00144         // Paranoia: check we agree with the mesh on who owns what
00145         assert(mpDistributedVectorFactory->GetLow()==mpMesh->GetDistributedVectorFactory()->GetLow());
00146         assert(mpDistributedVectorFactory->GetLocalOwnership()==mpMesh->GetDistributedVectorFactory()->GetLocalOwnership());
00147         // archive & mMeshUnarchived; Not archived since set to true when archiving constructor is called.
00148     }
00149 
00156     template<class Archive>
00157     void load(Archive & archive, const unsigned int version)
00158     {
00159         // archive & mpMesh; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
00160         // archive & mpIntracellularConductivityTensors; Loaded from HeartConfig every time constructor is called
00161 
00162         if (version >= 2)
00163         {
00164             archive & mExchangeHalos;
00165             if (mExchangeHalos)
00166             {
00167                 mpMesh->CalculateNodeExchange(mNodesToSendPerProcess, mNodesToReceivePerProcess);
00168                 CalculateHaloNodesFromNodeExchange();
00169                 unsigned num_halo_nodes = mHaloNodes.size();
00170                 mHaloCellsDistributed.resize( num_halo_nodes );
00171                 for (unsigned local_index = 0; local_index < num_halo_nodes; local_index++)
00172                 {
00173                     unsigned global_index = mHaloNodes[local_index];
00174                     mHaloGlobalToLocalIndexMap[global_index] = local_index;
00175                 }
00176             }
00177         }
00178 
00179         // mCellsDistributed & mHaloCellsDistributed:
00180         LoadCardiacCells(*ProcessSpecificArchive<Archive>::Get(), version);
00181 
00182         // archive & mIionicCacheReplicated; // will be regenerated
00183         // archive & mIntracellularStimulusCacheReplicated; // will be regenerated
00184         archive & mDoCacheReplication;
00185 
00186         // we no longer have a bool mDoOneCacheReplication, but to maintain backwards compatibility
00187         // we archive something if version==0
00188         if (version==0)
00189         {
00190             bool do_one_cache_replication = true;
00191             archive & do_one_cache_replication;
00192         }
00193 
00194         (*ProcessSpecificArchive<Archive>::Get()) & mpDistributedVectorFactory;
00195 
00196         // Paranoia: check we agree with the mesh on who owns what
00197         assert(mpDistributedVectorFactory->GetLow()==mpMesh->GetDistributedVectorFactory()->GetLow());
00198         assert(mpDistributedVectorFactory->GetLocalOwnership()==mpMesh->GetDistributedVectorFactory()->GetLocalOwnership());
00199         // archive & mMeshUnarchived; Not archived since set to true when archiving constructor is called.
00200 
00201         // not archiving mpConductivityModifier for the time being (mechanics simulations are only use-case at the moment, and they
00202         // do not get archived...). mpConductivityModifier has to be reset to NULL upon load.
00203         mpConductivityModifier = NULL;
00204     }
00205     BOOST_SERIALIZATION_SPLIT_MEMBER()
00206 
00207     
00210     void CreateIntracellularConductivityTensor();
00211 
00212 protected:
00213 
00215     AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* mpMesh;
00216 
00219     AbstractConductivityTensors<ELEMENT_DIM,SPACE_DIM>* mpIntracellularConductivityTensors;
00220 
00222     std::vector< AbstractCardiacCell* > mCellsDistributed;
00223 
00228     ReplicatableVector mIionicCacheReplicated;
00229 
00234     ReplicatableVector mIntracellularStimulusCacheReplicated;
00235 
00237     HeartConfig* mpConfig;
00238 
00248     bool mDoCacheReplication;
00249 
00258     DistributedVectorFactory* mpDistributedVectorFactory;
00259 
00263     bool mMeshUnarchived;
00264 
00268     std::string mFibreFilePathNoExtension;
00269 
00275     AbstractConductivityModifier<ELEMENT_DIM,SPACE_DIM>* mpConductivityModifier;
00276 
00281     bool mExchangeHalos;
00282 
00284     std::vector<unsigned> mHaloNodes;
00285 
00287     std::vector< AbstractCardiacCell* > mHaloCellsDistributed;
00288 
00290     std::map<unsigned, unsigned> mHaloGlobalToLocalIndexMap;
00291 
00297     std::vector<std::vector<unsigned> > mNodesToSendPerProcess;
00298 
00303     std::vector<std::vector<unsigned> > mNodesToReceivePerProcess;
00304 
00310     void CalculateHaloNodesFromNodeExchange();
00311 
00318     void SetUpHaloCells(AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory);
00319 
00320 public:
00330     AbstractCardiacTissue(AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory, bool exchangeHalos=false);
00331 
00337     AbstractCardiacTissue(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);
00338 
00340     virtual ~AbstractCardiacTissue();
00341 
00348     void SetCacheReplication(bool doCacheReplication);
00349 
00355     bool GetDoCacheReplication();
00356 
00360     const c_matrix<double, SPACE_DIM, SPACE_DIM>& rGetIntracellularConductivityTensor(unsigned elementIndex);
00361 
00369     virtual const c_matrix<double, SPACE_DIM, SPACE_DIM>& rGetExtracellularConductivityTensor(unsigned elementIndex);
00370 
00379     AbstractCardiacCell* GetCardiacCell( unsigned globalIndex );
00380 
00389     AbstractCardiacCell* GetCardiacCellOrHaloCell( unsigned globalIndex );
00390 
00400     virtual void SolveCellSystems(Vec existingSolution, double time, double nextTime, bool updateVoltage=false);
00401 
00403     ReplicatableVector& rGetIionicCacheReplicated();
00404 
00406     ReplicatableVector& rGetIntracellularStimulusCacheReplicated();
00407 
00408 
00416     void UpdateCaches(unsigned globalIndex, unsigned localIndex, double nextTime);
00417 
00421     void ReplicateCaches();
00422 
00426     const std::vector<AbstractCardiacCell*>& rGetCellsDistributed() const;
00427 
00433     const AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pGetMesh() const;
00434 
00440     void SetConductivityModifier(AbstractConductivityModifier<ELEMENT_DIM,SPACE_DIM>* pModifier);
00441 
00453     template<class Archive>
00454     void SaveCardiacCells(Archive & archive, const unsigned int version) const
00455     {
00456         const std::vector<AbstractCardiacCell*> & r_cells_distributed = rGetCellsDistributed();
00457         archive & mpDistributedVectorFactory; // Needed when loading
00458         const unsigned num_cells = r_cells_distributed.size();
00459         archive & num_cells;
00460         for (unsigned i=0; i<num_cells; i++)
00461         {
00462             AbstractDynamicallyLoadableEntity* p_entity = dynamic_cast<AbstractDynamicallyLoadableEntity*>(r_cells_distributed[i]);
00463             bool is_dynamic = (p_entity != NULL);
00464             archive & is_dynamic;
00465             if (is_dynamic)
00466             {
00467 #ifdef CHASTE_CAN_CHECKPOINT_DLLS
00468                 archive & p_entity->GetLoader()->GetLoadableModulePath();
00469 #else
00470                 // We should have thrown an exception before this point
00471                 NEVER_REACHED;
00472 #endif // CHASTE_CAN_CHECKPOINT_DLLS
00473             }
00474             archive & r_cells_distributed[i];
00475         }
00476     }
00477 
00490     template<class Archive>
00491     void LoadCardiacCells(Archive & archive, const unsigned int version)
00492     {
00493         DistributedVectorFactory* p_factory;
00494         DistributedVectorFactory* p_mesh_factory = this->mpMesh->GetDistributedVectorFactory();
00495         archive & p_factory;
00496         unsigned num_cells;
00497         archive & num_cells;
00498         if (mCellsDistributed.empty())
00499         {
00500             mCellsDistributed.resize(p_factory->GetLocalOwnership());
00501 #ifndef NDEBUG
00502             // Paranoia
00503             for (unsigned i=0; i<mCellsDistributed.size(); i++)
00504             {
00505                 assert(mCellsDistributed[i] == NULL);
00506             }
00507 #endif
00508         }
00509         else
00510         {
00511             assert(mCellsDistributed.size() == p_mesh_factory->GetLocalOwnership());
00512         }
00513 
00514         // We don't store a cell index in the archive, so need to work out what global index this tissue starts at.
00515         // If we have an original factory we use the original low index; otherwise we use the current low index.
00516         unsigned index_low = p_factory->GetOriginalFactory() ? p_factory->GetOriginalFactory()->GetLow() : p_mesh_factory->GetLow();
00517 
00518         // Track fake bath cells to make sure we only delete non-local ones
00519         std::set<FakeBathCell*> fake_bath_cells_non_local, fake_bath_cells_local;
00520 
00521         const std::vector<unsigned>& r_permutation = this->mpMesh->rGetNodePermutation();
00522         for (unsigned local_index=0; local_index<num_cells; local_index++)
00523         {
00524             // If we're permuting, figure out where this cell goes
00525             unsigned original_global_index = index_low + local_index;
00526             unsigned new_global_index;
00527             if (r_permutation.empty())
00528             {
00529                 new_global_index = original_global_index;
00530             }
00531             else
00532             {
00534 //                new_global_index = r_permutation[original_global_index];
00535                 NEVER_REACHED;
00536             }
00537             unsigned new_local_index = new_global_index - p_mesh_factory->GetLow();
00538             bool local = p_mesh_factory->IsGlobalIndexLocal(new_global_index);
00539 
00540             // Check if this will be a halo cell
00541             std::map<unsigned, unsigned>::const_iterator halo_position;
00542             bool halo = ((halo_position=mHaloGlobalToLocalIndexMap.find(new_global_index)) != mHaloGlobalToLocalIndexMap.end());
00543             // halo_position->second is local halo index
00544 
00545             bool is_dynamic;
00546             archive & is_dynamic;
00547             if (is_dynamic)
00548             {
00549 #ifdef CHASTE_CAN_CHECKPOINT_DLLS
00550                 // Ensure the shared object file for this cell model is loaded.
00551                 // We need to do this here, rather than in the class' serialization code,
00552                 // because that code won't be available until this is done...
00553                 std::string shared_object_path;
00554                 archive & shared_object_path;
00555                 DynamicModelLoaderRegistry::Instance()->GetLoader(shared_object_path);
00556 #else
00557                 // Since checkpoints with dynamically loadable cells can only be
00558                 // created on Boost>=1.37, trying to load such a checkpoint on an
00559                 // earlier Boost would give an error when first opening the archive.
00560                 NEVER_REACHED;
00561 #endif // CHASTE_CAN_CHECKPOINT_DLLS
00562             }
00563             AbstractCardiacCell* p_cell;
00564             archive & p_cell;
00565             // Check if it's a fake cell
00566             FakeBathCell* p_fake = dynamic_cast<FakeBathCell*>(p_cell);
00567             if (p_fake)
00568             {
00569                 if (halo || local)
00570                 {
00571                     fake_bath_cells_local.insert(p_fake);
00572                 }
00573                 else
00574                 {
00575                     fake_bath_cells_non_local.insert(p_fake);
00576                 }
00577             }
00578             // Add real cells to the local or halo vectors
00579             if (local)
00580             {
00581                 assert(mCellsDistributed[new_local_index] == NULL);
00582                 mCellsDistributed[new_local_index] = p_cell;
00583             }
00584             else if (halo)
00585             {
00586                 assert(mHaloCellsDistributed[halo_position->second] == NULL);
00587                 mHaloCellsDistributed[halo_position->second] = p_cell;
00588             }
00589             else if (!p_fake)
00590             {
00591                 // Non-local real cell, so free the memory.
00592                 delete p_cell;
00593             }
00594         }
00595 
00596         // Delete any unused fake cells
00597         if (!fake_bath_cells_non_local.empty())
00598         {
00599             for (std::set<FakeBathCell*>::iterator it = fake_bath_cells_non_local.begin();
00600                  it != fake_bath_cells_non_local.end();
00601                  ++it)
00602             {
00603                 if (fake_bath_cells_local.find(*it) == fake_bath_cells_local.end())
00604                 {
00605                     delete (*it);
00606                 }
00607             }
00608         }
00609     }
00610 };
00611 
00612 TEMPLATED_CLASS_IS_ABSTRACT_2_UNSIGNED(AbstractCardiacTissue)
00613 
00614 namespace boost {
00615 namespace serialization {
00622 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00623 struct version<AbstractCardiacTissue<ELEMENT_DIM, SPACE_DIM> >
00624 {
00626     CHASTE_VERSION_CONTENT(2);
00627 };
00628 } // namespace serialization
00629 } // namespace boost
00630 
00631 #endif /*ABSTRACTCARDIACTISSUE_HPP_*/
00632 

Generated on Tue May 31 14:31:42 2011 for Chaste by  doxygen 1.5.5