AbstractCardiacProblem.hpp

/*

Copyright (c) 2005-2014, University of Oxford.
All rights reserved.

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.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.
 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.
 * Neither the name of the University of Oxford nor the names of its
   contributors may be used to endorse or promote products derived from this
   software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/


#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"
#include "Warnings.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 : private 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)
        {
            Hdf5DataWriter writer(*mpMesh->GetDistributedVectorFactory(), ArchiveLocationInfo::GetArchiveRelativePath(), "AbstractCardiacProblem_mSolution", false);
            writer.DefineFixedDimension(mpMesh->GetDistributedVectorFactory()->GetProblemSize());
            writer.DefineUnlimitedDimension("Time", "msec", 1);

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

            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)
        {
            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();

            PetscTools::Destroy(vm);
            PetscTools::Destroy(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;

    virtual void CreateMeshFromHeartConfig();

    template<unsigned DIM, unsigned ELEC_PROB_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();

    bool 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;

    // How many processes were used by the saving simulation?
    DistributedVectorFactory* p_original_factory = p_mesh_factory->GetOriginalFactory();
    unsigned orig_num_procs = 1;
    if (p_original_factory)
    {
        orig_num_procs = p_original_factory->GetNumProcs();
    }

    // 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
        {
            // If the mesh which was archived was a TetrahedralMesh then we have all the boundary conditions
            // in every process-specific archive.  We no longer test for this.
#define COVERAGE_IGNORE
            if(!dynamic_cast<DistributedTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>*>(mpMesh) && orig_num_procs > 1)
            {
                // The correct way to do this should be:
                // p_bcc->LoadFromArchive(archive, mpMesh);
                WARNING("Loading from a parallel archive which used a non-distributed mesh.  This scenario should work but is not fully tested.");
            }
#undef COVERAGE_IGNORE
            mpBoundaryConditionsContainer->MergeFromArchive(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?
    BidomainProblem<ELEMENT_DIM>* p_bidomain_problem = dynamic_cast<BidomainProblem<ELEMENT_DIM>*>(this);
    if (p_bidomain_problem)
    {
        assert(ELEMENT_DIM == SPACE_DIM);
        p_bidomain_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_*/