AbstractCardiacTissue.hpp

/*

Copyright (C) University of Oxford, 2005-2010

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Chaste is free software: you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 2.1 of the License, or
(at your option) any later version.

Chaste is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details. The offer of Chaste under the terms of the
License is subject to the License being interpreted in accordance with
English Law and subject to any action against the University of Oxford
being under the jurisdiction of the English Courts.

You should have received a copy of the GNU Lesser General Public License
along with Chaste. If not, see <http://www.gnu.org/licenses/>.

*/
#ifndef ABSTRACTCARDIACTISSUE_HPP_
#define ABSTRACTCARDIACTISSUE_HPP_

#include <set>
#include <vector>
#include "ChasteSerialization.hpp"
#include "ClassIsAbstract.hpp"
#include <boost/serialization/base_object.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/serialization/shared_ptr.hpp>
#include <boost/serialization/vector.hpp>
#include <boost/serialization/string.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/serialization/version.hpp>
#include <boost/serialization/split_member.hpp>

#include "AbstractCardiacCell.hpp"
#include "FakeBathCell.hpp"
#include "AbstractCardiacCellFactory.hpp"
#include "AbstractConductivityTensors.hpp"

#include "ReplicatableVector.hpp"
#include "HeartConfig.hpp"
#include "ArchiveLocationInfo.hpp"
#include "AbstractDynamicallyLoadableEntity.hpp"
#include "DynamicModelLoaderRegistry.hpp"

template <unsigned ELEMENT_DIM, unsigned SPACE_DIM = ELEMENT_DIM>
class AbstractCardiacTissue : boost::noncopyable
{
private:

    friend class boost::serialization::access;
    template<class Archive>
    void save(Archive & archive, const unsigned int version) const
    {
        // archive & mpMesh; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
        // archive & mpIntracellularConductivityTensors; Loaded from HeartConfig every time constructor is called
        if (HeartConfig::Instance()->IsMeshProvided() && HeartConfig::Instance()->GetLoadMesh())
        {
            switch (HeartConfig::Instance()->GetConductivityMedia())
            {
                case cp::media_type::Orthotropic:
                {
                    FileFinder source_file(mFibreFilePathNoExtension + ".ortho", RelativeTo::AbsoluteOrCwd);
                    assert(source_file.Exists());
                    FileFinder dest_file(ArchiveLocationInfo::GetArchiveRelativePath() + ArchiveLocationInfo::GetMeshFilename() + ".ortho", RelativeTo::ChasteTestOutput);                    
                                       
                    if (PetscTools::AmMaster())
                    {
                        MPIABORTIFNON0(system,"cp " + source_file.GetAbsolutePath() + " " + dest_file.GetAbsolutePath());
                    }
                    PetscTools::Barrier();
                    break;
                }

                case cp::media_type::Axisymmetric:
                {
                    FileFinder source_file(mFibreFilePathNoExtension + ".axi", RelativeTo::AbsoluteOrCwd);
                    assert(source_file.Exists());
                    FileFinder dest_file(ArchiveLocationInfo::GetArchiveRelativePath() + ArchiveLocationInfo::GetMeshFilename() + ".axi", RelativeTo::ChasteTestOutput);

                    if (PetscTools::AmMaster())
                    {
                        MPIABORTIFNON0(system,"cp " + source_file.GetAbsolutePath() + " " + dest_file.GetAbsolutePath());
                    }
                    PetscTools::Barrier();
                    
                    break;
                }

                case cp::media_type::NoFibreOrientation:
                    break;

                default :
                    NEVER_REACHED;

            }
        }

        // archive & mCellsDistributed; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
        // archive & mIionicCacheReplicated; // will be regenerated
        // archive & mIntracellularStimulusCacheReplicated; // will be regenerated
        archive & mDoCacheReplication;

        (*ProcessSpecificArchive<Archive>::Get()) & mpDistributedVectorFactory;

        // Paranoia: check we agree with the mesh on who owns what
        assert(mpDistributedVectorFactory->GetLow()==mpMesh->GetDistributedVectorFactory()->GetLow());
        assert(mpDistributedVectorFactory->GetLocalOwnership()==mpMesh->GetDistributedVectorFactory()->GetLocalOwnership());
        // archive & mMeshUnarchived; Not archived since set to true when archiving constructor is called.
    }

