Chaste Release::3.1
LinearParabolicPdeSystemWithCoupledOdeSystemSolver.hpp
00001 /*
00002 
00003 Copyright (c) 2005-2012, University of Oxford.
00004 All rights reserved.
00005 
00006 University of Oxford means the Chancellor, Masters and Scholars of the
00007 University of Oxford, having an administrative office at Wellington
00008 Square, Oxford OX1 2JD, UK.
00009 
00010 This file is part of Chaste.
00011 
00012 Redistribution and use in source and binary forms, with or without
00013 modification, are permitted provided that the following conditions are met:
00014  * Redistributions of source code must retain the above copyright notice,
00015    this list of conditions and the following disclaimer.
00016  * Redistributions in binary form must reproduce the above copyright notice,
00017    this list of conditions and the following disclaimer in the documentation
00018    and/or other materials provided with the distribution.
00019  * Neither the name of the University of Oxford nor the names of its
00020    contributors may be used to endorse or promote products derived from this
00021    software without specific prior written permission.
00022 
00023 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00024 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00025 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00026 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
00027 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
00028 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
00029 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00030 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00031 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
00032 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00033 
00034 */
00035 #ifndef LINEARPARABOLICPDESYSTEMWITHCOUPLEDODESYSTEMSOLVER_HPP_
00036 #define LINEARPARABOLICPDESYSTEMWITHCOUPLEDODESYSTEMSOLVER_HPP_
00037 
00038 #include "AbstractAssemblerSolverHybrid.hpp"
00039 #include "AbstractDynamicLinearPdeSolver.hpp"
00040 #include "AbstractLinearParabolicPdeSystemForCoupledOdeSystem.hpp"
00041 #include "TetrahedralMesh.hpp"
00042 #include "BoundaryConditionsContainer.hpp"
00043 #include "AbstractOdeSystemForCoupledPdeSystem.hpp"
00044 #include "CvodeAdaptor.hpp"
00045 #include "BackwardEulerIvpOdeSolver.hpp"
00046 #include "Warnings.hpp"
00047 #include "VtkMeshWriter.hpp"
00048 
00049 #include <boost/shared_ptr.hpp>
00050 
00060 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM=ELEMENT_DIM, unsigned PROBLEM_DIM=1>
00061 class LinearParabolicPdeSystemWithCoupledOdeSystemSolver
00062     : public AbstractAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NORMAL>,
00063       public AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>
00064 {
00065 private:
00066 
00068     AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>* mpMesh;
00069 
00071     AbstractLinearParabolicPdeSystemForCoupledOdeSystem<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* mpPdeSystem;
00072 
00074     std::vector<AbstractOdeSystemForCoupledPdeSystem*> mOdeSystemsAtNodes;
00075 
00077     std::vector<double> mInterpolatedOdeStateVariables;
00078 
00080     boost::shared_ptr<AbstractIvpOdeSolver> mpOdeSolver;
00081 
00087     double mSamplingTimeStep;
00088 
00090     bool mOdeSystemsPresent;
00091 
00093     out_stream mpVtkMetaFile;
00094 
00098     void WriteVtkResultsToFile();
00099 
00110     c_matrix<double, PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1)> ComputeMatrixTerm(
00111         c_vector<double, ELEMENT_DIM+1>& rPhi,
00112         c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
00113         ChastePoint<SPACE_DIM>& rX,
00114         c_vector<double,PROBLEM_DIM>& rU,
00115         c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
00116         Element<ELEMENT_DIM, SPACE_DIM>* pElement);
00117 
00128     c_vector<double, PROBLEM_DIM*(ELEMENT_DIM+1)> ComputeVectorTerm(
00129         c_vector<double, ELEMENT_DIM+1>& rPhi,
00130         c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
00131         ChastePoint<SPACE_DIM>& rX,
00132         c_vector<double,PROBLEM_DIM>& rU,
00133         c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU,
00134         Element<ELEMENT_DIM, SPACE_DIM>* pElement);
00135 
00139     void ResetInterpolatedQuantities();
00140 
00148     void IncrementInterpolatedQuantities(double phiI, const Node<SPACE_DIM>* pNode);
00149 
00158     void InitialiseForSolve(Vec initialSolution=NULL);
00159 
00168     void SetupLinearSystem(Vec currentSolution, bool computeMatrix);
00169 
00170 public:
00171 
00181     LinearParabolicPdeSystemWithCoupledOdeSystemSolver(TetrahedralMesh<ELEMENT_DIM, SPACE_DIM>* pMesh,
00182                                                        AbstractLinearParabolicPdeSystemForCoupledOdeSystem<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* pPdeSystem,
00183                                                        BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* pBoundaryConditions,
00184                                                        std::vector<AbstractOdeSystemForCoupledPdeSystem*> odeSystemsAtNodes=std::vector<AbstractOdeSystemForCoupledPdeSystem*>(),
00185                                                        boost::shared_ptr<AbstractIvpOdeSolver> pOdeSolver=boost::shared_ptr<AbstractIvpOdeSolver>());
00186 
00190     ~LinearParabolicPdeSystemWithCoupledOdeSystemSolver();
00191 
00198     void PrepareForSetupLinearSystem(Vec currentPdeSolution);
00199 
00205     void SetOutputDirectory(std::string outputDirectory);
00206 
00212     void SetSamplingTimeStep(double samplingTimeStep);
00213 
00218     void SolveAndWriteResultsToFile();
00219 
00226     void WriteVtkResultsToFile(Vec solution, unsigned numTimeStepsElapsed);
00227 };
00228 
00230 // Implementation
00232 
00233 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00234 c_matrix<double, PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1)> LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ComputeMatrixTerm(
00235     c_vector<double, ELEMENT_DIM+1>& rPhi,
00236     c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
00237     ChastePoint<SPACE_DIM>& rX,
00238     c_vector<double,PROBLEM_DIM>& rU,
00239     c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
00240     Element<ELEMENT_DIM, SPACE_DIM>* pElement)
00241 {
00242     double timestep_inverse = PdeSimulationTime::GetPdeTimeStepInverse();
00243     c_matrix<double, PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1)> matrix_term = zero_matrix<double>(PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1));
00244 
00245     // Loop over PDEs and populate matrix_term
00246     for (unsigned pde_index=0; pde_index<PROBLEM_DIM; pde_index++)
00247     {
00248         double this_dudt_coefficient = mpPdeSystem->ComputeDuDtCoefficientFunction(rX, pde_index);
00249         c_matrix<double, SPACE_DIM, SPACE_DIM> this_pde_diffusion_term = mpPdeSystem->ComputeDiffusionTerm(rX, pde_index, pElement);
00250         c_matrix<double, 1*(ELEMENT_DIM+1), 1*(ELEMENT_DIM+1)> this_stiffness_matrix =
00251             prod(trans(rGradPhi), c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>(prod(this_pde_diffusion_term, rGradPhi)) )
00252                 + timestep_inverse * this_dudt_coefficient * outer_prod(rPhi, rPhi);
00253 
00254         for (unsigned i=0; i<ELEMENT_DIM+1; i++)
00255         {
00256             for (unsigned j=0; j<ELEMENT_DIM+1; j++)
00257             {
00258                 matrix_term(i*PROBLEM_DIM + pde_index, j*PROBLEM_DIM + pde_index) = this_stiffness_matrix(i,j);
00259             }
00260         }
00261     }
00262     return matrix_term;
00263 }
00264 
00265 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00266 c_vector<double, PROBLEM_DIM*(ELEMENT_DIM+1)> LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ComputeVectorTerm(
00267     c_vector<double, ELEMENT_DIM+1>& rPhi,
00268     c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
00269     ChastePoint<SPACE_DIM>& rX,
00270     c_vector<double,PROBLEM_DIM>& rU,
00271     c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU,
00272     Element<ELEMENT_DIM, SPACE_DIM>* pElement)
00273 {
00274     double timestep_inverse = PdeSimulationTime::GetPdeTimeStepInverse();
00275     c_vector<double, PROBLEM_DIM*(ELEMENT_DIM+1)> vector_term = zero_vector<double>(PROBLEM_DIM*(ELEMENT_DIM+1));
00276 
00277     // Loop over PDEs and populate vector_term
00278     for (unsigned pde_index=0; pde_index<PROBLEM_DIM; pde_index++)
00279     {
00280         double this_dudt_coefficient = mpPdeSystem->ComputeDuDtCoefficientFunction(rX, pde_index);
00281         double this_source_term = mpPdeSystem->ComputeSourceTerm(rX, rU, mInterpolatedOdeStateVariables, pde_index);
00282         c_vector<double, ELEMENT_DIM+1> this_vector_term = (this_source_term + timestep_inverse*this_dudt_coefficient*rU(pde_index))* rPhi;
00283 
00284         for (unsigned i=0; i<ELEMENT_DIM+1; i++)
00285         {
00286             vector_term(i*PROBLEM_DIM + pde_index) = this_vector_term(i);
00287         }
00288     }
00289 
00290     return vector_term;
00291 }
00292 
00293 
00294 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00295 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ResetInterpolatedQuantities()
00296 {
00297     mInterpolatedOdeStateVariables.clear();
00298 
00299     if (mOdeSystemsPresent)
00300     {
00301         unsigned num_state_variables = mOdeSystemsAtNodes[0]->GetNumberOfStateVariables();
00302         mInterpolatedOdeStateVariables.resize(num_state_variables, 0.0);
00303     }
00304 }
00305 
00306 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00307 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::IncrementInterpolatedQuantities(double phiI, const Node<SPACE_DIM>* pNode)
00308 {
00309     if (mOdeSystemsPresent)
00310     {
00311         unsigned num_state_variables = mOdeSystemsAtNodes[0]->GetNumberOfStateVariables();
00312 
00313         for (unsigned i=0; i<num_state_variables; i++)
00314         {
00315             mInterpolatedOdeStateVariables[i] += phiI * mOdeSystemsAtNodes[pNode->GetIndex()]->rGetStateVariables()[i];
00316         }
00317     }
00318 }
00319 
00320 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00321 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::InitialiseForSolve(Vec initialSolution)
00322 {
00323     if (this->mpLinearSystem == NULL)
00324     {
00325         unsigned preallocation = mpMesh->CalculateMaximumContainingElementsPerProcess() + ELEMENT_DIM;
00326         if (ELEMENT_DIM > 1)
00327         {
00328             // Highest connectivity is closed
00329             preallocation--;
00330         }
00331         preallocation *= PROBLEM_DIM;
00332 
00333         /*
00334          * Use the current solution (ie the initial solution) as the
00335          * template in the alternative constructor of LinearSystem.
00336          * This is to avoid problems with VecScatter.
00337          */
00338         this->mpLinearSystem = new LinearSystem(initialSolution, preallocation);
00339     }
00340 
00341     assert(this->mpLinearSystem);
00342     this->mpLinearSystem->SetMatrixIsSymmetric(true);
00343     this->mpLinearSystem->SetKspType("cg");
00344 }
00345 
00346 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00347 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetupLinearSystem(Vec currentSolution, bool computeMatrix)
00348 {
00349     SetupGivenLinearSystem(currentSolution, computeMatrix, this->mpLinearSystem);
00350 }
00351 
00352 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00353 LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::LinearParabolicPdeSystemWithCoupledOdeSystemSolver(
00354         TetrahedralMesh<ELEMENT_DIM, SPACE_DIM>* pMesh,
00355         AbstractLinearParabolicPdeSystemForCoupledOdeSystem<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* pPdeSystem,
00356         BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>* pBoundaryConditions,
00357         std::vector<AbstractOdeSystemForCoupledPdeSystem*> odeSystemsAtNodes,
00358         boost::shared_ptr<AbstractIvpOdeSolver> pOdeSolver)
00359     : AbstractAssemblerSolverHybrid<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, NORMAL>(pMesh, pBoundaryConditions),
00360       AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>(pMesh),
00361       mpMesh(pMesh),
00362       mpPdeSystem(pPdeSystem),
00363       mOdeSystemsAtNodes(odeSystemsAtNodes),
00364       mpOdeSolver(pOdeSolver),
00365       mSamplingTimeStep(DOUBLE_UNSET),
00366       mOdeSystemsPresent(false)
00367 {
00368     this->mpBoundaryConditions = pBoundaryConditions;
00369 
00370     /*
00371      * If any ODE systems are passed in to the constructor, then we aren't just
00372      * solving a coupled PDE system, in which case the number of ODE system objects
00373      * must match the number of nodes in the finite element mesh.
00374      */
00375     if (!mOdeSystemsAtNodes.empty())
00376     {
00377         mOdeSystemsPresent = true;
00378         assert(mOdeSystemsAtNodes.size() == mpMesh->GetNumNodes());
00379 
00380         /*
00381          * In this case, if an ODE solver is not explicitly passed into the
00382          * constructor, then we create a default solver.
00383          */
00384         if (!mpOdeSolver)
00385         {
00386 #ifdef CHASTE_CVODE
00387             mpOdeSolver.reset(new CvodeAdaptor);
00388 #else
00389             mpOdeSolver.reset(new BackwardEulerIvpOdeSolver(mOdeSystemsAtNodes[0]->GetNumberOfStateVariables()));
00390 #endif //CHASTE_CVODE
00391         }
00392     }
00393 }
00394 
00395 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00396 LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::~LinearParabolicPdeSystemWithCoupledOdeSystemSolver()
00397 {
00398     if (mOdeSystemsPresent)
00399     {
00400         for (unsigned i=0; i<mOdeSystemsAtNodes.size(); i++)
00401         {
00402             delete mOdeSystemsAtNodes[i];
00403         }
00404     }
00405 }
00406 
00407 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00408 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::PrepareForSetupLinearSystem(Vec currentPdeSolution)
00409 {
00410     if (mOdeSystemsPresent)
00411     {
00412         double time = PdeSimulationTime::GetTime();
00413         double dt = PdeSimulationTime::GetPdeTimeStep();
00414 
00415         ReplicatableVector soln_repl(currentPdeSolution);
00416         std::vector<double> current_soln_this_node(PROBLEM_DIM);
00417 
00418         // Loop over nodes
00419         for (unsigned node_index=0; node_index<mpMesh->GetNumNodes(); node_index++)
00420         {
00421             // Store the current solution to the PDE system at this node
00422             for (unsigned pde_index=0; pde_index<PROBLEM_DIM; pde_index++)
00423             {
00424                 double current_soln_this_pde_this_node = soln_repl[PROBLEM_DIM*node_index + pde_index];
00425                 current_soln_this_node[pde_index] = current_soln_this_pde_this_node;
00426             }
00427 
00428             // Pass it into the ODE system at this node
00429             mOdeSystemsAtNodes[node_index]->SetPdeSolution(current_soln_this_node);
00430 
00431             // Solve ODE system at this node
00432             mpOdeSolver->SolveAndUpdateStateVariable(mOdeSystemsAtNodes[node_index], time, time+dt, dt);
00433         }
00434     }
00435 }
00436 
00437 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00438 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetOutputDirectory(std::string outputDirectory)
00439 {
00440     this->mOutputDirectory = outputDirectory;
00441 }
00442 
00443 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00444 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetSamplingTimeStep(double samplingTimeStep)
00445 {
00446     assert(samplingTimeStep >= this->mIdealTimeStep);
00447     mSamplingTimeStep = samplingTimeStep;
00448 }
00449 
00450 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00451 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SolveAndWriteResultsToFile()
00452 {
00453     // A number of methods must have been called prior to this method
00454     if (this->mOutputDirectory == "")
00455     {
00456         EXCEPTION("SetOutputDirectory() must be called prior to SolveAndWriteResultsToFile()");
00457     }
00458     if (this->mTimesSet == false)
00459     {
00460         EXCEPTION("SetTimes() must be called prior to SolveAndWriteResultsToFile()");
00461     }
00462     if (this->mIdealTimeStep <= 0.0)
00463     {
00464         EXCEPTION("SetTimeStep() must be called prior to SolveAndWriteResultsToFile()");
00465     }
00466     if (mSamplingTimeStep == DOUBLE_UNSET)
00467     {
00468         EXCEPTION("SetSamplingTimeStep() must be called prior to SolveAndWriteResultsToFile()");
00469     }
00470     if (!this->mInitialCondition)
00471     {
00472         EXCEPTION("SetInitialCondition() must be called prior to SolveAndWriteResultsToFile()");
00473     }
00474 
00475 #ifdef CHASTE_VTK
00476     // Create a .pvd output file
00477     OutputFileHandler output_file_handler(this->mOutputDirectory, false);
00478     mpVtkMetaFile = output_file_handler.OpenOutputFile("results.pvd");
00479     *mpVtkMetaFile << "<?xml version=\"1.0\"?>\n";
00480     *mpVtkMetaFile << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n";
00481     *mpVtkMetaFile << "    <Collection>\n";
00482 
00483     // Write initial condition to VTK
00484     Vec initial_condition = this->mInitialCondition;
00485     WriteVtkResultsToFile(initial_condition, 0);
00486 
00487     // The helper class TimeStepper deals with issues such as small final timesteps so we don't have to
00488     TimeStepper stepper(this->mTstart, this->mTend, mSamplingTimeStep);
00489 
00490     // Main time loop
00491     while (!stepper.IsTimeAtEnd())
00492     {
00493         // Reset start and end times
00494         this->SetTimes(stepper.GetTime(), stepper.GetNextTime());
00495 
00496         // Solve the system up to the new end time
00497         Vec soln = this->Solve();
00498 
00499         // Reset the initial condition for the next timestep
00500         if (this->mInitialCondition != initial_condition)
00501         {
00502             PetscTools::Destroy(this->mInitialCondition);
00503         }
00504         this->mInitialCondition = soln;
00505 
00506         // Move forward in time
00507         stepper.AdvanceOneTimeStep();
00508 
00509         // Write solution to VTK
00510         WriteVtkResultsToFile(soln, stepper.GetTotalTimeStepsTaken());
00511     }
00512 
00513     // Restore saved initial condition to avoid user confusion!
00514     if (this->mInitialCondition != initial_condition)
00515     {
00516         PetscTools::Destroy(this->mInitialCondition);
00517     }
00518     this->mInitialCondition = initial_condition;
00519 
00520     // Close .pvd output file
00521     *mpVtkMetaFile << "    </Collection>\n";
00522     *mpVtkMetaFile << "</VTKFile>\n";
00523     mpVtkMetaFile->close();
00524 #else //CHASTE_VTK
00525 #define COVERAGE_IGNORE // We only test this in weekly builds
00526     WARNING("VTK is not installed and is required for this functionality");
00527 #undef COVERAGE_IGNORE
00528 #endif //CHASTE_VTK
00529 }
00530 
00531 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00532 void LinearParabolicPdeSystemWithCoupledOdeSystemSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::WriteVtkResultsToFile(Vec solution, unsigned numTimeStepsElapsed)
00533 {
00534 #ifdef CHASTE_VTK
00535 
00536     // Create a new VTK file for this time step
00537     std::stringstream time;
00538     time << numTimeStepsElapsed;
00539     VtkMeshWriter<ELEMENT_DIM, SPACE_DIM> mesh_writer(this->mOutputDirectory, "results_"+time.str(), false);
00540 
00541     /*
00542      * We first loop over PDEs. For each PDE we store the solution
00543      * at each node in a vector, then pass this vector to the mesh
00544      * writer.
00545      */
00546     ReplicatableVector solution_repl(solution);
00547     unsigned num_nodes = mpMesh->GetNumNodes();
00548     for (unsigned pde_index=0; pde_index<PROBLEM_DIM; pde_index++)
00549     {
00550         // Store the solution of this PDE at each node
00551         std::vector<double> pde_index_data;
00552         pde_index_data.resize(num_nodes, 0.0);
00553         for (unsigned node_index=0; node_index<num_nodes; node_index++)
00554         {
00555             pde_index_data[node_index] = solution_repl[PROBLEM_DIM*node_index + pde_index];
00556         }
00557 
00558         // Add this data to the mesh writer
00559         std::stringstream data_name;
00560         data_name << "PDE variable " << pde_index;
00561         mesh_writer.AddPointData(data_name.str(), pde_index_data);
00562     }
00563 
00564     if (mOdeSystemsPresent)
00565     {
00566         /*
00567          * We cannot loop over ODEs like PDEs, since the solutions are not
00568          * stored in one place. Therefore we build up a large 'vector of
00569          * vectors', then pass each component of this vector to the mesh
00570          * writer.
00571          */
00572         std::vector<std::vector<double> > ode_data;
00573         unsigned num_odes = mOdeSystemsAtNodes[0]->rGetStateVariables().size();
00574         ode_data.resize(num_odes);
00575         for (unsigned ode_index=0; ode_index<num_odes; ode_index++)
00576         {
00577             ode_data[ode_index].resize(num_nodes, 0.0);
00578         }
00579 
00580         for (unsigned node_index=0; node_index<num_nodes; node_index++)
00581         {
00582             std::vector<double> all_odes_this_node = mOdeSystemsAtNodes[node_index]->rGetStateVariables();
00583             for (unsigned i=0; i<num_odes; i++)
00584             {
00585                 ode_data[i][node_index] = all_odes_this_node[i];
00586             }
00587         }
00588 
00589         for (unsigned ode_index=0; ode_index<num_odes; ode_index++)
00590         {
00591             std::vector<double> ode_index_data = ode_data[ode_index];
00592 
00593             // Add this data to the mesh writer
00594             std::stringstream data_name;
00595             data_name << "ODE variable " << ode_index;
00596             mesh_writer.AddPointData(data_name.str(), ode_index_data);
00597         }
00598     }
00599 
00600     mesh_writer.WriteFilesUsingMesh(*mpMesh);
00601     *mpVtkMetaFile << "        <DataSet timestep=\"";
00602     *mpVtkMetaFile << numTimeStepsElapsed;
00603     *mpVtkMetaFile << "\" group=\"\" part=\"0\" file=\"results_";
00604     *mpVtkMetaFile << numTimeStepsElapsed;
00605     *mpVtkMetaFile << ".vtu\"/>\n";
00606 #endif // CHASTE_VTK
00607 }
00608 
00609 #endif /*LINEARPARABOLICPDESYSTEMWITHCOUPLEDODESYSTEMSOLVER_HPP_*/