AbstractCardiacProblem.hpp

/*

Copyright (C) University of Oxford, 2005-2011

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 ABSTRACTCARDIACPROBLEM_HPP_
#define ABSTRACTCARDIACPROBLEM_HPP_



#include <string>
#include <vector>
#include <cassert>
#include <climits>
#include <boost/shared_ptr.hpp>

#include "ChasteSerialization.hpp"
#include <boost/serialization/vector.hpp>
#include <boost/serialization/string.hpp>
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/shared_ptr.hpp>
#include "ClassIsAbstract.hpp"
#include "ChasteSerializationVersion.hpp"

#include "AbstractTetrahedralMesh.hpp"
#include "AbstractCardiacCell.hpp"
#include "AbstractCardiacCellFactory.hpp"
#include "AbstractCardiacTissue.hpp"
#include "AbstractDynamicLinearPdeSolver.hpp"
#include "BoundaryConditionsContainer.hpp"
#include "DistributedVectorFactory.hpp"
#include "Hdf5DataReader.hpp"
#include "Hdf5DataWriter.hpp"

/*
 * Archiving extravaganza:
 *
 * We archive mesh and tissue through a pointer to an abstract class.  All the potential concrete
 * classes need to be included here, so they are registered with boost.  If not, boost won't be
 * able to find the archiving methods of the concrete class and will throw the following
 * exception:
 *
 *       terminate called after throwing an instance of 'boost::archive::archive_exception'
 *       what():  unregistered class
 *
 * No member variable is defined to be of any of these clases, so removing them won't
 * produce any compiler error.  The exception above will occur at runtime.
 *
 * This might not be even necessary in certain cases, if the file is included implicitly by another
 * header file or by the test itself.  It's safer though.
 */
#include "DistributedTetrahedralMesh.hpp"
#include "TetrahedralMesh.hpp" //May be needed for unarchiving a mesh
#include "MonodomainTissue.hpp"
#include "BidomainTissue.hpp"

class AbstractUntemplatedCardiacProblem : boost::noncopyable
{
public:
    virtual ~AbstractUntemplatedCardiacProblem()
    {}
};

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
class AbstractCardiacProblem : public AbstractUntemplatedCardiacProblem
{
    friend class TestBidomainWithBath;
    friend class TestCardiacSimulationArchiver;

    typedef typename boost::shared_ptr<BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM> >
        BccType;

private:
    friend class boost::serialization::access;

    template<class Archive>
    void save(Archive & archive, const unsigned int version) const
    {
        if (version >= 1)
        {
            const unsigned element_dim=ELEMENT_DIM;
            archive & element_dim;
            const unsigned space_dim=SPACE_DIM;
            archive & space_dim;
            const unsigned problem_dim=PROBLEM_DIM;
            archive & problem_dim;
        }
        archive & mMeshFilename;
        archive & mpMesh;
        //archive & mAllocatedMemoryForMesh; // Mesh is deleted by AbstractCardiacTissue

        // We shouldn't ever have to save the old version
        assert(version >= 2);
//        {
//            bool use_matrix_based_assembly = true;
//            archive & use_matrix_based_assembly;
//        }

        archive & mWriteInfo;
        archive & mPrintOutput;
        archive & mNodesToOutput;
        //archive & mVoltageColumnId; // Created by InitialiseWriter, called from Solve
        //archive & mExtraVariablesId; // Created by InitialiseWriter, called from Solve
        //archive & mTimeColumnId; // Created by InitialiseWriter, called from Solve
        //archive & mNodeColumnId; // Created by InitialiseWriter, called from Solve
        //archive & mpWriter; // Created by InitialiseWriter, called from Solve
        archive & mpCardiacTissue;
        //archive & mpSolver; // Only exists during calls to the Solve method
        bool has_solution = (mSolution != NULL);
        archive & has_solution;
        if (has_solution)
        {
            std::string filename = ArchiveLocationInfo::GetArchiveDirectory() + "AbstractCardiacProblem_mSolution.vec";

            Hdf5DataWriter writer(*mpMesh->GetDistributedVectorFactory(), ArchiveLocationInfo::GetArchiveRelativePath(), "AbstractCardiacProblem_mSolution", false);
            writer.DefineFixedDimension(mpMesh->GetDistributedVectorFactory()->GetProblemSize());

            writer.DefineUnlimitedDimension("Time", "msec", 1);

            // Make sure the file does not take more disc space than really needed (related to #1200)
            writer.SetFixedChunkSize(1);

            int vm_col = writer.DefineVariable("Vm","mV");

            assert(HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering() == false );

            if (PROBLEM_DIM==1)
            {
                writer.EndDefineMode();
                writer.PutUnlimitedVariable(0.0);
                writer.PutVector(vm_col, mSolution);
            }

            if (PROBLEM_DIM==2)
            {
                int phie_col = writer.DefineVariable("Phie","mV");
                std::vector<int> variable_ids;
                variable_ids.push_back(vm_col);
                variable_ids.push_back(phie_col);
                writer.EndDefineMode();
                writer.PutUnlimitedVariable(0.0);
                writer.PutStripedVector(variable_ids, mSolution);
            }

            writer.Close();

        }
        archive & mCurrentTime;

        // Save boundary conditions
        SaveBoundaryConditions(archive, mpMesh, mpBoundaryConditionsContainer);
        SaveBoundaryConditions(archive, mpMesh, mpDefaultBoundaryConditionsContainer);
    }

