CellBasedSimulationWithPdes.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 
00029 
00030 #include "CellBasedSimulationWithPdes.hpp"
00031 #include "MeshBasedCellPopulationWithGhostNodes.hpp"
00032 #include "SimpleDataWriter.hpp"
00033 #include "BoundaryConditionsContainer.hpp"
00034 #include "ConstBoundaryCondition.hpp"
00035 #include "SimpleLinearEllipticSolver.hpp"
00036 #include "CellBasedSimulationWithPdesAssembler.hpp"
00037 #include "CellwiseData.hpp"
00038 #include "AbstractTwoBodyInteractionForce.hpp"
00039 #include "TrianglesMeshReader.hpp"
00040 
00041 template<unsigned DIM>
00042 CellBasedSimulationWithPdes<DIM>::CellBasedSimulationWithPdes(AbstractCellPopulation<DIM>& rCellPopulation,
00043                                                         std::vector<PdeAndBoundaryConditions<DIM>*> pdeAndBcCollection,
00044                                                         bool deleteCellPopulationAndForceCollection,
00045                                                         bool initialiseCells)
00046     : CellBasedSimulation<DIM>(rCellPopulation,
00047                                deleteCellPopulationAndForceCollection,
00048                                initialiseCells),
00049       mPdeAndBcCollection(pdeAndBcCollection),
00050       mWriteAverageRadialPdeSolution(false),
00051       mWriteDailyAverageRadialPdeSolution(false),
00052       mNumRadialIntervals(0), // 'unset' value
00053       mpCoarsePdeMesh(NULL)
00054 {
00055     // We must be using a mesh-based cell population
00056     assert(dynamic_cast<MeshBasedCellPopulation<DIM>*>(&(this->mrCellPopulation)) != NULL);
00057 
00058     // We must not have any ghost nodes
00059     assert(dynamic_cast<MeshBasedCellPopulationWithGhostNodes<DIM>*>(&(this->mrCellPopulation)) == NULL);
00060 }
00061 
00062 template<unsigned DIM>
00063 CellBasedSimulationWithPdes<DIM>::~CellBasedSimulationWithPdes()
00064 {
00065     if (mpCoarsePdeMesh)
00066     {
00067         delete mpCoarsePdeMesh;
00068     }
00069 }
00070 
00071 template<unsigned DIM>
00072 void CellBasedSimulationWithPdes<DIM>::SetPdeAndBcCollection(std::vector<PdeAndBoundaryConditions<DIM>*> pdeAndBcCollection)
00073 {
00074     mPdeAndBcCollection = pdeAndBcCollection;
00075 }
00076 
00077 template<unsigned DIM>
00078 Vec CellBasedSimulationWithPdes<DIM>::GetCurrentPdeSolution(unsigned pdeIndex)
00079 {
00080     assert(pdeIndex<mPdeAndBcCollection.size());
00081     return mPdeAndBcCollection[pdeIndex]->GetSolution();
00082 }
00083 
00085 //                          Setup/AfterSolve methods                        //
00087 
00088 template<unsigned DIM>
00089 void CellBasedSimulationWithPdes<DIM>::WriteVisualizerSetupFile()
00090 {
00091     for (unsigned i=0; i<this->mForceCollection.size(); i++)
00092     {
00093         if (dynamic_cast<AbstractTwoBodyInteractionForce<DIM>*>(this->mForceCollection[i]))
00094         {
00095             double cutoff = (static_cast<AbstractTwoBodyInteractionForce<DIM>*>(this->mForceCollection[i]))->GetCutOffLength();
00096             *(this->mpVizSetupFile) << "Cutoff\t" << cutoff << "\n";
00097         }
00098     }
00099 }
00100 
00101 template<unsigned DIM>
00102 void CellBasedSimulationWithPdes<DIM>::SetupSolve()
00103 {
00104     if (mpCoarsePdeMesh != NULL)
00105     {
00106         InitialiseCoarsePdeMesh();
00107     }
00108 
00109     // We must initially have at least one cell in the cell-based simulation
00110     assert(this->mrCellPopulation.GetNumRealCells() != 0);
00111 
00112     SetupWritePdeSolution();
00113     double current_time = SimulationTime::Instance()->GetTime();
00114     WritePdeSolution(current_time);
00115 }
00116 
00117 template<unsigned DIM>
00118 void CellBasedSimulationWithPdes<DIM>::SetupWritePdeSolution()
00119 {
00120     OutputFileHandler output_file_handler(this->mSimulationOutputDirectory+"/", false);
00121     if (PetscTools::AmMaster())
00122     {
00123         mpVizPdeSolutionResultsFile = output_file_handler.OpenOutputFile("results.vizpdesolution");
00124         *this->mpVizSetupFile << "PDE \n";
00125         if (mWriteAverageRadialPdeSolution)
00126         {
00127             mpAverageRadialPdeSolutionResultsFile = output_file_handler.OpenOutputFile("radial_dist.dat");
00128         }
00129     }
00130 }
00131 
00132 template<unsigned DIM>
00133 void CellBasedSimulationWithPdes<DIM>::UseCoarsePdeMesh(double coarseGrainScaleFactor)
00134 {
00135     assert(!mPdeAndBcCollection.empty());
00136     for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00137     {
00138         assert(mPdeAndBcCollection[pde_index]->HasAveragedSourcePde());
00139     }
00140 
00141     CreateCoarsePdeMesh(coarseGrainScaleFactor);
00142 }
00143 
00144 template<unsigned DIM>
00145 void CellBasedSimulationWithPdes<DIM>::CreateCoarsePdeMesh(double coarseGrainScaleFactor)
00146 {
00147     EXCEPTION("This method is only implemented in 2D");
00148 }
00149 
00156 template<>
00157 void CellBasedSimulationWithPdes<2>::CreateCoarsePdeMesh(double coarseGrainScaleFactor)
00158 {
00159     // Create coarse PDE mesh (can use a larger mesh if required, e.g. disk_984_elements)
00160     TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/disk_522_elements");
00161     mpCoarsePdeMesh = new TetrahedralMesh<2,2>;
00162     mpCoarsePdeMesh->ConstructFromMeshReader(mesh_reader);
00163 
00164     // Find centre of cell population
00165     c_vector<double,2> centre_of_cell_population = zero_vector<double>(2);
00166     for (AbstractCellPopulation<2>::Iterator cell_iter = this->mrCellPopulation.Begin();
00167         cell_iter != this->mrCellPopulation.End();
00168         ++cell_iter)
00169     {
00170         centre_of_cell_population += this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
00171     }
00172     centre_of_cell_population /= this->mrCellPopulation.GetNumRealCells();
00173 
00174     // Find max radius of cell population
00175     double max_cell_population_radius = 0.0;
00176     for (AbstractCellPopulation<2>::Iterator cell_iter = this->mrCellPopulation.Begin();
00177         cell_iter != this->mrCellPopulation.End();
00178         ++cell_iter)
00179     {
00180         double radius = norm_2(centre_of_cell_population - this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter));
00181         if (radius > max_cell_population_radius)
00182         {
00183             max_cell_population_radius = radius;
00184         }
00185     }
00186 
00187     // Find centre of coarse PDE mesh
00188     c_vector<double,2> centre_of_coarse_mesh = zero_vector<double>(2);
00189     for (unsigned i=0; i<mpCoarsePdeMesh->GetNumNodes(); i++)
00190     {
00191         centre_of_coarse_mesh += mpCoarsePdeMesh->GetNode(i)->rGetLocation();
00192     }
00193     centre_of_coarse_mesh /= mpCoarsePdeMesh->GetNumNodes();
00194 
00195     // Find max radius of coarse PDE mesh
00196     double max_mesh_radius = 0.0;
00197     for (unsigned i=0; i<mpCoarsePdeMesh->GetNumNodes(); i++)
00198     {
00199         double radius = norm_2(centre_of_coarse_mesh - mpCoarsePdeMesh->GetNode(i)->rGetLocation());
00200         if (radius > max_mesh_radius)
00201         {
00202             max_mesh_radius = radius;
00203         }
00204     }
00205 
00206     // Translate centre of coarse PDE mesh to the origin
00207     mpCoarsePdeMesh->Translate(-centre_of_coarse_mesh[0], -centre_of_coarse_mesh[1]);
00208 
00209     // Scale coarse PDE mesh
00210     double scale_factor = (max_cell_population_radius/max_mesh_radius)*coarseGrainScaleFactor;
00211     mpCoarsePdeMesh->Scale(scale_factor, scale_factor);
00212 
00213     // Translate centre of coarse PDE mesh to centre of the cell population
00214     mpCoarsePdeMesh->Translate(centre_of_cell_population[0], centre_of_cell_population[1]);
00215 }
00216 
00217 template<unsigned DIM>
00218 void CellBasedSimulationWithPdes<DIM>::InitialiseCoarsePdeMesh()
00219 {
00220     mCellPdeElementMap.clear();
00221 
00222     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mrCellPopulation.Begin();
00223         cell_iter != this->mrCellPopulation.End();
00224         ++cell_iter)
00225     {
00226         // Find the element of mpCoarsePdeMesh that contains this cell
00227         const ChastePoint<DIM>& r_position_of_cell = this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
00228         unsigned elem_index = mpCoarsePdeMesh->GetContainingElementIndex(r_position_of_cell);
00229         mCellPdeElementMap[*cell_iter] = elem_index;
00230     }
00231 }
00232 
00233 template<unsigned DIM>
00234 void CellBasedSimulationWithPdes<DIM>::AfterSolve()
00235 {
00236     if (this->mrCellPopulation.GetNumRealCells() != 0 && PetscTools::AmMaster())
00237     {
00238         mpVizPdeSolutionResultsFile->close();
00239 
00240         if (mWriteAverageRadialPdeSolution)
00241         {
00242             WriteAverageRadialPdeSolution(SimulationTime::Instance()->GetTime(), mNumRadialIntervals);
00243             mpAverageRadialPdeSolutionResultsFile->close();
00244         }
00245     }
00246 }
00247 
00249 //                             PostSolve methods                            //
00251 
00252 template<unsigned DIM>
00253 void CellBasedSimulationWithPdes<DIM>::SolvePde()
00254 {
00255     if (mpCoarsePdeMesh != NULL)
00256     {
00257         SolvePdeUsingCoarseMesh();
00258         return;
00259     }
00260 
00261     assert(!mPdeAndBcCollection.empty());
00262 
00263     for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00264     {
00265         assert(mPdeAndBcCollection[pde_index]);
00266         assert(mPdeAndBcCollection[pde_index]->HasAveragedSourcePde() == false);
00267     }
00268 
00269     // Note: If not using a coarse PDE mesh, we MUST be using a MeshBasedCellPopulation
00270     // Make sure the mesh is in a nice state
00271     this->mrCellPopulation.Update();
00272 
00273     TetrahedralMesh<DIM,DIM>& r_mesh = static_cast<MeshBasedCellPopulation<DIM>*>(&(this->mrCellPopulation))->rGetMesh();
00274     CellwiseData<DIM>::Instance()->ReallocateMemory();
00275 
00276     // Loop over elements of mPdeAndBcCollection
00277     for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00278     {
00279         // Get pointer to this PdeAndBoundaryConditions object
00280         PdeAndBoundaryConditions<DIM>* p_pde_and_bc = mPdeAndBcCollection[pde_index];
00281 
00282         // Set up boundary conditions
00283         AbstractBoundaryCondition<DIM>* p_bc = p_pde_and_bc->GetBoundaryCondition();
00284 
00285         BoundaryConditionsContainer<DIM,DIM,1> bcc(false);
00286         if (p_pde_and_bc->IsNeumannBoundaryCondition())
00287         {
00288             for (typename TetrahedralMesh<DIM,DIM>::BoundaryElementIterator elem_iter = r_mesh.GetBoundaryElementIteratorBegin();
00289                  elem_iter != r_mesh.GetBoundaryElementIteratorEnd();
00290                  ++elem_iter)
00291             {
00292                 bcc.AddNeumannBoundaryCondition(*elem_iter, p_bc);
00293             }
00294         }
00295         else // assuming that if it's not Neumann, then it's Dirichlet
00296         {
00297             for (typename TetrahedralMesh<DIM,DIM>::BoundaryNodeIterator node_iter = r_mesh.GetBoundaryNodeIteratorBegin();
00298                  node_iter != r_mesh.GetBoundaryNodeIteratorEnd();
00299                  ++node_iter)
00300             {
00301                 bcc.AddDirichletBoundaryCondition(*node_iter, p_bc);
00302             }
00303         }
00304 
00305         /*
00306          * Set up solver. This is a purpose-made elliptic solver which must
00307          * interpolate contributions to source terms from nodes onto Gauss points,
00308          * because the PDE solution is only stored at the cells (nodes).
00309          */
00310         CellBasedSimulationWithPdesAssembler<DIM> solver(&r_mesh, p_pde_and_bc->GetPde(), &bcc);
00311 
00312         PetscInt size_of_soln_previous_step = 0;
00313 
00314         if (p_pde_and_bc->GetSolution())
00315         {
00316             VecGetSize(p_pde_and_bc->GetSolution(), &size_of_soln_previous_step);
00317         }
00318         if (size_of_soln_previous_step == (int)r_mesh.GetNumNodes())
00319         {
00320             // We make an initial guess which gets copied by the Solve method of
00321             // SimpleLinearSolver, so we need to delete it too.
00322             Vec initial_guess;
00323             VecDuplicate(p_pde_and_bc->GetSolution(), &initial_guess);
00324             VecCopy(p_pde_and_bc->GetSolution(), initial_guess);
00325 
00326             // Use current solution as the initial guess
00327             p_pde_and_bc->DestroySolution(); // Solve method makes its own mCurrentPdeSolution
00328             p_pde_and_bc->SetSolution(solver.Solve(initial_guess));
00329             VecDestroy(initial_guess);
00330         }
00331         else
00332         {
00333             if (p_pde_and_bc->GetSolution())
00334             {
00335                 assert(size_of_soln_previous_step != 0);
00336                 p_pde_and_bc->DestroySolution();
00337             }
00338             p_pde_and_bc->SetSolution(solver.Solve());
00339         }
00340 
00341         ReplicatableVector solution_repl(p_pde_and_bc->GetSolution());
00342 
00343         // Update cellwise data
00344         for (unsigned i=0; i<r_mesh.GetNumNodes(); i++)
00345         {
00346             double solution = solution_repl[i];
00347             unsigned index = r_mesh.GetNode(i)->GetIndex();
00348             CellwiseData<DIM>::Instance()->SetValue(solution, index, pde_index);
00349         }
00350     }
00351 }
00352 
00353 template<unsigned DIM>
00354 void CellBasedSimulationWithPdes<DIM>::SolvePdeUsingCoarseMesh()
00355 {
00356     assert(!mPdeAndBcCollection.empty());
00357 
00358     for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00359     {
00360         assert(mPdeAndBcCollection[pde_index]);
00361         assert(mPdeAndBcCollection[pde_index]->HasAveragedSourcePde() == true);
00362     }
00363 
00364     TetrahedralMesh<DIM,DIM>& r_mesh = *mpCoarsePdeMesh;
00365     CellwiseData<DIM>::Instance()->ReallocateMemory();
00366 
00367     // Loop over cells and calculate centre of distribution
00368     c_vector<double, DIM> centre = zero_vector<double>(DIM);
00369     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mrCellPopulation.Begin();
00370         cell_iter != this->mrCellPopulation.End();
00371         ++cell_iter)
00372     {
00373         centre += this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
00374     }
00375     centre /= this->mrCellPopulation.GetNumRealCells();
00376 
00377     // Find max radius
00378     double max_radius = 0.0;
00379     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mrCellPopulation.Begin();
00380         cell_iter != this->mrCellPopulation.End();
00381         ++cell_iter)
00382     {
00383         double radius = norm_2(centre - this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter));
00384         if (radius > max_radius)
00385         {
00386             max_radius = radius;
00387         }
00388     }
00389 
00390     // Loop over elements of mPdeAndBcCollection
00391     for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
00392     {
00393         // Get pointer to this PdeAndBoundaryConditions object
00394         PdeAndBoundaryConditions<DIM>* p_pde_and_bc = mPdeAndBcCollection[pde_index];
00395 
00396         // Set up boundary conditions
00397         AbstractBoundaryCondition<DIM>* p_bc = p_pde_and_bc->GetBoundaryCondition();
00398 
00399         BoundaryConditionsContainer<DIM,DIM,1> bcc(false);
00400         if (p_pde_and_bc->IsNeumannBoundaryCondition())
00401         {
00402             EXCEPTION("Neumann BCs not yet implemented when using a coarse PDE mesh");
00403         }
00404         else
00405         {
00406             // Get the set of coarse element indices that contain cells
00407             std::set<unsigned> coarse_element_indices_in_map;
00408             for (typename MeshBasedCellPopulation<DIM>::Iterator cell_iter = this->mrCellPopulation.Begin();
00409                 cell_iter != this->mrCellPopulation.End();
00410                 ++cell_iter)
00411             {
00412                 coarse_element_indices_in_map.insert(mCellPdeElementMap[*cell_iter]);
00413             }
00414 
00415             // Find the node indices that associated with elements whose
00416             // indices are NOT in the set coarse_element_indices_in_map
00417             std::set<unsigned> coarse_mesh_boundary_node_indices;
00418 
00419             for (unsigned i=0; i<r_mesh.GetNumElements(); i++)
00420             {
00421                 // If the element index is NOT in the set...
00422                 if (coarse_element_indices_in_map.find(i) == coarse_element_indices_in_map.end())
00423                 {
00424                     // ... then get the element...
00425                     Element<DIM,DIM>* p_element = r_mesh.GetElement(i);
00426 
00427                     // ... and add its associated nodes to coarse_mesh_boundary_node_indices
00428                     for (unsigned local_index=0; local_index<DIM+1; local_index++)
00429                     {
00430                         unsigned node_index = p_element->GetNodeGlobalIndex(local_index);
00431                         coarse_mesh_boundary_node_indices.insert(node_index);
00432                     }
00433                 }
00434             }
00435 
00436             // Apply boundary condition to the nodes in the set coarse_mesh_boundary_node_indices
00437             for (std::set<unsigned>::iterator iter = coarse_mesh_boundary_node_indices.begin();
00438                  iter != coarse_mesh_boundary_node_indices.end();
00439                  ++iter)
00440             {
00441                 bcc.AddDirichletBoundaryCondition(r_mesh.GetNode(*iter), p_bc, 0, false);
00442             }
00443         }
00444 
00445         PetscInt size_of_soln_previous_step = 0;
00446 
00447         if (p_pde_and_bc->GetSolution())
00448         {
00449             VecGetSize(p_pde_and_bc->GetSolution(), &size_of_soln_previous_step);
00450         }
00451 
00452         p_pde_and_bc->SetUpSourceTermsForAveragedSourcePde(mpCoarsePdeMesh);
00453 
00454         SimpleLinearEllipticSolver<DIM,DIM> solver(mpCoarsePdeMesh, p_pde_and_bc->GetPde(), &bcc);
00455 
00456         if (size_of_soln_previous_step == (int)r_mesh.GetNumNodes())
00457         {
00458             // We make an initial guess which gets copied by the Solve method of
00459             // SimpleLinearSolver, so we need to delete it too.
00460             Vec initial_guess;
00461             VecDuplicate(p_pde_and_bc->GetSolution(), &initial_guess);
00462             VecCopy(p_pde_and_bc->GetSolution(), initial_guess);
00463 
00464             // Use current solution as the initial guess
00465             p_pde_and_bc->DestroySolution(); // Solve method makes its own mCurrentPdeSolution
00466             p_pde_and_bc->SetSolution(solver.Solve(initial_guess));
00467             VecDestroy(initial_guess);
00468         }
00469         else
00470         {
00471             /*
00472              * Eventually we will enable the coarse PDE mesh to change size, for example
00473              * in the case of a spheroid that grows a lot (see #630). In this case we should
00474              * uncomment the following code.
00475              *
00476             if (mCurrentPdeSolution)
00477             {
00478                 assert(0);
00479                 VecDestroy(mCurrentPdeSolution);
00480             }
00481             *
00482             */
00483             p_pde_and_bc->SetSolution(solver.Solve());
00484         }
00485 
00486         /*
00487          * Update cellwise data - since the cells are not nodes on the coarse
00488          * mesh, we have to interpolate from the nodes of the coarse mesh onto
00489          * the cell locations.
00490          */
00491         ReplicatableVector solution_repl(p_pde_and_bc->GetSolution());
00492 
00493         for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mrCellPopulation.Begin();
00494             cell_iter != this->mrCellPopulation.End();
00495             ++cell_iter)
00496         {
00497             // Find coarse mesh element containing cell
00498             unsigned elem_index = FindCoarseElementContainingCell(*cell_iter);
00499 
00500             Element<DIM,DIM>* p_element = mpCoarsePdeMesh->GetElement(elem_index);
00501 
00502             const ChastePoint<DIM>& r_position_of_cell = this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
00503 
00504             c_vector<double,DIM+1> weights = p_element->CalculateInterpolationWeights(r_position_of_cell);
00505 
00506             double interpolated_solution = 0.0;
00507             for (unsigned i=0; i<DIM+1/*num_nodes*/; i++)
00508             {
00509                 double nodal_value = solution_repl[ p_element->GetNodeGlobalIndex(i) ];
00510                 interpolated_solution += nodal_value*weights(i);
00511             }
00512 
00513             unsigned index = this->mrCellPopulation.GetLocationIndexUsingCell(*cell_iter);
00514             CellwiseData<DIM>::Instance()->SetValue(interpolated_solution, index,pde_index);
00515         }
00516     }
00517 }
00518 
00519 template<unsigned DIM>
00520 unsigned CellBasedSimulationWithPdes<DIM>::FindCoarseElementContainingCell(CellPtr pCell)
00521 {
00522     // Get containing element at last timestep from mCellPdeElementMap
00523     unsigned old_element_index = mCellPdeElementMap[pCell];
00524 
00525     // Create a std::set of guesses for the current containing element
00526     std::set<unsigned> test_elements;
00527     test_elements.insert(old_element_index);
00528 
00529     Element<DIM,DIM>* p_element = mpCoarsePdeMesh->GetElement(old_element_index);
00530 
00531     for (unsigned local_index=0; local_index<DIM+1; local_index++)
00532     {
00533         std::set<unsigned> element_indices = p_element->GetNode(local_index)->rGetContainingElementIndices();
00534 
00535         for (std::set<unsigned>::iterator iter = element_indices.begin();
00536              iter != element_indices.end();
00537              ++iter)
00538         {
00539             test_elements.insert(*iter);
00540         }
00541     }
00542 
00543     // Find new element, using the previous one as a guess
00544     const ChastePoint<DIM>& r_cell_position = this->mrCellPopulation.GetLocationOfCellCentre(pCell);
00545     unsigned new_element_index = mpCoarsePdeMesh->GetContainingElementIndex(r_cell_position, false, test_elements);
00546 
00547     // Update mCellPdeElementMap
00548     mCellPdeElementMap[pCell] = new_element_index;
00549 
00550     return new_element_index;
00551 }
00552 
00553 template<unsigned DIM>
00554 void CellBasedSimulationWithPdes<DIM>::PostSolve()
00555 {
00556     SolvePde();
00557 
00558     // Save results to file
00559     SimulationTime* p_time = SimulationTime::Instance();
00560 
00561     double time_next_step = p_time->GetTime() + p_time->GetTimeStep();
00562 
00563     if ((p_time->GetTimeStepsElapsed()+1)%this->mSamplingTimestepMultiple == 0)
00564     {
00565         WritePdeSolution(time_next_step);
00566     }
00567 
00568 #define COVERAGE_IGNORE
00569     if (mWriteDailyAverageRadialPdeSolution)
00570     {
00572         unsigned num_timesteps_per_day = (unsigned) (DBL_EPSILON + 24/SimulationTime::Instance()->GetTimeStep());
00573 
00574         if ((p_time->GetTimeStepsElapsed()+1) % num_timesteps_per_day == 0)
00575         {
00576             WriteAverageRadialPdeSolution(time_next_step, mNumRadialIntervals);
00577         }
00578     }
00579 #undef COVERAGE_IGNORE
00580 
00581 }
00582 
00584 //                             Output methods                               //
00586 
00587 template<unsigned DIM>
00588 void CellBasedSimulationWithPdes<DIM>::WritePdeSolution(double time)
00589 {
00590     if (PetscTools::AmMaster())
00591     {
00592         // Since there are no ghost nodes, the number of nodes must equal the number of real cells
00593         assert(this->mrCellPopulation.GetNumNodes() == this->mrCellPopulation.GetNumRealCells());
00594 
00595         (*mpVizPdeSolutionResultsFile) << time << "\t";
00596 
00597         for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mrCellPopulation.Begin();
00598              cell_iter != this->mrCellPopulation.End();
00599              ++cell_iter)
00600         {
00601             unsigned global_index = this->mrCellPopulation.GetLocationIndexUsingCell(*cell_iter);
00602             (*mpVizPdeSolutionResultsFile) << global_index << " ";
00603 
00604             const c_vector<double,DIM>& position = this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
00605             for (unsigned i=0; i<DIM; i++)
00606             {
00607                 (*mpVizPdeSolutionResultsFile) << position[i] << " ";
00608             }
00609 
00610             double solution = CellwiseData<DIM>::Instance()->GetValue(*cell_iter);
00611             (*mpVizPdeSolutionResultsFile) << solution << " ";
00612         }
00613         (*mpVizPdeSolutionResultsFile) << "\n";
00614     }
00615 }
00616 
00617 template<unsigned DIM>
00618 void CellBasedSimulationWithPdes<DIM>::SetWriteAverageRadialPdeSolution(unsigned numRadialIntervals, bool writeDailyResults)
00619 {
00620     mWriteAverageRadialPdeSolution = true;
00621     mNumRadialIntervals = numRadialIntervals;
00622     mWriteDailyAverageRadialPdeSolution = writeDailyResults;
00623 }
00624 
00625 template<unsigned DIM>
00626 void CellBasedSimulationWithPdes<DIM>::WriteAverageRadialPdeSolution(double time, unsigned numRadialIntervals)
00627 {
00628     (*mpAverageRadialPdeSolutionResultsFile) << time << " ";
00629 
00630     // Calculate the centre of the cell population
00631     c_vector<double,DIM> centre = zero_vector<double>(DIM);
00632     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mrCellPopulation.Begin();
00633          cell_iter != this->mrCellPopulation.End();
00634          ++cell_iter)
00635     {
00636        centre += (this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter));
00637     }
00638     centre /= ((double) this->mrCellPopulation.GetNumNodes());
00639 
00640     // Calculate the distance between each node and the centre of the cell population, as well as the maximum of these
00641     std::map<double, CellPtr> distance_cell_map;
00642     double max_distance_from_centre = 0.0;
00643     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mrCellPopulation.Begin();
00644          cell_iter != this->mrCellPopulation.End();
00645          ++cell_iter)
00646     {
00647         double distance = norm_2(this->mrCellPopulation.GetLocationOfCellCentre(*cell_iter) - centre);
00648         distance_cell_map[distance] = *cell_iter;
00649 
00650         if (distance > max_distance_from_centre)
00651         {
00652             max_distance_from_centre = distance;
00653         }
00654     }
00655 
00656     // Create vector of radius intervals
00657     std::vector<double> radius_intervals;
00658     for (unsigned i=0; i<numRadialIntervals; i++)
00659     {
00660         double upper_radius = max_distance_from_centre*((double) i+1)/((double) numRadialIntervals);
00661         radius_intervals.push_back(upper_radius);
00662     }
00663 
00664     // Calculate PDE solution in each radial interval
00665     double lower_radius = 0.0;
00666     for (unsigned i=0; i<numRadialIntervals; i++)
00667     {
00668         unsigned counter = 0;
00669         double average_solution = 0.0;
00670 
00671         for (std::map<double, CellPtr>::iterator iter = distance_cell_map.begin();
00672              iter != distance_cell_map.end();
00673              ++iter)
00674         {
00675             if (iter->first > lower_radius && iter->first <= radius_intervals[i])
00676             {
00677                 average_solution += CellwiseData<DIM>::Instance()->GetValue(iter->second);
00678                 counter++;
00679             }
00680         }
00681         if (counter > 0)
00682         {
00683             average_solution /= (double) counter;
00684         }
00685 
00686         // Write results to file
00687         (*mpAverageRadialPdeSolutionResultsFile) << radius_intervals[i] << " " << average_solution << " ";
00688         lower_radius = radius_intervals[i];
00689     }
00690     (*mpAverageRadialPdeSolutionResultsFile) << "\n";
00691 }
00692 
00693 template<unsigned DIM>
00694 void CellBasedSimulationWithPdes<DIM>::OutputSimulationParameters(out_stream& rParamsFile)
00695 {
00696     *rParamsFile << "\t\t<WriteAverageRadialPdeSolution>" << mWriteAverageRadialPdeSolution << "</WriteAverageRadialPdeSolution>\n";
00697     *rParamsFile << "\t\t<WriteDailyAverageRadialPdeSolution>" << mWriteDailyAverageRadialPdeSolution << "</WriteDailyAverageRadialPdeSolution>\n";
00698 
00699     // Call method on direct parent class
00700     CellBasedSimulation<DIM>::OutputSimulationParameters(rParamsFile);
00701 }
00702 
00704 // Explicit instantiation
00706 
00707 template class CellBasedSimulationWithPdes<1>;
00708 template class CellBasedSimulationWithPdes<2>;
00709 template class CellBasedSimulationWithPdes<3>;
00710 
00711 // Serialization for Boost >= 1.36
00712 #include "SerializationExportWrapperForCpp.hpp"
00713 EXPORT_TEMPLATE_CLASS_SAME_DIMS(CellBasedSimulationWithPdes)

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5