    template<class Archive>
    void load(Archive & archive, const unsigned int version)
    {
        // archive & mpMesh; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
        // archive & mpIntracellularConductivityTensors; Loaded from HeartConfig every time constructor is called
        // archive & mCellsDistributed; Archived in save/load_constructs at the bottom of Mono/BidomainTissue.hpp
        // archive & mIionicCacheReplicated; // will be regenerated
        // archive & mIntracellularStimulusCacheReplicated; // will be regenerated
        archive & mDoCacheReplication;

        // we no longer have a bool mDoOneCacheReplication, but to maintain backwards compatibility
        // we archive something if version==0
        if(version==0)
        {
            bool do_one_cache_replication = true;
            archive & do_one_cache_replication;
        }

        (*ProcessSpecificArchive<Archive>::Get()) & mpDistributedVectorFactory;

        // Paranoia: check we agree with the mesh on who owns what
        assert(mpDistributedVectorFactory->GetLow()==mpMesh->GetDistributedVectorFactory()->GetLow());
        assert(mpDistributedVectorFactory->GetLocalOwnership()==mpMesh->GetDistributedVectorFactory()->GetLocalOwnership());
        // archive & mMeshUnarchived; Not archived since set to true when archiving constructor is called.
    }
    BOOST_SERIALIZATION_SPLIT_MEMBER()

    
    void CreateIntracellularConductivityTensor();

protected:

    AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* mpMesh;

    AbstractConductivityTensors<ELEMENT_DIM,SPACE_DIM>* mpIntracellularConductivityTensors;

    std::vector< AbstractCardiacCell* > mCellsDistributed;

    ReplicatableVector mIionicCacheReplicated;

    ReplicatableVector mIntracellularStimulusCacheReplicated;

    HeartConfig* mpConfig;

    bool mDoCacheReplication;

    DistributedVectorFactory* mpDistributedVectorFactory;

    bool mMeshUnarchived;
    
    std::string mFibreFilePathNoExtension;

public:
    AbstractCardiacTissue(AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory);

    AbstractCardiacTissue(std::vector<AbstractCardiacCell*>& rCellsDistributed,
                          AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);

    virtual ~AbstractCardiacTissue();

    void MergeCells(const std::vector<AbstractCardiacCell*>& rOtherCells);

    void SetCacheReplication(bool doCacheReplication);

    bool GetDoCacheReplication();

    const c_matrix<double, SPACE_DIM, SPACE_DIM>& rGetIntracellularConductivityTensor(unsigned elementIndex);

    AbstractCardiacCell* GetCardiacCell( unsigned globalIndex );

    virtual void SolveCellSystems(Vec existingSolution, double time, double nextTime);

    ReplicatableVector& rGetIionicCacheReplicated();

    ReplicatableVector& rGetIntracellularStimulusCacheReplicated();


    void UpdateCaches(unsigned globalIndex, unsigned localIndex, double nextTime);

    void ReplicateCaches();

    const std::vector<AbstractCardiacCell*>& rGetCellsDistributed() const;

    const AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pGetMesh() const;

    template<class Archive>
    void SaveCardiacCells(Archive & archive, const unsigned int version) const
    {
        Archive& r_archive = *ProcessSpecificArchive<Archive>::Get();
        const std::vector<AbstractCardiacCell*> & r_cells_distributed = rGetCellsDistributed();
        r_archive & mpDistributedVectorFactory; // Needed when loading
        const unsigned num_cells = r_cells_distributed.size();
        r_archive & num_cells;
        for (unsigned i=0; i<num_cells; i++)
        {
            AbstractDynamicallyLoadableEntity* p_entity = dynamic_cast<AbstractDynamicallyLoadableEntity*>(r_cells_distributed[i]);
            bool is_dynamic = (p_entity != NULL);
            r_archive & is_dynamic;
            if (is_dynamic)
            {
#ifdef CHASTE_CAN_CHECKPOINT_DLLS
                r_archive & p_entity->GetLoader()->GetLoadableModulePath();
#else
                // We should have thrown an exception before this point
                NEVER_REACHED;
#endif // CHASTE_CAN_CHECKPOINT_DLLS
            }
            r_archive & r_cells_distributed[i];
        }
    }

