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

#include <string>
#include <vector>
#include <cassert>
#include <climits> // Work around a boost bug - see #1024.
#include <boost/serialization/vector.hpp>
#include <boost/serialization/string.hpp>
#include "ChasteSerialization.hpp"
#include <boost/serialization/split_member.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/serialization/shared_ptr.hpp>
#include "ClassIsAbstract.hpp"
#include "AbstractCardiacCell.hpp"
#include "AbstractCardiacCellFactory.hpp"
#include "AbstractCardiacPde.hpp"
#include "AbstractDynamicAssemblerMixin.hpp"
#include "AbstractTetrahedralMesh.hpp"
#include "BoundaryConditionsContainer.hpp"
#include "DistributedVectorFactory.hpp"
#include "Hdf5DataReader.hpp"
#include "Hdf5DataWriter.hpp"

/*
 * Archiving extravaganza:
 *
 * We archive mesh and pde 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, 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 implicitely 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 "MonodomainPde.hpp"
#include "BidomainPde.hpp"


template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
class AbstractCardiacProblem
{
    friend class TestBidomainWithBathAssembler;
    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
    {
        archive & mMeshFilename;
        archive & mpMesh;
        //archive & mAllocatedMemoryForMesh; // Mesh is deleted by AbstractCardiacPde
        archive & mUseMatrixBasedRhsAssembly;
        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 & mpCardiacPde;
        //archive & mpAssembler; // 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");
            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");
                writer.EndDefineMode();
                writer.PutUnlimitedVariable(0.0);
                writer.PutStripedVector(vm_col, phie_col, 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)
    {
        archive & mMeshFilename;
        archive & mpMesh;
        //archive & mAllocatedMemoryForMesh; // Will always be true after a load
        archive & mUseMatrixBasedRhsAssembly;
        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 & mpCardiacPde;
        //archive & mpAssembler; // Only exists during calls to the Solve method
        bool has_solution;
        archive & has_solution;
        if (has_solution)
        {
            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();

            Hdf5DataReader reader(ArchiveLocationInfo::GetArchiveRelativePath(), "AbstractCardiacProblem_mSolution", ArchiveLocationInfo::GetIsDirRelativeToChasteTestOutput());
            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 mUseMatrixBasedRhsAssembly;
    bool mAllocatedMemoryForMesh;
    bool mWriteInfo;
    bool mPrintOutput;

    std::vector<unsigned> mNodesToOutput;

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

    AbstractCardiacPde<ELEMENT_DIM,SPACE_DIM>* mpCardiacPde;

    BccType mpBoundaryConditionsContainer;
    BccType mpDefaultBoundaryConditionsContainer;
    AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* mpAssembler;
    AbstractCardiacCellFactory<ELEMENT_DIM,SPACE_DIM>* mpCellFactory;
    AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* mpMesh;

    Vec mSolution;

    double mCurrentTime;

    virtual AbstractCardiacPde<ELEMENT_DIM,SPACE_DIM>* CreateCardiacPde() =0;

    virtual AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* CreateAssembler() =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();

    AbstractCardiacPde<ELEMENT_DIM,SPACE_DIM>* GetPde();

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

    void UseMatrixBasedRhsAssembly(bool usematrix=true);

    virtual void AtBeginningOfTimestep(double time)
    {}

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

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

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
    std::vector<AbstractCardiacCell*> cells;
    // Load only the cells we actually own
    AbstractCardiacPde<ELEMENT_DIM,SPACE_DIM>::LoadCardiacCells(archive, version, cells, this->mpMesh);
    mpCardiacPde->MergeCells(cells);

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


#endif /*ABSTRACTCARDIACPROBLEM_HPP_*/