Chaste Release::3.1
|
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_*/