    template<class Archive>
    void load(Archive & archive, const unsigned int version)
    {
        if (version >= 1)
        {
            unsigned element_dim;
            unsigned space_dim;
            unsigned problem_dim;
            archive & element_dim;
            archive & space_dim;
            archive & problem_dim;
            if ( (element_dim != ELEMENT_DIM) ||(space_dim != SPACE_DIM) ||(problem_dim != PROBLEM_DIM) )
            {
                /*If we carry on from this point then the mesh produced by unarchiving from the
                 * archive is templated as AbstractTetrahedralMesh<element_dim, space_dim>
                 * which doesn't match AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>*  mpMesh.
                 * Boost will through away the unarchived one, without deleting it properly and
                 * then set mpMesh=NULL.  We need to avoid this happening by bailing out.
                 */
                EXCEPTION("Failed to load from checkpoint because the dimensions of the archive do not match the object it's being read into.");
            }
        }
        archive & mMeshFilename;
        archive & mpMesh;
        assert(mpMesh != NULL); //If NULL then loading mesh has failed without an exception so Boost has given up on the mesh.  This would happen if a 2-dimensional mesh was successfully unarchived but mpMesh was expecting a 3-d mesh etc.
        //archive & mAllocatedMemoryForMesh; // Will always be true after a load

        if (version < 2)
        {
            bool use_matrix_based_assembly;
            archive & use_matrix_based_assembly;
        }

        archive & mWriteInfo;
        archive & mPrintOutput;
        archive & mNodesToOutput;
        //archive & mVoltageColumnId; // Created by InitialiseWriter, called from Solve
        //archive & mExtraVariablesId; // Created by InitialiseWriter, called from Solve
        //archive & mTimeColumnId; // Created by InitialiseWriter, called from Solve
        //archive & mNodeColumnId; // Created by InitialiseWriter, called from Solve
        //archive & mpWriter; // Created by InitialiseWriter, called from Solve
        archive & mpCardiacTissue;
        //archive & mpSolver; // Only exists during calls to the Solve method
        bool has_solution;
        archive & has_solution;
        if ((has_solution) && PROBLEM_DIM < 3)
        {
            std::string filename = ArchiveLocationInfo::GetArchiveDirectory() + "AbstractCardiacProblem_mSolution.vec";

            mSolution = mpMesh->GetDistributedVectorFactory()->CreateVec(PROBLEM_DIM);
            DistributedVector mSolution_distri = mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(mSolution);

            Vec vm = mpMesh->GetDistributedVectorFactory()->CreateVec();
            Vec phie = mpMesh->GetDistributedVectorFactory()->CreateVec();

            std::string archive_dir = ArchiveLocationInfo::GetArchiveRelativePath();
            Hdf5DataReader reader(archive_dir, "AbstractCardiacProblem_mSolution", !FileFinder::IsAbsolutePath(archive_dir));
            reader.GetVariableOverNodes(vm, "Vm", 0);

            if (PROBLEM_DIM==1)
            {
                //reader.Close(); // no need to call close explicitly, done in the destructor

                DistributedVector vm_distri = mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(vm);

                DistributedVector::Stripe mSolution_vm(mSolution_distri,0);

                for (DistributedVector::Iterator index = mSolution_distri.Begin();
                     index != mSolution_distri.End();
                     ++index)
                {
                    mSolution_vm[index] = vm_distri[index];
                }
            }

            if (PROBLEM_DIM==2)
            {
                reader.GetVariableOverNodes(phie, "Phie", 0);
                //reader.Close(); // no need to call close explicitly, done in the destructor

                DistributedVector vm_distri = mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(vm);
                DistributedVector phie_distri = mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(phie);

                DistributedVector::Stripe mSolution_vm(mSolution_distri,0);
                DistributedVector::Stripe mSolution_phie(mSolution_distri,1);

                for (DistributedVector::Iterator index = mSolution_distri.Begin();
                     index != mSolution_distri.End();
                     ++index)
                {
                    mSolution_vm[index] = vm_distri[index];
                    mSolution_phie[index] = phie_distri[index];
                }
            }

            mSolution_distri.Restore();

            VecDestroy(vm);
            VecDestroy(phie);

        }
        archive & mCurrentTime;

        // Load boundary conditions
        mpBoundaryConditionsContainer = LoadBoundaryConditions(archive, mpMesh);
        mpDefaultBoundaryConditionsContainer = LoadBoundaryConditions(archive, mpMesh);
    }

