00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "MeshBasedCellPopulation.hpp"
00030 #include "CellwiseData.hpp"
00031 #include "TrianglesMeshWriter.hpp"
00032 #include "VtkMeshWriter.hpp"
00033 #include "CellBasedEventHandler.hpp"
00034 #include "ApoptoticCellProperty.hpp"
00035 #include "Cylindrical2dMesh.hpp"
00036 #include "NodesOnlyMesh.hpp"
00037
00038
00039 template<unsigned DIM>
00040 MeshBasedCellPopulation<DIM>::MeshBasedCellPopulation(MutableMesh<DIM, DIM>& rMesh,
00041 std::vector<CellPtr>& rCells,
00042 const std::vector<unsigned> locationIndices,
00043 bool deleteMesh,
00044 bool validate)
00045 : AbstractCentreBasedCellPopulation<DIM>(rCells, locationIndices),
00046 mrMesh(rMesh),
00047 mpVoronoiTessellation(NULL),
00048 mDeleteMesh(deleteMesh),
00049 mUseAreaBasedDampingConstant(false),
00050 mAreaBasedDampingConstantParameter(0.1),
00051 mOutputVoronoiData(false),
00052 mOutputCellPopulationVolumes(false),
00053 mWriteVtkAsPoints(false)
00054 {
00055
00056 assert(this->mCells.size() <= mrMesh.GetNumNodes());
00057
00058 this->mCellPopulationContainsMesh = true;
00059
00060 if (validate)
00061 {
00062 Validate();
00063 }
00064 }
00065
00066 template<unsigned DIM>
00067 MeshBasedCellPopulation<DIM>::MeshBasedCellPopulation(MutableMesh<DIM, DIM>& rMesh)
00068 : mrMesh(rMesh)
00069 {
00070 this->mCellPopulationContainsMesh = true;
00071 mpVoronoiTessellation = NULL;
00072 mDeleteMesh = true;
00073 }
00074
00075 template<unsigned DIM>
00076 MeshBasedCellPopulation<DIM>::~MeshBasedCellPopulation()
00077 {
00078 delete mpVoronoiTessellation;
00079
00080 if (mDeleteMesh)
00081 {
00082 delete &mrMesh;
00083 }
00084 }
00085
00086 template<unsigned DIM>
00087 bool MeshBasedCellPopulation<DIM>::UseAreaBasedDampingConstant()
00088 {
00089 return mUseAreaBasedDampingConstant;
00090 }
00091
00092 template<unsigned DIM>
00093 void MeshBasedCellPopulation<DIM>::SetAreaBasedDampingConstant(bool useAreaBasedDampingConstant)
00094 {
00095 assert(DIM==2);
00096 mUseAreaBasedDampingConstant = useAreaBasedDampingConstant;
00097 }
00098
00099 template<unsigned DIM>
00100 unsigned MeshBasedCellPopulation<DIM>::AddNode(Node<DIM>* pNewNode)
00101 {
00102 return mrMesh.AddNode(pNewNode);
00103 }
00104
00105 template<unsigned DIM>
00106 void MeshBasedCellPopulation<DIM>::SetNode(unsigned nodeIndex, ChastePoint<DIM>& rNewLocation)
00107 {
00108 mrMesh.SetNode(nodeIndex, rNewLocation, false);
00109 }
00110
00111 template<unsigned DIM>
00112 double MeshBasedCellPopulation<DIM>::GetDampingConstant(unsigned nodeIndex)
00113 {
00114 double damping_multiplier = AbstractCentreBasedCellPopulation<DIM>::GetDampingConstant(nodeIndex);
00115
00116 if (mUseAreaBasedDampingConstant)
00117 {
00127 assert(DIM==2);
00128
00129 double rest_length = 1.0;
00130 double d0 = mAreaBasedDampingConstantParameter;
00131
00137 double d1 = 2.0*(1.0 - d0)/(sqrt(3)*rest_length*rest_length);
00138
00139 double area_cell = GetVolumeOfVoronoiElement(nodeIndex);
00140
00146 assert(area_cell < 1000);
00147
00148 damping_multiplier = d0 + area_cell*d1;
00149 }
00150
00151 return damping_multiplier;
00152 }
00153
00154 template<unsigned DIM>
00155 void MeshBasedCellPopulation<DIM>::Validate()
00156 {
00157 std::vector<bool> validated_node = std::vector<bool>(this->GetNumNodes(), false);
00158
00159 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
00160 {
00161 unsigned node_index = GetLocationIndexUsingCell(*cell_iter);
00162 validated_node[node_index] = true;
00163 }
00164
00165 for (unsigned i=0; i<validated_node.size(); i++)
00166 {
00167 if (!validated_node[i])
00168 {
00169 std::stringstream ss;
00170 ss << "Node " << i << " does not appear to have a cell associated with it";
00171 EXCEPTION(ss.str());
00172 }
00173 }
00174 }
00175
00176 template<unsigned DIM>
00177 MutableMesh<DIM, DIM>& MeshBasedCellPopulation<DIM>::rGetMesh()
00178 {
00179 return mrMesh;
00180 }
00181
00182 template<unsigned DIM>
00183 const MutableMesh<DIM, DIM>& MeshBasedCellPopulation<DIM>::rGetMesh() const
00184 {
00185 return mrMesh;
00186 }
00187
00188 template<unsigned DIM>
00189 unsigned MeshBasedCellPopulation<DIM>::RemoveDeadCells()
00190 {
00191 unsigned num_removed = 0;
00192 for (std::list<CellPtr>::iterator it = this->mCells.begin();
00193 it != this->mCells.end();
00194 ++it)
00195 {
00196 if ((*it)->IsDead())
00197 {
00198
00199 std::vector<const std::pair<CellPtr,CellPtr>*> pairs_to_remove;
00200 for (std::set<std::pair<CellPtr,CellPtr> >::iterator it1 = mMarkedSprings.begin();
00201 it1 != mMarkedSprings.end();
00202 ++it1)
00203 {
00204 const std::pair<CellPtr,CellPtr>& r_pair = *it1;
00205
00206 for (unsigned i=0; i<2; i++)
00207 {
00208 CellPtr p_cell = (i==0 ? r_pair.first : r_pair.second);
00209
00210 if (p_cell == *it)
00211 {
00212
00213 pairs_to_remove.push_back(&r_pair);
00214 break;
00215 }
00216 }
00217 }
00218
00219
00220 for (std::vector<const std::pair<CellPtr,CellPtr>* >::iterator pair_it = pairs_to_remove.begin();
00221 pair_it != pairs_to_remove.end();
00222 ++pair_it)
00223 {
00224 mMarkedSprings.erase(**pair_it);
00225 }
00226
00227
00228 num_removed++;
00229 mrMesh.DeleteNodePriorToReMesh(this->mCellLocationMap[(*it).get()]);
00230
00231
00232 unsigned location_index_of_removed_node = this->mCellLocationMap[(*it).get()];
00233 this->mCellLocationMap.erase((*it).get());
00234 this->mLocationCellMap.erase(location_index_of_removed_node);
00235
00236
00237 it = this->mCells.erase(it);
00238 --it;
00239 }
00240 }
00241
00242 return num_removed;
00243 }
00244
00245
00246 template<unsigned DIM>
00247 void MeshBasedCellPopulation<DIM>::Update(bool hasHadBirthsOrDeaths)
00248 {
00249 NodeMap map(mrMesh.GetNumAllNodes());
00250 mrMesh.ReMesh(map);
00251
00252 if (!map.IsIdentityMap())
00253 {
00254 UpdateGhostNodesAfterReMesh(map);
00255
00256
00257 std::map<Cell*, unsigned> old_map = this->mCellLocationMap;
00258
00259
00260 this->mLocationCellMap.clear();
00261 this->mCellLocationMap.clear();
00262
00263 for (std::list<CellPtr>::iterator it = this->mCells.begin();
00264 it != this->mCells.end();
00265 ++it)
00266 {
00267 unsigned old_node_index = old_map[(*it).get()];
00268
00269
00270 assert(!map.IsDeleted(old_node_index));
00271
00272 unsigned new_node_index = map.GetNewIndex(old_node_index);
00273 this->mLocationCellMap[new_node_index] = *it;
00274 this->mCellLocationMap[(*it).get()] = new_node_index;
00275 }
00276
00277 this->Validate();
00278 }
00279
00280
00281 std::vector<const std::pair<CellPtr,CellPtr>*> springs_to_remove;
00282 for (std::set<std::pair<CellPtr,CellPtr> >::iterator spring_it = mMarkedSprings.begin();
00283 spring_it != mMarkedSprings.end();
00284 ++spring_it)
00285 {
00286 CellPtr p_cell_1 = spring_it->first;
00287 CellPtr p_cell_2 = spring_it->second;
00288 Node<DIM>* p_node_1 = this->GetNodeCorrespondingToCell(p_cell_1);
00289 Node<DIM>* p_node_2 = this->GetNodeCorrespondingToCell(p_cell_2);
00290
00291 bool joined = false;
00292
00293
00294 std::set<unsigned> node2_elements = p_node_2->rGetContainingElementIndices();
00295 for (typename Node<DIM>::ContainingElementIterator elem_iter = p_node_1->ContainingElementsBegin();
00296 elem_iter != p_node_1->ContainingElementsEnd();
00297 ++elem_iter)
00298 {
00299 if (node2_elements.find(*elem_iter) != node2_elements.end())
00300 {
00301 joined = true;
00302 break;
00303 }
00304 }
00305
00306
00307 if (!joined)
00308 {
00309 springs_to_remove.push_back(&(*spring_it));
00310 }
00311 }
00312
00313
00314 for (std::vector<const std::pair<CellPtr,CellPtr>* >::iterator spring_it = springs_to_remove.begin();
00315 spring_it != springs_to_remove.end();
00316 ++spring_it)
00317 {
00318 mMarkedSprings.erase(**spring_it);
00319 }
00320
00321
00322 if (DIM==2 || DIM==3)
00323 {
00324 CellBasedEventHandler::BeginEvent(CellBasedEventHandler::TESSELLATION);
00325 if (mUseAreaBasedDampingConstant || mOutputVoronoiData || mOutputCellPopulationVolumes || this->mOutputCellVolumes)
00326 {
00327 CreateVoronoiTessellation();
00328 }
00329 CellBasedEventHandler::EndEvent(CellBasedEventHandler::TESSELLATION);
00330 }
00331
00332 mrMesh.SetMeshHasChangedSinceLoading();
00333 }
00334
00335 template<unsigned DIM>
00336 Node<DIM>* MeshBasedCellPopulation<DIM>::GetNode(unsigned index)
00337 {
00338 return mrMesh.GetNode(index);
00339 }
00340
00341 template<unsigned DIM>
00342 unsigned MeshBasedCellPopulation<DIM>::GetNumNodes()
00343 {
00344 return mrMesh.GetNumAllNodes();
00345 }
00346
00347 template<unsigned DIM>
00348 void MeshBasedCellPopulation<DIM>::UpdateGhostNodesAfterReMesh(NodeMap& rMap)
00349 {
00350 }
00351
00352 template<unsigned DIM>
00353 CellPtr MeshBasedCellPopulation<DIM>::AddCell(CellPtr pNewCell, const c_vector<double,DIM>& rCellDivisionVector, CellPtr pParentCell)
00354 {
00355 assert(pNewCell);
00356 assert(pParentCell);
00357
00358
00359 CellPtr p_created_cell = AbstractCentreBasedCellPopulation<DIM>::AddCell(pNewCell, rCellDivisionVector, pParentCell);
00360 assert(p_created_cell == pNewCell);
00361
00362
00363 std::pair<CellPtr,CellPtr> cell_pair = CreateCellPair(pParentCell, p_created_cell);
00364 MarkSpring(cell_pair);
00365
00366
00367 return p_created_cell;
00368 }
00369
00371
00373
00374 template<unsigned DIM>
00375 void MeshBasedCellPopulation<DIM>::CreateOutputFiles(const std::string& rDirectory, bool cleanOutputDirectory)
00376 {
00377 AbstractCentreBasedCellPopulation<DIM>::CreateOutputFiles(rDirectory, cleanOutputDirectory);
00378
00379 OutputFileHandler output_file_handler(rDirectory, cleanOutputDirectory);
00380 mpVizElementsFile = output_file_handler.OpenOutputFile("results.vizelements");
00381
00382 if (mOutputVoronoiData)
00383 {
00384 mpVoronoiFile = output_file_handler.OpenOutputFile("voronoi.dat");
00385 }
00386 if (mOutputCellPopulationVolumes)
00387 {
00388 mpCellPopulationVolumesFile = output_file_handler.OpenOutputFile("cellpopulationareas.dat");
00389 }
00390 if (this->mOutputCellVolumes)
00391 {
00392 mpCellVolumesFile = output_file_handler.OpenOutputFile("cellareas.dat");
00393 }
00394 }
00395
00396 template<unsigned DIM>
00397 void MeshBasedCellPopulation<DIM>::CloseOutputFiles()
00398 {
00399 AbstractCentreBasedCellPopulation<DIM>::CloseOutputFiles();
00400
00401 mpVizElementsFile->close();
00402
00403 if (mOutputVoronoiData)
00404 {
00405 mpVoronoiFile->close();
00406 }
00407 if (mOutputCellPopulationVolumes)
00408 {
00409 mpCellPopulationVolumesFile->close();
00410 }
00411 if (this->mOutputCellVolumes)
00412 {
00413 mpCellVolumesFile->close();
00414 }
00415 }
00416
00417 template<unsigned DIM>
00418 void MeshBasedCellPopulation<DIM>::WriteResultsToFiles()
00419 {
00420 AbstractCentreBasedCellPopulation<DIM>::WriteResultsToFiles();
00421
00422
00423 *mpVizElementsFile << SimulationTime::Instance()->GetTime() << "\t";
00424
00425 for (typename MutableMesh<DIM,DIM>::ElementIterator elem_iter = mrMesh.GetElementIteratorBegin();
00426 elem_iter != mrMesh.GetElementIteratorEnd();
00427 ++elem_iter)
00428 {
00429 bool element_contains_dead_cells_or_deleted_nodes = false;
00430
00431
00432 for (unsigned i=0; i<DIM+1; i++)
00433 {
00434 unsigned node_index = elem_iter->GetNodeGlobalIndex(i);
00435
00436 if (this->GetNode(node_index)->IsDeleted())
00437 {
00438 element_contains_dead_cells_or_deleted_nodes = true;
00439 break;
00440 }
00441 else if (this->mLocationCellMap[node_index])
00442 {
00443 if (this->mLocationCellMap[node_index]->IsDead())
00444 {
00445 element_contains_dead_cells_or_deleted_nodes = true;
00446 break;
00447 }
00448 }
00449 }
00450 if (!element_contains_dead_cells_or_deleted_nodes)
00451 {
00452 for (unsigned i=0; i<DIM+1; i++)
00453 {
00454 *mpVizElementsFile << elem_iter->GetNodeGlobalIndex(i) << " ";
00455 }
00456 }
00457 }
00458 *mpVizElementsFile << "\n";
00459
00460 if (mpVoronoiTessellation != NULL)
00461 {
00462 if (mOutputVoronoiData)
00463 {
00464 WriteVoronoiResultsToFile();
00465 }
00466 if (mOutputCellPopulationVolumes)
00467 {
00468 WriteCellPopulationVolumeResultsToFile();
00469 }
00470 if (this->mOutputCellVolumes)
00471 {
00472 WriteCellVolumeResultsToFile();
00473 }
00474 }
00475 }
00476
00477 template<unsigned DIM>
00478 void MeshBasedCellPopulation<DIM>::WriteVtkResultsToFile()
00479 {
00480 #ifdef CHASTE_VTK
00481
00482 std::stringstream time;
00483 time << SimulationTime::Instance()->GetTimeStepsElapsed();
00484
00485 unsigned num_points = GetNumNodes();
00486 if (!mWriteVtkAsPoints && (mpVoronoiTessellation != NULL))
00487 {
00488 num_points = mpVoronoiTessellation->GetNumElements();
00489 }
00490
00491 std::vector<double> cell_types(num_points);
00492 std::vector<double> cell_ancestors(num_points);
00493 std::vector<double> cell_mutation_states(num_points);
00494 std::vector<double> cell_ages(num_points);
00495 std::vector<double> cell_cycle_phases(num_points);
00496 std::vector<std::vector<double> > cellwise_data;
00497
00498 if (CellwiseData<DIM>::Instance()->IsSetUp())
00499 {
00500 CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00501 for (unsigned var=0; var<p_data->GetNumVariables(); var++)
00502 {
00503 std::vector<double> cellwise_data_var(num_points);
00504 cellwise_data.push_back(cellwise_data_var);
00505 }
00506 }
00507
00508 if (mWriteVtkAsPoints)
00509 {
00510 VtkMeshWriter<DIM, DIM> mesh_writer(this->mDirPath, "results_"+time.str(), false);
00511
00512
00513 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
00514 cell_iter != this->End();
00515 ++cell_iter)
00516 {
00517
00518 unsigned node_index = this->GetLocationIndexUsingCell(*cell_iter);
00519
00520
00521 AbstractCellCycleModel* p_model = cell_iter->GetCellCycleModel();
00522
00523 if (this->mOutputCellAncestors)
00524 {
00525 double ancestor_index = (cell_iter->GetAncestor() == UNSIGNED_UNSET) ? (-1.0) : (double)cell_iter->GetAncestor();
00526 cell_ancestors[node_index] = ancestor_index;
00527 }
00528 if (this->mOutputCellProliferativeTypes)
00529 {
00530 cell_types[node_index] = p_model->GetCellProliferativeType();
00531 }
00532 if (this->mOutputCellMutationStates)
00533 {
00534 cell_mutation_states[node_index] = cell_iter->GetMutationState()->GetColour();
00535 }
00536 if (this->mOutputCellAges)
00537 {
00538 cell_ages[node_index] = cell_iter->GetAge();
00539 }
00540 if (this->mOutputCellCyclePhases)
00541 {
00542 cell_cycle_phases[node_index] = p_model->GetCurrentCellCyclePhase();
00543 }
00544 if (CellwiseData<DIM>::Instance()->IsSetUp())
00545 {
00546 CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00547 for (unsigned var=0; var<p_data->GetNumVariables(); var++)
00548 {
00549 cellwise_data[var][node_index] = p_data->GetValue(*cell_iter, var);
00550 }
00551 }
00552 }
00553
00554 if (this->mOutputCellProliferativeTypes)
00555 {
00556 mesh_writer.AddPointData("Cell types", cell_types);
00557 }
00558 if (this->mOutputCellAncestors)
00559 {
00560 mesh_writer.AddPointData("Ancestors", cell_ancestors);
00561 }
00562 if (this->mOutputCellMutationStates)
00563 {
00564 mesh_writer.AddPointData("Mutation states", cell_mutation_states);
00565 }
00566 if (this->mOutputCellAges)
00567 {
00568 mesh_writer.AddPointData("Ages", cell_ages);
00569 }
00570 if (this->mOutputCellCyclePhases)
00571 {
00572 mesh_writer.AddPointData("Cycle phases", cell_cycle_phases);
00573 }
00574 if (CellwiseData<DIM>::Instance()->IsSetUp())
00575 {
00576 for (unsigned var=0; var<cellwise_data.size(); var++)
00577 {
00578 std::stringstream data_name;
00579 data_name << "Cellwise data " << var;
00580 std::vector<double> cellwise_data_var = cellwise_data[var];
00581 mesh_writer.AddPointData(data_name.str(), cellwise_data_var);
00582 }
00583 }
00584
00585 {
00586
00587 std::vector<Node<DIM>* > nodes;
00588 for (unsigned index=0; index<mrMesh.GetNumNodes(); index++)
00589 {
00590 Node<DIM>* p_node = mrMesh.GetNode(index);
00591 nodes.push_back(p_node);
00592 }
00593
00594 NodesOnlyMesh<DIM> mesh;
00595 mesh.ConstructNodesWithoutMesh(nodes);
00596 mesh_writer.WriteFilesUsingMesh(mesh);
00597 }
00598
00599 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
00600 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
00601 *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
00602 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
00603 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
00604 }
00605 else if (mpVoronoiTessellation != NULL)
00606 {
00607 VertexMeshWriter<DIM, DIM> mesh_writer(this->mDirPath, "results", false);
00608 std::vector<double> cell_volumes(num_points);
00609
00610
00611 for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = mpVoronoiTessellation->GetElementIteratorBegin();
00612 elem_iter != mpVoronoiTessellation->GetElementIteratorEnd();
00613 ++elem_iter)
00614 {
00615
00616 unsigned elem_index = elem_iter->GetIndex();
00617
00618
00619 unsigned node_index = mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
00620
00621
00622 assert(!this->IsGhostNode(node_index));
00623
00624
00625 CellPtr p_cell = this->mLocationCellMap[node_index];
00626
00627
00628 AbstractCellCycleModel* p_model = p_cell->GetCellCycleModel();
00629
00630 if (this->mOutputCellAncestors)
00631 {
00632 double ancestor_index = (p_cell->GetAncestor() == UNSIGNED_UNSET) ? (-1.0) : (double)p_cell->GetAncestor();
00633 cell_ancestors[elem_index] = ancestor_index;
00634 }
00635 if (this->mOutputCellProliferativeTypes)
00636 {
00637 cell_types[elem_index] = p_model->GetCellProliferativeType();
00638 }
00639 if (this->mOutputCellMutationStates)
00640 {
00641 cell_mutation_states[elem_index] = p_cell->GetMutationState()->GetColour();
00642 }
00643 if (this->mOutputCellAges)
00644 {
00645 cell_ages[elem_index] = p_cell->GetAge();
00646 }
00647 if (this->mOutputCellCyclePhases)
00648 {
00649 cell_cycle_phases[elem_index] = p_model->GetCurrentCellCyclePhase();
00650 }
00651 if (this->mOutputCellVolumes)
00652 {
00653 cell_volumes[elem_index] = mpVoronoiTessellation->GetVolumeOfElement(elem_index);
00654 }
00655 if (CellwiseData<DIM>::Instance()->IsSetUp())
00656 {
00657 CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00658 for (unsigned var=0; var<p_data->GetNumVariables(); var++)
00659 {
00660 cellwise_data[var][elem_index] = p_data->GetValue(p_cell, var);
00661 }
00662 }
00663 }
00664
00665 if (this->mOutputCellProliferativeTypes)
00666 {
00667 mesh_writer.AddCellData("Cell types", cell_types);
00668 }
00669 if (this->mOutputCellAncestors)
00670 {
00671 mesh_writer.AddCellData("Ancestors", cell_ancestors);
00672 }
00673 if (this->mOutputCellMutationStates)
00674 {
00675 mesh_writer.AddCellData("Mutation states", cell_mutation_states);
00676 }
00677 if (this->mOutputCellAges)
00678 {
00679 mesh_writer.AddCellData("Ages", cell_ages);
00680 }
00681 if (this->mOutputCellCyclePhases)
00682 {
00683 mesh_writer.AddCellData("Cycle phases", cell_cycle_phases);
00684 }
00685 if (this->mOutputCellVolumes)
00686 {
00687 mesh_writer.AddCellData("Cell volumes", cell_volumes);
00688 }
00689 if (CellwiseData<DIM>::Instance()->IsSetUp())
00690 {
00691 for (unsigned var=0; var<cellwise_data.size(); var++)
00692 {
00693 std::stringstream data_name;
00694 data_name << "Cellwise data " << var;
00695 std::vector<double> cellwise_data_var = cellwise_data[var];
00696 mesh_writer.AddCellData(data_name.str(), cellwise_data_var);
00697 }
00698 }
00699
00700 mesh_writer.WriteVtkUsingMesh(*mpVoronoiTessellation, time.str());
00701 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
00702 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
00703 *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
00704 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
00705 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
00706 }
00707 #endif //CHASTE_VTK
00708 }
00709
00710 template<unsigned DIM>
00711 void MeshBasedCellPopulation<DIM>::WriteVoronoiResultsToFile()
00712 {
00713 assert(DIM==2 || DIM==3);
00714
00715
00716 *mpVoronoiFile << SimulationTime::Instance()->GetTime() << " ";
00717
00718
00719 for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = mpVoronoiTessellation->GetElementIteratorBegin();
00720 elem_iter != mpVoronoiTessellation->GetElementIteratorEnd();
00721 ++elem_iter)
00722 {
00723
00724 unsigned elem_index = elem_iter->GetIndex();
00725
00726
00727 unsigned node_index = mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
00728
00729
00730 *mpVoronoiFile << node_index << " ";
00731 c_vector<double, DIM> node_location = this->GetNode(node_index)->rGetLocation();
00732 for (unsigned i=0; i<DIM; i++)
00733 {
00734 *mpVoronoiFile << node_location[i] << " ";
00735 }
00736
00737 double cell_volume = mpVoronoiTessellation->GetVolumeOfElement(elem_index);
00738 double cell_surface_area = mpVoronoiTessellation->GetSurfaceAreaOfElement(elem_index);
00739 *mpVoronoiFile << cell_volume << " " << cell_surface_area << " ";
00740 }
00741 *mpVoronoiFile << "\n";
00742 }
00743
00744 template<unsigned DIM>
00745 void MeshBasedCellPopulation<DIM>::WriteCellPopulationVolumeResultsToFile()
00746 {
00747 assert(DIM==2 || DIM==3);
00748
00749
00750 *mpCellPopulationVolumesFile << SimulationTime::Instance()->GetTime() << " ";
00751
00752
00753 double total_area = mrMesh.GetVolume();
00754 double apoptotic_area = 0.0;
00755
00756
00757 for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = mpVoronoiTessellation->GetElementIteratorBegin();
00758 elem_iter != mpVoronoiTessellation->GetElementIteratorEnd();
00759 ++elem_iter)
00760 {
00761
00762 unsigned elem_index = elem_iter->GetIndex();
00763
00764
00765 unsigned node_index = mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
00766
00767
00768 if (!this->IsGhostNode(node_index))
00769 {
00770
00771 CellPtr p_cell = this->mLocationCellMap[node_index];
00772
00773
00774 bool cell_is_apoptotic = p_cell->HasCellProperty<ApoptoticCellProperty>();
00775 if (cell_is_apoptotic)
00776 {
00777 double cell_volume = mpVoronoiTessellation->GetVolumeOfElement(elem_index);
00778 apoptotic_area += cell_volume;
00779 }
00780 }
00781 }
00782 *mpCellPopulationVolumesFile << total_area << " " << apoptotic_area << "\n";
00783 }
00784
00785 template<unsigned DIM>
00786 void MeshBasedCellPopulation<DIM>::WriteCellVolumeResultsToFile()
00787 {
00788 assert(DIM==2 || DIM==3);
00789
00790
00791 *mpCellVolumesFile << SimulationTime::Instance()->GetTime() << " ";
00792
00793
00794 for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = mpVoronoiTessellation->GetElementIteratorBegin();
00795 elem_iter != mpVoronoiTessellation->GetElementIteratorEnd();
00796 ++elem_iter)
00797 {
00798
00799 unsigned elem_index = elem_iter->GetIndex();
00800
00801
00802 unsigned node_index = mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
00803
00804
00805 if (!this->IsGhostNode(node_index))
00806 {
00807
00808 *mpCellVolumesFile << node_index << " ";
00809
00810
00811 CellPtr p_cell = this->mLocationCellMap[node_index];
00812
00813
00814 unsigned cell_index = p_cell->GetCellId();
00815 *mpCellVolumesFile << cell_index << " ";
00816
00817
00818 c_vector<double, DIM> node_location = this->GetNode(node_index)->rGetLocation();
00819 for (unsigned i=0; i<DIM; i++)
00820 {
00821 *mpCellVolumesFile << node_location[i] << " ";
00822 }
00823
00824
00825 double cell_volume = mpVoronoiTessellation->GetVolumeOfElement(elem_index);
00826 *mpCellVolumesFile << cell_volume << " ";
00827 }
00828 }
00829 *mpCellVolumesFile << "\n";
00830 }
00831
00832 template<unsigned DIM>
00833 void MeshBasedCellPopulation<DIM>::SetWriteVtkAsPoints(bool writeVtkAsPoints)
00834 {
00835 mWriteVtkAsPoints = writeVtkAsPoints;
00836 }
00837
00838 template<unsigned DIM>
00839 bool MeshBasedCellPopulation<DIM>::GetWriteVtkAsPoints()
00840 {
00841 return mWriteVtkAsPoints;
00842 }
00843
00845
00847
00848 template<unsigned DIM>
00849 Node<DIM>* MeshBasedCellPopulation<DIM>::SpringIterator::GetNodeA()
00850 {
00851 return mEdgeIter.GetNodeA();
00852 }
00853
00854 template<unsigned DIM>
00855 Node<DIM>* MeshBasedCellPopulation<DIM>::SpringIterator::GetNodeB()
00856 {
00857 return mEdgeIter.GetNodeB();
00858 }
00859
00860 template<unsigned DIM>
00861 CellPtr MeshBasedCellPopulation<DIM>::SpringIterator::GetCellA()
00862 {
00863 assert((*this) != mrCellPopulation.SpringsEnd());
00864 return mrCellPopulation.GetCellUsingLocationIndex(mEdgeIter.GetNodeA()->GetIndex());
00865 }
00866
00867 template<unsigned DIM>
00868 CellPtr MeshBasedCellPopulation<DIM>::SpringIterator::GetCellB()
00869 {
00870 assert((*this) != mrCellPopulation.SpringsEnd());
00871 return mrCellPopulation.GetCellUsingLocationIndex(mEdgeIter.GetNodeB()->GetIndex());
00872 }
00873
00874 template<unsigned DIM>
00875 bool MeshBasedCellPopulation<DIM>::SpringIterator::operator!=(const MeshBasedCellPopulation<DIM>::SpringIterator& rOther)
00876 {
00877 return (mEdgeIter != rOther.mEdgeIter);
00878 }
00879
00880 template<unsigned DIM>
00881 typename MeshBasedCellPopulation<DIM>::SpringIterator& MeshBasedCellPopulation<DIM>::SpringIterator::operator++()
00882 {
00883 bool edge_is_ghost = false;
00884
00885 do
00886 {
00887 ++mEdgeIter;
00888 if (*this != mrCellPopulation.SpringsEnd())
00889 {
00890 bool a_is_ghost = mrCellPopulation.IsGhostNode(mEdgeIter.GetNodeA()->GetIndex());
00891 bool b_is_ghost = mrCellPopulation.IsGhostNode(mEdgeIter.GetNodeB()->GetIndex());
00892
00893 edge_is_ghost = (a_is_ghost || b_is_ghost);
00894 }
00895 }
00896 while (*this!=mrCellPopulation.SpringsEnd() && edge_is_ghost);
00897
00898 return (*this);
00899 }
00900
00901 template<unsigned DIM>
00902 MeshBasedCellPopulation<DIM>::SpringIterator::SpringIterator(
00903 MeshBasedCellPopulation<DIM>& rCellPopulation,
00904 typename MutableMesh<DIM,DIM>::EdgeIterator edgeIter)
00905 : mrCellPopulation(rCellPopulation),
00906 mEdgeIter(edgeIter)
00907 {
00908 if (mEdgeIter!=mrCellPopulation.mrMesh.EdgesEnd())
00909 {
00910 bool a_is_ghost = mrCellPopulation.IsGhostNode(mEdgeIter.GetNodeA()->GetIndex());
00911 bool b_is_ghost = mrCellPopulation.IsGhostNode(mEdgeIter.GetNodeB()->GetIndex());
00912
00913 if (a_is_ghost || b_is_ghost)
00914 {
00915 ++(*this);
00916 }
00917 }
00918 }
00919
00920 template<unsigned DIM>
00921 typename MeshBasedCellPopulation<DIM>::SpringIterator MeshBasedCellPopulation<DIM>::SpringsBegin()
00922 {
00923 return SpringIterator(*this, mrMesh.EdgesBegin());
00924 }
00925
00926 template<unsigned DIM>
00927 typename MeshBasedCellPopulation<DIM>::SpringIterator MeshBasedCellPopulation<DIM>::SpringsEnd()
00928 {
00929 return SpringIterator(*this, mrMesh.EdgesEnd());
00930 }
00931
00935 template<>
00936 void MeshBasedCellPopulation<2>::CreateVoronoiTessellation()
00937 {
00938 delete mpVoronoiTessellation;
00939
00940
00941 bool is_mesh_periodic = false;
00942 if (dynamic_cast<Cylindrical2dMesh*>(&mrMesh))
00943 {
00944 is_mesh_periodic = true;
00945 }
00946
00947 mpVoronoiTessellation = new VertexMesh<2, 2>(mrMesh, is_mesh_periodic);
00948 }
00949
00955 template<>
00956 void MeshBasedCellPopulation<3>::CreateVoronoiTessellation()
00957 {
00958 delete mpVoronoiTessellation;
00959 mpVoronoiTessellation = new VertexMesh<3, 3>(mrMesh);
00960 }
00961
00966 template<>
00967 void MeshBasedCellPopulation<1>::CreateVoronoiTessellation()
00968 {
00969
00970 NEVER_REACHED;
00971 }
00972
00973 template<unsigned DIM>
00974 VertexMesh<DIM, DIM>* MeshBasedCellPopulation<DIM>::GetVoronoiTessellation()
00975 {
00976 assert(mpVoronoiTessellation!=NULL);
00977 return mpVoronoiTessellation;
00978 }
00979
00980 template<unsigned DIM>
00981 double MeshBasedCellPopulation<DIM>::GetVolumeOfVoronoiElement(unsigned index)
00982 {
00983 unsigned element_index = mpVoronoiTessellation->GetVoronoiElementIndexCorrespondingToDelaunayNodeIndex(index);
00984 double volume = mpVoronoiTessellation->GetVolumeOfElement(element_index);
00985 return volume;
00986 }
00987
00988 template<unsigned DIM>
00989 double MeshBasedCellPopulation<DIM>::GetSurfaceAreaOfVoronoiElement(unsigned index)
00990 {
00991 unsigned element_index = mpVoronoiTessellation->GetVoronoiElementIndexCorrespondingToDelaunayNodeIndex(index);
00992 double surface_area = mpVoronoiTessellation->GetSurfaceAreaOfElement(element_index);
00993 return surface_area;
00994 }
00995
00996 template<unsigned DIM>
00997 double MeshBasedCellPopulation<DIM>::GetVoronoiEdgeLength(unsigned index1, unsigned index2)
00998 {
00999 unsigned element_index1 = mpVoronoiTessellation->GetVoronoiElementIndexCorrespondingToDelaunayNodeIndex(index1);
01000 unsigned element_index2 = mpVoronoiTessellation->GetVoronoiElementIndexCorrespondingToDelaunayNodeIndex(index2);
01001 double edge_length = mpVoronoiTessellation->GetEdgeLength(element_index1, element_index2);
01002 return edge_length;
01003 }
01004
01005 template<unsigned DIM>
01006 void MeshBasedCellPopulation<DIM>::CheckCellPointers()
01007 {
01008 bool res = true;
01009 for (std::list<CellPtr>::iterator it=this->mCells.begin();
01010 it!=this->mCells.end();
01011 ++it)
01012 {
01013 CellPtr p_cell = *it;
01014 assert(p_cell);
01015 AbstractCellCycleModel* p_model = p_cell->GetCellCycleModel();
01016 assert(p_model);
01017
01018
01019 unsigned node_index = this->mCellLocationMap[p_cell.get()];
01020 std::cout << "Cell at node " << node_index << " addr " << p_cell << std::endl << std::flush;
01021 CellPtr p_cell_in_cell_population = this->GetCellUsingLocationIndex(node_index);
01022 #define COVERAGE_IGNORE //Debugging code. Shouldn't fail under normal conditions
01023 if (p_cell_in_cell_population != p_cell)
01024 {
01025 std::cout << " Mismatch with cell population" << std::endl << std::flush;
01026 res = false;
01027 }
01028
01029
01030 if (p_model->GetCell() != p_cell)
01031 {
01032 std::cout << " Mismatch with cycle model" << std::endl << std::flush;
01033 res = false;
01034 }
01035 }
01036 assert(res);
01037 #undef COVERAGE_IGNORE
01038
01039 res = true;
01040 for (std::set<std::pair<CellPtr,CellPtr> >::iterator it1 = mMarkedSprings.begin();
01041 it1 != mMarkedSprings.end();
01042 ++it1)
01043 {
01044 const std::pair<CellPtr,CellPtr>& r_pair = *it1;
01045
01046 for (unsigned i=0; i<2; i++)
01047 {
01048 CellPtr p_cell = (i==0 ? r_pair.first : r_pair.second);
01049
01050 assert(p_cell);
01051 AbstractCellCycleModel* p_model = p_cell->GetCellCycleModel();
01052 assert(p_model);
01053 unsigned node_index = this->mCellLocationMap[p_cell.get()];
01054 std::cout << "Cell at node " << node_index << " addr " << p_cell << std::endl << std::flush;
01055
01056 #define COVERAGE_IGNORE //Debugging code. Shouldn't fail under normal conditions
01057
01058 if (p_cell->IsDead())
01059 {
01060 std::cout << " Cell is dead" << std::endl << std::flush;
01061 res = false;
01062 }
01063
01064
01065 CellPtr p_cell_in_cell_population = this->GetCellUsingLocationIndex(node_index);
01066 if (p_cell_in_cell_population != p_cell)
01067 {
01068 std::cout << " Mismatch with cell population" << std::endl << std::flush;
01069 res = false;
01070 }
01071
01072
01073 if (p_model->GetCell() != p_cell)
01074 {
01075 std::cout << " Mismatch with cycle model" << std::endl << std::flush;
01076 res = false;
01077 }
01078 }
01079 #undef COVERAGE_IGNORE
01080 }
01081 assert(res);
01082 }
01083
01084 template<unsigned DIM>
01085 std::pair<CellPtr,CellPtr> MeshBasedCellPopulation<DIM>::CreateCellPair(CellPtr pCell1, CellPtr pCell2)
01086 {
01087 assert(pCell1);
01088 assert(pCell2);
01089
01090 std::pair<CellPtr,CellPtr> cell_pair;
01091
01092 if (pCell1->GetCellId() < pCell2->GetCellId())
01093 {
01094 cell_pair.first = pCell1;
01095 cell_pair.second = pCell2;
01096 }
01097 else
01098 {
01099 cell_pair.first = pCell2;
01100 cell_pair.second = pCell1;
01101 }
01102 return cell_pair;
01103 }
01104
01105 template<unsigned DIM>
01106 bool MeshBasedCellPopulation<DIM>::IsMarkedSpring(const std::pair<CellPtr,CellPtr>& rCellPair)
01107 {
01108
01109 assert(rCellPair.first->GetCellId() < rCellPair.second->GetCellId());
01110
01111 return mMarkedSprings.find(rCellPair) != mMarkedSprings.end();
01112 }
01113
01114 template<unsigned DIM>
01115 void MeshBasedCellPopulation<DIM>::MarkSpring(std::pair<CellPtr,CellPtr>& rCellPair)
01116 {
01117
01118 assert(rCellPair.first->GetCellId() < rCellPair.second->GetCellId());
01119
01120 mMarkedSprings.insert(rCellPair);
01121 }
01122
01123 template<unsigned DIM>
01124 void MeshBasedCellPopulation<DIM>::UnmarkSpring(std::pair<CellPtr,CellPtr>& rCellPair)
01125 {
01126
01127 assert(rCellPair.first->GetCellId() < rCellPair.second->GetCellId());
01128
01129 mMarkedSprings.erase(rCellPair);
01130 }
01131
01132 template<unsigned DIM>
01133 double MeshBasedCellPopulation<DIM>::GetAreaBasedDampingConstantParameter()
01134 {
01135 return mAreaBasedDampingConstantParameter;
01136 }
01137
01138 template<unsigned DIM>
01139 void MeshBasedCellPopulation<DIM>::SetAreaBasedDampingConstantParameter(double areaBasedDampingConstantParameter)
01140 {
01141 assert(areaBasedDampingConstantParameter >= 0.0);
01142 mAreaBasedDampingConstantParameter = areaBasedDampingConstantParameter;
01143 }
01144
01145 template<unsigned DIM>
01146 bool MeshBasedCellPopulation<DIM>::GetOutputVoronoiData()
01147 {
01148 return mOutputVoronoiData;
01149 }
01150
01151 template<unsigned DIM>
01152 void MeshBasedCellPopulation<DIM>::SetOutputVoronoiData(bool outputVoronoiData)
01153 {
01154 mOutputVoronoiData = outputVoronoiData;
01155 }
01156
01157 template<unsigned DIM>
01158 bool MeshBasedCellPopulation<DIM>::GetOutputCellPopulationVolumes()
01159 {
01160 return mOutputCellPopulationVolumes;
01161 }
01162
01163 template<unsigned DIM>
01164 void MeshBasedCellPopulation<DIM>::SetOutputCellPopulationVolumes(bool outputCellPopulationVolumes)
01165 {
01166 mOutputCellPopulationVolumes = outputCellPopulationVolumes;
01167 }
01168
01169 template<unsigned DIM>
01170 void MeshBasedCellPopulation<DIM>::OutputCellPopulationParameters(out_stream& rParamsFile)
01171 {
01172 *rParamsFile << "\t\t<UseAreaBasedDampingConstant>" << mUseAreaBasedDampingConstant << "</UseAreaBasedDampingConstant> \n";
01173 *rParamsFile << "\t\t<AreaBasedDampingConstantParameter>" << mAreaBasedDampingConstantParameter << "</AreaBasedDampingConstantParameter> \n";
01174 *rParamsFile << "\t\t<OutputVoronoiData>" << mOutputVoronoiData << "</OutputVoronoiData> \n";
01175 *rParamsFile << "\t\t<OutputCellPopulationVolumes>" << mOutputCellPopulationVolumes << "</OutputCellPopulationVolumes> \n";
01176
01177
01178 AbstractCentreBasedCellPopulation<DIM>::OutputCellPopulationParameters(rParamsFile);
01179 }
01180
01181 template<unsigned DIM>
01182 double MeshBasedCellPopulation<DIM>::GetWidth(const unsigned& rDimension)
01183 {
01184
01185 double width = mrMesh.GetWidth(rDimension);
01186 return width;
01187 }
01188
01190
01192
01193 template class MeshBasedCellPopulation<1>;
01194 template class MeshBasedCellPopulation<2>;
01195 template class MeshBasedCellPopulation<3>;
01196
01197
01198 #include "SerializationExportWrapperForCpp.hpp"
01199 EXPORT_TEMPLATE_CLASS_SAME_DIMS(MeshBasedCellPopulation)