    template<class Archive>
    static void LoadCardiacCells(Archive & archive, const unsigned int version,
                                 std::vector<AbstractCardiacCell*>& rCells,
                                 AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
    {
        DistributedVectorFactory* p_factory;
        archive & p_factory;
        unsigned num_cells;
        archive & num_cells;
        rCells.resize(p_factory->GetLocalOwnership());
#ifndef NDEBUG
        // Paranoia
        for (unsigned i=0; i<rCells.size(); i++)
        {
            assert(rCells[i] == NULL);
        }
#endif

        // We don't store a cell index in the archive, so need to work out what global
        // index this tissue starts up.  If we're migrating (so have an
        // original factory) we use the original low index; otherwise we use the current
        // low index.
        unsigned index_low = p_factory->GetOriginalFactory() ? p_factory->GetOriginalFactory()->GetLow() : p_factory->GetLow();

        // Track fake bath cells to make sure we only delete non-local ones
        std::set<FakeBathCell*> fake_bath_cells_non_local, fake_bath_cells_local;

        for (unsigned local_index=0; local_index<num_cells; local_index++)
        {
            // If we're permuting, figure out where this cell goes
            unsigned original_global_index = index_low + local_index;
            const std::vector<unsigned>& r_permutation = pMesh->rGetNodePermutation();
            unsigned new_global_index;
            if (r_permutation.empty())
            {
                new_global_index = original_global_index;
            }
            else
            {
//                new_global_index = r_permutation[original_global_index];
                NEVER_REACHED;
            }
            unsigned new_local_index = new_global_index - p_factory->GetLow();
            bool local = p_factory->IsGlobalIndexLocal(new_global_index);

            bool is_dynamic;
            archive & is_dynamic;
            if (is_dynamic)
            {
#ifdef CHASTE_CAN_CHECKPOINT_DLLS
                // Ensure the shared object file for this cell model is loaded.
                // We need to do this here, rather than in the class' serialization code,
                // because that code won't be available until this is done...
                std::string shared_object_path;
                archive & shared_object_path;
                DynamicModelLoaderRegistry::Instance()->GetLoader(shared_object_path);
#else
                // Since checkpoints with dynamically loadable cells can only be
                // created on Boost>=1.37, trying to load such a checkpoint on an
                // earlier Boost would give an error when first opening the archive.
                NEVER_REACHED;
#endif // CHASTE_CAN_CHECKPOINT_DLLS
            }
            AbstractCardiacCell* p_cell;
            archive & p_cell;
            // Check if it's a fake cell
            FakeBathCell* p_fake = dynamic_cast<FakeBathCell*>(p_cell);
            if (local)
            {
                rCells[new_local_index] = p_cell; // Add to local cells 
                if (p_fake)
                {
                    fake_bath_cells_local.insert(p_fake);
                }
            }
            else
            {
                if (p_fake)
                {
                    fake_bath_cells_non_local.insert(p_fake);
                }
                else
                {
                    // Non-local real cell, so free the memory.
                    delete p_cell;
                }
            }
        }

        // Delete any unused fake cells
        if (!fake_bath_cells_non_local.empty())
        {
            for (std::set<FakeBathCell*>::iterator it = fake_bath_cells_non_local.begin();
                 it != fake_bath_cells_non_local.end();
                 ++it)
            {
                if (fake_bath_cells_local.find(*it) == fake_bath_cells_local.end())
                {
                    delete (*it);
                }
            }
        }
    }
};

TEMPLATED_CLASS_IS_ABSTRACT_2_UNSIGNED(AbstractCardiacTissue)

namespace boost {
namespace serialization {
template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
struct version<AbstractCardiacTissue<ELEMENT_DIM, SPACE_DIM> >
{
    BOOST_STATIC_CONSTANT(unsigned, value = 1);
};
} // namespace serialization
} // namespace boost

#endif /*ABSTRACTCARDIACTISSUE_HPP_*/

---