AbstractCardiacPde.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 ABSTRACTCARDIACPDE_HPP_
#define ABSTRACTCARDIACPDE_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 "AbstractCardiacCell.hpp"
#include "AbstractCardiacCellFactory.hpp"
#include "AbstractConductivityTensors.hpp"

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

// IMPORTANT NOTE: the inheritance of AbstractPde has to be 'virtual'
// ie "class AbstractCardiacPde : public virtual AbstractPde"
// because AbstractPde will be the top class in a 'dreaded diamond':
//      A
//     / \     A = AbstractPde, B = AbstractCardiac, C = AbstractLinearParabolic (etc)
//    B   C    D = MonodomainPde
//     \ /
//      D
//
// B and C must use virtual inheritence of A in order for D to only contain 1 instance
// of the member variables in A

template <unsigned ELEMENT_DIM, unsigned SPACE_DIM = ELEMENT_DIM>
class AbstractCardiacPde
{
private:

    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive & archive, const unsigned int version)
    {
        // archive & mpMesh; Archived in save/load_constructs at the bottom of mono/bidomainPde.hpp
        // archive & mpIntracellularConductivityTensors; Loaded from HeartConfig every time constructor is called
        // archive & mCellsDistributed; Archived in save/load_constructs at the bottom of mono/bidomainPde.hpp
        // archive & mIionicCacheReplicated; // will be regenerated
        // archive & mIntracellularStimulusCacheReplicated; // will be regenerated
        // archive & mStride; // archiving constructor sets this.
        archive & mDoCacheReplication;
        archive & mDoOneCacheReplication;
        (*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.
    }

    void CreateIntracellularConductivityTensor();

protected:

    AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* mpMesh;

    AbstractConductivityTensors<SPACE_DIM> *mpIntracellularConductivityTensors;

    std::vector< AbstractCardiacCell* > mCellsDistributed;

    ReplicatableVector mIionicCacheReplicated;

    ReplicatableVector mIntracellularStimulusCacheReplicated;

    const unsigned mStride;

    HeartConfig* mpConfig;

    bool mDoCacheReplication;

    bool mDoOneCacheReplication;

    DistributedVectorFactory* mpDistributedVectorFactory;

    bool mMeshUnarchived;

    void DeleteCells(bool deleteFakeCells);

public:
    AbstractCardiacPde(AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* pCellFactory,
                       const unsigned stride=1);

    AbstractCardiacPde(std::vector<AbstractCardiacCell*>& rCellsDistributed,
                       AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
                       const unsigned stride);

    virtual ~AbstractCardiacPde();

    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 collection of cells 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 collection
                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(AbstractCardiacPde);

#endif /*ABSTRACTCARDIACPDE_HPP_*/