AbstractDynamicAssemblerMixin.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 #ifndef _ABSTRACTDYNAMICASSEMBLERMIXIN_HPP_
00029 #define _ABSTRACTDYNAMICASSEMBLERMIXIN_HPP_
00030 
00031 #include <vector>
00032 #include <petscvec.h>
00033 
00034 #include "AbstractAssembler.hpp"
00035 #include "TimeStepper.hpp"
00036 #include "PdeSimulationTime.hpp"
00037 
00045 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00046 class AbstractDynamicAssemblerMixin : virtual public AbstractAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>
00047 {
00048 protected:
00049 
00051     double mTstart;
00052 
00054     double mTend;
00055 
00057     double mDt;
00058 
00060     double mDtInverse;
00061 
00063     bool mTimesSet;
00064 
00066     Vec mInitialCondition;
00067 
00069     bool mMatrixIsAssembled;
00070 
00073     bool mMatrixIsConstant;
00074 
00076     bool mUseMatrixBasedRhsAssembly;
00077 
00079     Mat* mpMatrixForMatrixBasedRhsAssembly;
00080 
00082     Vec mVectorForMatrixBasedRhsAssembly;
00083 
00094     void DoMatrixBasedRhsAssembly(Vec currentSolution, double time);
00095 
00096 public:
00097 
00102     AbstractDynamicAssemblerMixin();
00103 
00111     void SetTimes(double tStart, double tEnd, double dt);
00112 
00118     void SetInitialCondition(Vec initialCondition);
00119 
00125     void SetMatrixIsConstant(bool matrixIsConstant=true);
00126 
00130     void SetMatrixIsNotAssembled();
00131 
00146     Vec Solve(Vec currentSolutionOrGuess=NULL, double currentTime=0.0);
00147 
00154     virtual void ConstructVectorForMatrixBasedRhsAssembly(Vec currentSolution);
00155 
00156 };
00157 
00158 
00160 // Implementation
00162 
00163 
00164 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00165 void AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::DoMatrixBasedRhsAssembly(Vec currentSolution, double time)
00166 {
00167     assert(mpMatrixForMatrixBasedRhsAssembly!=NULL);
00168 
00169     // as bypassing AssembleSystem, need to make sure we call
00170     // Prepare and Finalize
00171     this->PrepareForAssembleSystem(currentSolution, time);
00172 
00173     HeartEventHandler::BeginEvent(HeartEventHandler::ASSEMBLE_RHS);
00174 
00175     (*(this->GetLinearSystem()))->ZeroRhsVector();
00176 
00177     // construct z
00178     ConstructVectorForMatrixBasedRhsAssembly(currentSolution);
00179 
00180     // b = Bz
00181     MatMult(*mpMatrixForMatrixBasedRhsAssembly, mVectorForMatrixBasedRhsAssembly, (*(this->GetLinearSystem()))->rGetRhsVector());
00182 
00183     // apply boundary conditions
00184     this->ApplyNeummanBoundaryConditions();
00185     (*(this->GetLinearSystem()))->AssembleRhsVector();
00186 
00187     this->ApplyDirichletConditions(currentSolution, false);
00188 
00189     // as bypassing AssembleSystem, need to make sure we call
00190     // Prepare and Finalise
00191     this->FinaliseAssembleSystem(currentSolution, time);
00192     (*(this->GetLinearSystem()))->AssembleRhsVector();
00193 
00194     HeartEventHandler::EndEvent(HeartEventHandler::ASSEMBLE_RHS);
00195 }
00196 
00197 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00198 AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::AbstractDynamicAssemblerMixin()
00199     : mTimesSet(false),
00200       mInitialCondition(NULL),
00201       mMatrixIsAssembled(false),
00202       mMatrixIsConstant(false),
00203       mUseMatrixBasedRhsAssembly(false),
00204       mpMatrixForMatrixBasedRhsAssembly(NULL)
00205 {
00206 }
00207 
00208 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00209 void AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetTimes(double tStart, double tEnd, double dt)
00210 {
00211     mTstart = tStart;
00212     mTend   = tEnd;
00213 
00214     if (mTstart >= mTend)
00215     {
00216         EXCEPTION("Starting time has to less than ending time");
00217     }
00218 
00219     mDt     = dt;
00220 
00221     if (mDt <= 0)
00222     {
00223         EXCEPTION("Time step has to be greater than zero");
00224     }
00225 
00226     mDtInverse = 1/dt;
00227     mTimesSet = true;
00228 }
00229 
00230 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00231 void AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetInitialCondition(Vec initialCondition)
00232 {
00233     assert(initialCondition!=NULL);
00234     mInitialCondition = initialCondition;
00235 }
00236 
00237 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00238 void AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetMatrixIsConstant(bool matrixIsConstant)
00239 {
00240     mMatrixIsConstant = matrixIsConstant;
00241     this->SetMatrixIsConst(mMatrixIsConstant);
00242 }
00243 
00244 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00245 void AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetMatrixIsNotAssembled()
00246 {
00247     mMatrixIsAssembled = false;
00248 }
00249 
00250 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00251 Vec AbstractDynamicAssemblerMixin<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::Solve(Vec currentSolutionOrGuess, double currentTime)
00252 {
00253     assert(mTimesSet);
00254     assert(mInitialCondition != NULL);
00255 
00256     this->PrepareForSolve();
00257     this->InitialiseForSolve(mInitialCondition);
00258 
00259     TimeStepper stepper(mTstart, mTend, mDt, mMatrixIsConstant);
00260 
00261     Vec current_solution = mInitialCondition;
00262     Vec next_solution;
00263 
00264     while ( !stepper.IsTimeAtEnd() )
00265     {
00266         mDt = stepper.GetNextTimeStep();
00267         mDtInverse = 1.0/mDt;
00268 
00269         PdeSimulationTime::SetTime(stepper.GetTime());
00270 
00271         // NOTE: even if mUseMatrixBasedRhsAssembly==true,
00272         // the RHS is assembled without using matrix-based assembly
00273         // in the first timestep (when the LHS matrix is set up) - no
00274         // easy way around this
00275         if (!mUseMatrixBasedRhsAssembly || !mMatrixIsAssembled)
00276         {
00277             // matrix is constant case: only assembler matrix the first time
00278             // matrix is not constant case: always assemble
00279             bool assemble_matrix = (!mMatrixIsConstant || !mMatrixIsAssembled);
00280 
00281             next_solution = this->StaticSolve(current_solution, stepper.GetTime(), assemble_matrix);
00282         }
00283         else
00284         {
00285             DoMatrixBasedRhsAssembly(current_solution, stepper.GetTime());
00286             next_solution = (*(this->GetLinearSystem()))->Solve(current_solution);
00287         }
00288 
00289         if (mMatrixIsConstant)
00290         {
00291             mMatrixIsAssembled = true;
00292         }
00293 
00294         stepper.AdvanceOneTimeStep();
00295 
00296         // Avoid memory leaks
00297         if (current_solution != mInitialCondition)
00298         {
00299             HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00300             VecDestroy(current_solution);
00301             HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00302         }
00303         current_solution = next_solution;
00304     }
00305     return current_solution;
00306 }
00307 
00308 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00309 void AbstractDynamicAssemblerMixin<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM>::ConstructVectorForMatrixBasedRhsAssembly(Vec currentSolution)
00310 {
00311     #define COVERAGE_IGNORE
00312     EXCEPTION("mUseMatrixBasedRhsAssembly=true but ConstructVectorForMatrixBasedRhsAssembly() has not been overloaded");
00313     #undef COVERAGE_IGNORE
00314 }
00315 
00316 #endif //_ABSTRACTDYNAMICASSEMBLERMIXIN_HPP_

Generated by  doxygen 1.6.2