    BOOST_SERIALIZATION_SPLIT_MEMBER()

    
    template<class Archive>
    void SaveBoundaryConditions(Archive & archive,
                                AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
                                BccType pBcc) const
    {
        (*ProcessSpecificArchive<Archive>::Get()) & pBcc;
    }

    template<class Archive>
    BccType LoadBoundaryConditions(
            Archive & archive,
            AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
    {
        // Load pointer from archive
        BccType p_bcc;
        (*ProcessSpecificArchive<Archive>::Get()) & p_bcc;

        // Fill in the conditions, if we have a container and it's not already full
        if (p_bcc)
        {
            p_bcc->LoadFromArchive(*ProcessSpecificArchive<Archive>::Get(), pMesh);
        }

        return p_bcc;
    }

protected:
    std::string mMeshFilename;

    bool mAllocatedMemoryForMesh;
    bool mWriteInfo;
    bool mPrintOutput;

    std::vector<unsigned> mNodesToOutput;

    unsigned mVoltageColumnId;
    std::vector<unsigned> mExtraVariablesId;
    unsigned mTimeColumnId;
    unsigned mNodeColumnId;

    AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>* mpCardiacTissue;

    BccType mpBoundaryConditionsContainer;
    BccType mpDefaultBoundaryConditionsContainer;
    AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* mpSolver;
    AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* mpCellFactory;
    AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* mpMesh;

    Vec mSolution;

    double mCurrentTime;

    AbstractTimeAdaptivityController* mpTimeAdaptivityController;



    virtual AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>* CreateCardiacTissue() =0;

    virtual AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* CreateSolver() =0;

protected:
    template<unsigned DIM>
    friend class CardiacElectroMechanicsProblem;
    Hdf5DataWriter* mpWriter;

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

    AbstractCardiacProblem();

    virtual ~AbstractCardiacProblem();

    void Initialise();

    void SetNodesPerProcessorFilename(const std::string& rFilename);

    void SetBoundaryConditionsContainer(BccType pBcc);

    virtual void PreSolveChecks();

    virtual Vec CreateInitialCondition();

    void SetMesh(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);

    void PrintOutput(bool rPrintOutput);

    void SetWriteInfo(bool writeInfo = true);

    Vec GetSolution();

    DistributedVector GetSolutionDistributedVector();

    double GetCurrentTime();

    AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM> & rGetMesh();

    AbstractCardiacTissue<ELEMENT_DIM,SPACE_DIM>* GetTissue();

    void Solve();

    void CloseFilesAndPostProcess();

    virtual void WriteInfo(double time)=0;

    virtual void DefineWriterColumns(bool extending);

    void DefineExtraVariablesWriterColumns(bool extending);

    virtual void WriteOneStep(double time, Vec voltageVec) = 0;

    void WriteExtraVariablesOneStep();

    void InitialiseWriter();

    void SetOutputNodes(std::vector<unsigned> & rNodesToOutput);

    Hdf5DataReader GetDataReader();

    virtual void AtBeginningOfTimestep(double time)
    {}

    virtual void OnEndOfTimestep(double time)
    {}

    virtual void SetUpAdditionalStoppingTimes(std::vector<double>& rAdditionalStoppingTimes)
    {}



    // and no controller, in which case the default is used.

    void SetUseTimeAdaptivityController(bool useAdaptivity,
                                        AbstractTimeAdaptivityController* pController = NULL);

    template<class Archive>
    void LoadExtraArchive(Archive & archive, unsigned version);

    virtual bool GetHasBath();

    virtual void SetElectrodes();
};

TEMPLATED_CLASS_IS_ABSTRACT_3_UNSIGNED(AbstractCardiacProblem)


template<unsigned DIM>
class BidomainProblem;

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
template<class Archive>
void AbstractCardiacProblem<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::LoadExtraArchive(Archive & archive, unsigned version)
{
    // The vector factory must be loaded, but isn't needed for anything.
    DistributedVectorFactory* p_mesh_factory;
    archive >> p_mesh_factory;

    // The cardiac cells - load only the cells we actually own
    mpCardiacTissue->LoadCardiacCells(archive, version);

    {
        DistributedVectorFactory* p_pde_factory;
        archive >> p_pde_factory;
        assert(p_pde_factory == p_mesh_factory); // Paranoia...
    }
    // We no longer need this vector factory, since we already have our own.
    delete p_mesh_factory;

    // The boundary conditions
    BccType p_bcc;
    archive >> p_bcc;
    if (p_bcc)
    {
        if (!mpBoundaryConditionsContainer)
        {
            mpBoundaryConditionsContainer = p_bcc;
            mpBoundaryConditionsContainer->LoadFromArchive(archive, mpMesh);
        }
        else
        {
            // The BCs will only actually be different if using a distributed tetrahedral mesh
            DistributedTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* p_dist_mesh = dynamic_cast<DistributedTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>*>(mpMesh);
            if (p_dist_mesh)
            {
                mpBoundaryConditionsContainer->MergeFromArchive(archive, mpMesh);
            }
            else
            {
                // Load into the temporary container, which will get thrown away shortly
                p_bcc->LoadFromArchive(archive, mpMesh);
            }
        }
    }
    BccType p_default_bcc;
    archive >> p_default_bcc;
    if (p_default_bcc)
    {
        // This always holds, so we never need to load the default BCs
        assert(p_bcc == p_default_bcc);
    }

    // Are we a bidomain problem?
    if (PROBLEM_DIM == 2)
    {
        assert(ELEMENT_DIM == SPACE_DIM);
        BidomainProblem<ELEMENT_DIM>* p_problem = dynamic_cast<BidomainProblem<ELEMENT_DIM>*>(this);
        assert(p_problem);
        p_problem->LoadExtraArchiveForBidomain(archive, version);
    }
}

namespace boost {
namespace serialization {
template <unsigned ELEMENT_DIM, unsigned SPACE_DIM,  unsigned PROBLEM_DIM>
struct version<AbstractCardiacProblem<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM> >
{
    CHASTE_VERSION_CONTENT(2);
};
} // namespace serialization
} // namespace boost
#endif /*ABSTRACTCARDIACPROBLEM_HPP_*/