AbstractLinearAssembler.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 _ABSTRACTLINEARASSEMBLER_HPP_
#define _ABSTRACTLINEARASSEMBLER_HPP_


#include <vector>
#include <iostream>
#include <petscvec.h>

#include "AbstractStaticAssembler.hpp"


template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
class AbstractLinearAssembler : public AbstractStaticAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>
{
private:

    bool mMatrixIsConstant;

protected:

    void SetMatrixIsConst(bool matrixIsConstant=true);

    void ApplyDirichletConditions(Vec unusedVector, bool applyToMatrix);

    virtual void InitialiseForSolve(Vec initialSolution);

    bool ProblemIsNonlinear();

    virtual Vec StaticSolve(Vec currentSolutionOrGuess=NULL,
                            double currentTime=0.0,
                            bool assembleMatrix=true);

public:

    AbstractLinearAssembler(unsigned numQuadPoints=2);

    ~AbstractLinearAssembler();

    virtual Vec Solve(Vec currentSolutionOrGuess=NULL, double currentTime=0.0);

    /*
    void DebugWithSolution(Vec sol)
    {
        std::cout<<"\n\nWS: This is the matrix:\n";
        mpLinearSystem->DisplayMatrix();
        std::cout<<"\n\nWS: This is the righthand side:\n";
        mpLinearSystem->DisplayRhs();
        std::cout<<"\n\nWS: This is the solution:\n";
        VecView(sol, PETSC_VIEWER_STDOUT_WORLD);
    }

    void Debug()
    {
        std::cout<<"\n\nThis is the matrix:\n";
        mpLinearSystem->DisplayMatrix();
        std::cout<<"\n\nThis is the righthand side:\n";
        mpLinearSystem->DisplayRhs();
    }
    */
};


// Implementation


template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
void AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>::SetMatrixIsConst(bool matrixIsConstant)
{
     mMatrixIsConstant = matrixIsConstant;
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
void AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>::ApplyDirichletConditions(Vec unusedVector, bool applyToMatrix)
{
    this->mpBoundaryConditions->ApplyDirichletToLinearProblem(*(this->mpLinearSystem), applyToMatrix);
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
void AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>::InitialiseForSolve(Vec initialSolution)
{
    if (this->mpLinearSystem == NULL)
    {
        HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
        if (initialSolution == NULL)
        {
            // Static problem, create linear system
            // The following ensures all the unknowns for a particular node
            // are on the same processor
            Vec template_vec = this->mpMesh->GetDistributedVectorFactory()->CreateVec(PROBLEM_DIM);

            this->mpLinearSystem = new LinearSystem(template_vec);

            VecDestroy(template_vec);
        }
        else
        {
            // Use the currrent solution (ie the initial solution)
            // as the template in the alternative constructor of
            // LinearSystem. This is to avoid problems with VecScatter.
            this->mpLinearSystem = new LinearSystem(initialSolution);
        }

        this->mpLinearSystem->SetMatrixIsConstant(mMatrixIsConstant);
        HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
    }
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
bool AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>::ProblemIsNonlinear()
{
    return false;
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
Vec AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>::StaticSolve(Vec currentSolutionOrGuess,
                        double currentTime,
                        bool assembleMatrix)
{
    this->AssembleSystem(true, assembleMatrix, currentSolutionOrGuess, currentTime);
    return this->mpLinearSystem->Solve(currentSolutionOrGuess);
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>::AbstractLinearAssembler(unsigned numQuadPoints)
    : AbstractStaticAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM, NON_HEART, CONCRETE>(numQuadPoints),
      mMatrixIsConstant(true)
{
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>::~AbstractLinearAssembler()
{
}

template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, bool NON_HEART, class CONCRETE>
Vec AbstractLinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NON_HEART, CONCRETE>::Solve(Vec currentSolutionOrGuess, double currentTime)
{
    // have to comment this out because the flagged mesh assembler uses a different
    // bcc and so this is null.
    //assert(this->mpBoundaryConditions!=NULL);

    this->PrepareForSolve();
    this->InitialiseForSolve(currentSolutionOrGuess);
    return this->StaticSolve(currentSolutionOrGuess, currentTime);
}

#endif //_ABSTRACTLINEARASSEMBLER_HPP_