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 #include "NodeBasedCellPopulation.hpp"
00029 #include "CellwiseData.hpp"
00030 #include "VtkMeshWriter.hpp"
00031
00032 template<unsigned DIM>
00033 NodeBasedCellPopulation<DIM>::NodeBasedCellPopulation(NodesOnlyMesh<DIM>& rMesh,
00034 std::vector<CellPtr>& rCells,
00035 const std::vector<unsigned> locationIndices,
00036 bool deleteMesh)
00037 : AbstractCentreBasedCellPopulation<DIM>(rCells, locationIndices),
00038 mrMesh(rMesh),
00039 mpBoxCollection(NULL),
00040 mDeleteMesh(deleteMesh),
00041 mMechanicsCutOffLength(DBL_MAX)
00042 {
00043 Validate();
00044 }
00045
00046 template<unsigned DIM>
00047 NodeBasedCellPopulation<DIM>::NodeBasedCellPopulation(NodesOnlyMesh<DIM>& rMesh)
00048 : AbstractCentreBasedCellPopulation<DIM>(),
00049 mrMesh(rMesh),
00050 mpBoxCollection(NULL),
00051 mDeleteMesh(true),
00052 mMechanicsCutOffLength(DBL_MAX)
00053 {
00054
00055 }
00056
00057 template<unsigned DIM>
00058 NodeBasedCellPopulation<DIM>::~NodeBasedCellPopulation()
00059 {
00060 Clear();
00061 if (mDeleteMesh)
00062 {
00063 delete &mrMesh;
00064 }
00065 }
00066
00067 template<unsigned DIM>
00068 NodesOnlyMesh<DIM>& NodeBasedCellPopulation<DIM>::rGetMesh()
00069 {
00070 return mrMesh;
00071 }
00072
00073 template<unsigned DIM>
00074 const NodesOnlyMesh<DIM>& NodeBasedCellPopulation<DIM>::rGetMesh() const
00075 {
00076 return mrMesh;
00077 }
00078
00079 template<unsigned DIM>
00080 void NodeBasedCellPopulation<DIM>::Clear()
00081 {
00082 delete mpBoxCollection;
00083 mpBoxCollection = NULL;
00084 mNodePairs.clear();
00085 }
00086
00087 template<unsigned DIM>
00088 void NodeBasedCellPopulation<DIM>::Validate()
00089 {
00090 std::vector<bool> validated_node(GetNumNodes());
00091 for (unsigned i=0; i<validated_node.size(); i++)
00092 {
00093 validated_node[i] = false;
00094 }
00095
00096 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
00097 {
00098 unsigned node_index = this->mCellLocationMap[(*cell_iter).get()];
00099 validated_node[node_index] = true;
00100 }
00101
00102 for (unsigned i=0; i<validated_node.size(); i++)
00103 {
00104 if (!validated_node[i])
00105 {
00106 std::stringstream ss;
00107 ss << "Node " << i << " does not appear to have a cell associated with it";
00108 EXCEPTION(ss.str());
00109 }
00110 }
00111 }
00112
00113 template<unsigned DIM>
00114 void NodeBasedCellPopulation<DIM>::SplitUpIntoBoxes(double cutOffLength, c_vector<double, 2*DIM> domainSize)
00115 {
00116 mpBoxCollection = new BoxCollection<DIM>(cutOffLength, domainSize);
00117 mpBoxCollection->SetupLocalBoxesHalfOnly();
00118
00119 for (unsigned i=0; i<mrMesh.GetNumNodes(); i++)
00120 {
00121 unsigned box_index = mpBoxCollection->CalculateContainingBox(this->GetNode(i));
00122 mpBoxCollection->rGetBox(box_index).AddNode(this->GetNode(i));
00123 }
00124 }
00125
00126 template<unsigned DIM>
00127 void NodeBasedCellPopulation<DIM>::FindMaxAndMin()
00128 {
00129 c_vector<double, DIM> min_posn;
00130 c_vector<double, DIM> max_posn;
00131 for (unsigned i=0; i<DIM; i++)
00132 {
00133 min_posn(i) = DBL_MAX;
00134 max_posn(i) = -DBL_MAX;
00135 }
00136
00137 for (unsigned i=0; i<mrMesh.GetNumNodes(); i++)
00138 {
00139 c_vector<double, DIM> posn = this->GetNode(i)->rGetLocation();
00140
00141 for (unsigned j=0; j<DIM; j++)
00142 {
00143 if (posn(j) > max_posn(j))
00144 {
00145 max_posn(j) = posn(j);
00146 }
00147 if (posn(j) < min_posn(j))
00148 {
00149 min_posn(j) = posn(j);
00150 }
00151 }
00152 }
00153
00154 for (unsigned i=0; i<DIM; i++)
00155 {
00156 assert(min_posn(i) != DBL_MAX);
00157 mMinSpatialPositions(i) = min_posn(i);
00158
00159 assert(max_posn(i) != -DBL_MAX);
00160 mMaxSpatialPositions(i) = max_posn(i);
00161 }
00162 }
00163
00164 template<unsigned DIM>
00165 Node<DIM>* NodeBasedCellPopulation<DIM>::GetNode(unsigned index)
00166 {
00167 return mrMesh.GetNode(index);
00168 }
00169
00170 template<unsigned DIM>
00171 void NodeBasedCellPopulation<DIM>::SetNode(unsigned nodeIndex, ChastePoint<DIM>& rNewLocation)
00172 {
00173 mrMesh.GetNode(nodeIndex)->SetPoint(rNewLocation);
00174 }
00175
00176 template<unsigned DIM>
00177 void NodeBasedCellPopulation<DIM>::Update(bool hasHadBirthsOrDeaths)
00178 {
00179 NodeMap map(mrMesh.GetNumAllNodes());
00180 mrMesh.ReMesh(map);
00181
00182 if (!map.IsIdentityMap())
00183 {
00184
00185 std::map<Cell*, unsigned> old_map = this->mCellLocationMap;
00186
00187
00188 this->mLocationCellMap.clear();
00189 this->mCellLocationMap.clear();
00190
00191 for (std::list<CellPtr>::iterator it = this->mCells.begin();
00192 it != this->mCells.end();
00193 ++it)
00194 {
00195 unsigned old_node_index = old_map[(*it).get()];
00196
00197
00198 assert(!map.IsDeleted(old_node_index));
00199
00200 unsigned new_node_index = map.GetNewIndex(old_node_index);
00201 this->mLocationCellMap[new_node_index] = *it;
00202 this->mCellLocationMap[(*it).get()] = new_node_index;
00203 }
00204
00205 this->Validate();
00206 }
00207
00208 mrMesh.SetMeshHasChangedSinceLoading();
00209
00210 if (mpBoxCollection != NULL)
00211 {
00212 delete mpBoxCollection;
00213 }
00214
00215 FindMaxAndMin();
00216
00217
00218 c_vector<double, 2*DIM> domain_size;
00219 for (unsigned i=0; i<DIM; i++)
00220 {
00221 domain_size(2*i) = mMinSpatialPositions(i);
00222 domain_size(2*i+1) = mMaxSpatialPositions(i);
00223 }
00224
00225 if (mMechanicsCutOffLength == DBL_MAX)
00226 {
00227 std::string error = std::string("NodeBasedCellPopulation cannot create boxes if the cut-off length has not been set - ")
00228 + std::string("Call SetMechanicsCutOffLength on the CellPopulation ensuring it is larger than GetCutOffLength() on the force law");
00229 EXCEPTION(error);
00230 }
00231
00232
00233
00234 SplitUpIntoBoxes(mMechanicsCutOffLength, domain_size);
00235
00237 std::vector<Node<DIM>*> nodes;
00238 for (unsigned index=0; index<mrMesh.GetNumNodes(); index++)
00239 {
00240 Node<DIM>* p_node = mrMesh.GetNode(index);
00241 nodes.push_back(p_node);
00242 }
00243 mpBoxCollection->CalculateNodePairs(nodes, mNodePairs);
00244 }
00245
00246 template<unsigned DIM>
00247 unsigned NodeBasedCellPopulation<DIM>::RemoveDeadCells()
00248 {
00249 unsigned num_removed = 0;
00250 for (std::list<CellPtr>::iterator it = this->mCells.begin();
00251 it != this->mCells.end();
00252 ++it)
00253 {
00254 if ((*it)->IsDead())
00255 {
00256
00257 num_removed++;
00258 mrMesh.DeleteNodePriorToReMesh(this->mCellLocationMap[(*it).get()]);
00259
00260
00261 unsigned location_index_of_removed_node = this->mCellLocationMap[(*it).get()];
00262 this->mCellLocationMap.erase((*it).get());
00263 this->mLocationCellMap.erase(location_index_of_removed_node);
00264
00265
00266 it = this->mCells.erase(it);
00267 --it;
00268 }
00269 }
00270
00271 return num_removed;
00272 }
00273
00274 template<unsigned DIM>
00275 unsigned NodeBasedCellPopulation<DIM>::AddNode(Node<DIM>* pNewNode)
00276 {
00277 return mrMesh.AddNode(pNewNode);
00278 }
00279
00280 template<unsigned DIM>
00281 unsigned NodeBasedCellPopulation<DIM>::GetNumNodes()
00282 {
00283 return mrMesh.GetNumAllNodes();
00284 }
00285
00286 template<unsigned DIM>
00287 BoxCollection<DIM>* NodeBasedCellPopulation<DIM>::GetBoxCollection()
00288 {
00289 return mpBoxCollection;
00290 }
00291
00292 template<unsigned DIM>
00293 std::set< std::pair<Node<DIM>*, Node<DIM>* > >& NodeBasedCellPopulation<DIM>::rGetNodePairs()
00294 {
00295 if (mNodePairs.empty())
00296 {
00297 EXCEPTION("No node pairs set up, rGetNodePairs probably called before Update");
00298 }
00299 return mNodePairs;
00300 }
00301
00302 template<unsigned DIM>
00303 void NodeBasedCellPopulation<DIM>::OutputCellPopulationParameters(out_stream& rParamsFile)
00304 {
00305 *rParamsFile << "\t\t<MechanicsCutOffLength>" << mMechanicsCutOffLength << "</MechanicsCutOffLength>\n";
00306
00307
00308
00309
00310 AbstractCentreBasedCellPopulation<DIM>::OutputCellPopulationParameters(rParamsFile);
00311 }
00312
00313 template<unsigned DIM>
00314 void NodeBasedCellPopulation<DIM>::SetMechanicsCutOffLength(double mechanicsCutOffLength)
00315 {
00316 assert(mechanicsCutOffLength > 0.0);
00317 mMechanicsCutOffLength = mechanicsCutOffLength;
00318 }
00319
00320 template<unsigned DIM>
00321 double NodeBasedCellPopulation<DIM>::GetMechanicsCutOffLength()
00322 {
00323 return mMechanicsCutOffLength;
00324 }
00325
00326 template<unsigned DIM>
00327 double NodeBasedCellPopulation<DIM>::GetWidth(const unsigned& rDimension)
00328 {
00329
00330 FindMaxAndMin();
00331
00332
00333 double width = mMaxSpatialPositions(rDimension) - mMinSpatialPositions(rDimension);
00334
00335 return width;
00336 }
00337
00338 template<unsigned DIM>
00339 void NodeBasedCellPopulation<DIM>::WriteVtkResultsToFile()
00340 {
00341 #ifdef CHASTE_VTK
00342 std::stringstream time;
00343 time << SimulationTime::Instance()->GetTimeStepsElapsed();
00344 VtkMeshWriter<DIM, DIM> mesh_writer(this->mDirPath, "results_"+time.str(), false);
00345
00346 unsigned num_nodes = GetNumNodes();
00347 std::vector<double> cell_types(num_nodes);
00348 std::vector<double> cell_ancestors(num_nodes);
00349 std::vector<double> cell_mutation_states(num_nodes);
00350 std::vector<double> cell_ages(num_nodes);
00351 std::vector<double> cell_cycle_phases(num_nodes);
00352 std::vector<std::vector<double> > cellwise_data;
00353
00354 if (CellwiseData<DIM>::Instance()->IsSetUp())
00355 {
00356 CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00357 unsigned num_variables = p_data->GetNumVariables();
00358 for (unsigned var=0; var<num_variables; var++)
00359 {
00360 std::vector<double> cellwise_data_var(num_nodes);
00361 cellwise_data.push_back(cellwise_data_var);
00362 }
00363 }
00364
00365
00366 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
00367 cell_iter != this->End();
00368 ++cell_iter)
00369 {
00370
00371 unsigned node_index = this->GetLocationIndexUsingCell(*cell_iter);
00372
00373 if (this->mOutputCellAncestors)
00374 {
00375 double ancestor_index = (cell_iter->GetAncestor() == UNSIGNED_UNSET) ? (-1.0) : (double)cell_iter->GetAncestor();
00376 cell_ancestors[node_index] = ancestor_index;
00377 }
00378 if (this->mOutputCellProliferativeTypes)
00379 {
00380 double cell_type = cell_iter->GetCellCycleModel()->GetCellProliferativeType();
00381 cell_types[node_index] = cell_type;
00382 }
00383 if (this->mOutputCellMutationStates)
00384 {
00385 double mutation_state = cell_iter->GetMutationState()->GetColour();
00386 cell_mutation_states[node_index] = mutation_state;
00387 }
00388 if (this->mOutputCellAges)
00389 {
00390 double age = cell_iter->GetAge();
00391 cell_ages[node_index] = age;
00392 }
00393 if (this->mOutputCellCyclePhases)
00394 {
00395 double cycle_phase = cell_iter->GetCellCycleModel()->GetCurrentCellCyclePhase();
00396 cell_cycle_phases[node_index] = cycle_phase;
00397 }
00398 if (CellwiseData<DIM>::Instance()->IsSetUp())
00399 {
00400 CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00401 unsigned num_variables = p_data->GetNumVariables();
00402 for (unsigned var=0; var<num_variables; var++)
00403 {
00404 cellwise_data[var][node_index] = p_data->GetValue(*cell_iter, var);
00405 }
00406 }
00407 }
00408
00409 if (this->mOutputCellProliferativeTypes)
00410 {
00411 mesh_writer.AddPointData("Cell types", cell_types);
00412 }
00413 if (this->mOutputCellAncestors)
00414 {
00415 mesh_writer.AddPointData("Ancestors", cell_ancestors);
00416 }
00417 if (this->mOutputCellMutationStates)
00418 {
00419 mesh_writer.AddPointData("Mutation states", cell_mutation_states);
00420 }
00421 if (this->mOutputCellAges)
00422 {
00423 mesh_writer.AddPointData("Ages", cell_ages);
00424 }
00425 if (this->mOutputCellCyclePhases)
00426 {
00427 mesh_writer.AddPointData("Cycle phases", cell_cycle_phases);
00428 }
00429 if (CellwiseData<DIM>::Instance()->IsSetUp())
00430 {
00431 for (unsigned var=0; var<cellwise_data.size(); var++)
00432 {
00433 std::stringstream data_name;
00434 data_name << "Cellwise data " << var;
00435 std::vector<double> cellwise_data_var = cellwise_data[var];
00436 mesh_writer.AddPointData(data_name.str(), cellwise_data_var);
00437 }
00438 }
00439
00440 mesh_writer.WriteFilesUsingMesh(mrMesh);
00441
00442 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
00443 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
00444 *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
00445 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
00446 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
00447 #endif //CHASTE_VTK
00448 }
00449
00450 template<unsigned DIM>
00451 void NodeBasedCellPopulation<DIM>::SetOutputCellVolumes(bool outputCellVolumes)
00452 {
00453 if (outputCellVolumes)
00454 {
00455 EXCEPTION("This method currently not implemented for a NodeBasedCellPopulation");
00456 }
00457 }
00458
00460
00462
00463 template class NodeBasedCellPopulation<1>;
00464 template class NodeBasedCellPopulation<2>;
00465 template class NodeBasedCellPopulation<3>;
00466
00467
00468 #include "SerializationExportWrapperForCpp.hpp"
00469 EXPORT_TEMPLATE_CLASS_SAME_DIMS(NodeBasedCellPopulation)