Chaste  Release::2018.1
AbstractCellPopulation.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include <boost/bind.hpp>
37 
38 #include <algorithm>
39 #include <functional>
40 
41 #include "AbstractCellPopulation.hpp"
42 #include "AbstractPhaseBasedCellCycleModel.hpp"
43 #include "SmartPointers.hpp"
44 #include "CellAncestor.hpp"
45 #include "ApoptoticCellProperty.hpp"
46 
47 // Cell writers
48 #include "BoundaryNodeWriter.hpp"
49 #include "CellProliferativeTypesWriter.hpp"
50 
51 // Cell population writers
52 #include "CellMutationStatesCountWriter.hpp"
53 #include "CellProliferativePhasesCountWriter.hpp"
54 #include "CellProliferativeTypesCountWriter.hpp"
55 #include "NodeLocationWriter.hpp"
56 
57 // These #includes are needed for SetDefaultCellMutationStateAndProliferativeTypeOrdering()
58 #include "WildTypeCellMutationState.hpp"
59 #include "ApcOneHitCellMutationState.hpp"
60 #include "ApcTwoHitCellMutationState.hpp"
61 #include "BetaCateninOneHitCellMutationState.hpp"
62 #include "DefaultCellProliferativeType.hpp"
63 #include "StemCellProliferativeType.hpp"
64 #include "TransitCellProliferativeType.hpp"
65 #include "DifferentiatedCellProliferativeType.hpp"
66 
67 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
69  std::vector<CellPtr>& rCells,
70  const std::vector<unsigned> locationIndices)
71  : mrMesh(rMesh),
72  mCells(rCells.begin(), rCells.end()),
73  mCentroid(zero_vector<double>(SPACE_DIM)),
74  mpCellPropertyRegistry(CellPropertyRegistry::Instance()->TakeOwnership()),
75  mOutputResultsForChasteVisualizer(true)
76 {
77  /*
78  * To avoid double-counting problems, clear the passed-in cells vector.
79  * We force a reallocation of memory so that subsequent usage of the
80  * vector is more likely to give an error.
81  */
82  std::vector<CellPtr>().swap(rCells);
83 
84  // There must be a one-one correspondence between cells and location indices
85  if (!locationIndices.empty())
86  {
87  if (mCells.size() != locationIndices.size())
88  {
89  EXCEPTION("There is not a one-one correspondence between cells and location indices");
90  }
91  }
92 
93  // Set up the map between location indices and cells
94  mLocationCellMap.clear();
95  mCellLocationMap.clear();
96 
97  std::list<CellPtr>::iterator it = mCells.begin();
98  for (unsigned i=0; it != mCells.end(); ++it, ++i)
99  {
100  // Give each cell a pointer to the property registry (we have taken ownership in this constructor)
101  (*it)->rGetCellPropertyCollection().SetCellPropertyRegistry(mpCellPropertyRegistry.get());
102  }
103 }
104 
105 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
107  : mrMesh(rMesh)
108 {
109 }
110 
111 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
113 {
114 }
115 
116 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
118 {
119  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter=this->Begin();
120  cell_iter!=this->End();
121  ++cell_iter)
122  {
123  cell_iter->InitialiseCellCycleModel();
124  cell_iter->InitialiseSrnModel();
125  }
126 }
127 
128 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
129 void AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::SetDataOnAllCells(const std::string& rDataName, double dataValue)
130 {
131  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter=this->Begin();
132  cell_iter!=this->End();
133  ++cell_iter)
134  {
135  cell_iter->GetCellData()->SetItem(rDataName, dataValue);
136  }
137 }
138 
139 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
141 {
142  return mrMesh;
143 }
144 
145 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
147 {
148  return mCells;
149 }
150 
151 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
153 {
154  unsigned counter = 0;
155  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter=this->Begin();
156  cell_iter!=this->End();
157  ++cell_iter)
158  {
159  counter++;
160  }
161  return counter;
162 }
163 
164 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
166 {
167  return mCells.size();
168 }
169 
170 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
172 {
173  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
174  {
175  MAKE_PTR_ARGS(CellAncestor, p_cell_ancestor, (mCellLocationMap[(*cell_iter).get()]));
176  cell_iter->SetAncestor(p_cell_ancestor);
177  }
178 }
179 
180 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
182 {
183  std::set<unsigned> remaining_ancestors;
184  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
185  {
186  remaining_ancestors.insert(cell_iter->GetAncestor());
187  }
188  return remaining_ancestors;
189 }
190 
191 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
193 {
194  std::vector<unsigned> mutation_state_count;
195  const std::vector<boost::shared_ptr<AbstractCellProperty> >& r_cell_properties
196  = mpCellPropertyRegistry->rGetAllCellProperties();
197 
198  // Calculate mutation states count
199  for (unsigned i=0; i<r_cell_properties.size(); i++)
200  {
201  if (r_cell_properties[i]->IsSubType<AbstractCellMutationState>())
202  {
203  mutation_state_count.push_back(r_cell_properties[i]->GetCellCount());
204  }
205  }
206 
207  // Reduce results onto all processes
209  {
210  // Make sure the vector on each process has the same size
211  unsigned local_size = mutation_state_count.size();
212  unsigned global_size;
213  MPI_Allreduce(&local_size, &global_size, 1, MPI_UNSIGNED, MPI_MAX, PetscTools::GetWorld());
214  assert(local_size == global_size);
215 
216  std::vector<unsigned> mutation_counts(global_size);
217  MPI_Allreduce(&mutation_state_count[0], &mutation_counts[0], mutation_counts.size(), MPI_UNSIGNED, MPI_SUM, PetscTools::GetWorld());
218 
219  mutation_state_count = mutation_counts;
220  }
221 
222  return mutation_state_count;
223 }
224 
225 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
227 {
228  std::vector<unsigned> proliferative_type_count;
229  const std::vector<boost::shared_ptr<AbstractCellProperty> >& r_cell_properties
230  = mpCellPropertyRegistry->rGetAllCellProperties();
231 
232  // Calculate proliferative types count
233  for (unsigned i=0; i<r_cell_properties.size(); i++)
234  {
235  if (r_cell_properties[i]->IsSubType<AbstractCellProliferativeType>())
236  {
237  proliferative_type_count.push_back(r_cell_properties[i]->GetCellCount());
238  }
239  }
240 
241  // Reduce results onto all processes
243  {
244  // Make sure the vector on each process has the same size
245  unsigned local_size = proliferative_type_count.size();
246  unsigned global_size;
247 
248  MPI_Allreduce(&local_size, &global_size, 1, MPI_UNSIGNED, MPI_MAX, PetscTools::GetWorld());
249  assert(local_size == global_size);
250 
251  std::vector<unsigned> total_types_counts(global_size);
252  MPI_Allreduce(&proliferative_type_count[0], &total_types_counts[0], total_types_counts.size(), MPI_UNSIGNED, MPI_SUM, PetscTools::GetWorld());
253 
254  proliferative_type_count = total_types_counts;
255  }
256 
257  return proliferative_type_count;
258 }
259 
260 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
262 {
263  std::vector<unsigned> cell_cycle_phase_count(5);
264  for (unsigned i=0; i<5; i++)
265  {
266  cell_cycle_phase_count[i] = 0;
267  }
268 
269  /*
270  * Note that in parallel with a poor partition a process could end up with zero cells
271  * in which case the calculation should be skipped since `this->Begin()` is not defined.
272  */
273  if (GetNumAllCells() > 0u)
274  {
275  if (dynamic_cast<AbstractPhaseBasedCellCycleModel*>((*(this->Begin()))->GetCellCycleModel()) == nullptr)
276  {
277  EXCEPTION("You are trying to record the cell cycle phase of cells with a non phase based cell cycle model.");
278  }
279 
280  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter = this->Begin();
281  cell_iter != this->End();
282  ++cell_iter)
283  {
284  switch (static_cast<AbstractPhaseBasedCellCycleModel*>((*cell_iter)->GetCellCycleModel())->GetCurrentCellCyclePhase())
285  {
286  case G_ZERO_PHASE:
287  cell_cycle_phase_count[0]++;
288  break;
289  case G_ONE_PHASE:
290  cell_cycle_phase_count[1]++;
291  break;
292  case S_PHASE:
293  cell_cycle_phase_count[2]++;
294  break;
295  case G_TWO_PHASE:
296  cell_cycle_phase_count[3]++;
297  break;
298  case M_PHASE:
299  cell_cycle_phase_count[4]++;
300  break;
301  default:
303  }
304  }
305  }
306 
307  // Reduce results onto all processes
309  {
310  std::vector<unsigned> phase_counts(cell_cycle_phase_count.size(), 0u);
311  MPI_Allreduce(&cell_cycle_phase_count[0], &phase_counts[0], phase_counts.size(), MPI_UNSIGNED, MPI_SUM, PetscTools::GetWorld());
312 
313  cell_cycle_phase_count = phase_counts;
314  }
315 
316  return cell_cycle_phase_count;
317 }
318 
319 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
321 {
322  // Get the set of pointers to cells corresponding to this location index
323  std::set<CellPtr> cells = mLocationCellMap[index];
324 
325  // If there is only one cell attached return the cell. Note currently only one cell per index.
326  if (cells.size() == 1)
327  {
328  return *(cells.begin());
329  }
330  if (cells.empty())
331  {
332  EXCEPTION("Location index input argument does not correspond to a Cell");
333  }
334  else
335  {
336  EXCEPTION("Multiple cells are attached to a single location index.");
337  }
338 }
339 
340 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
342 {
343  // Return the set of pointers to cells corresponding to this location index, note the set may be empty.
344  return mLocationCellMap[index];
345 }
346 
347 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
349 {
350  // Get the set of pointers to cells corresponding to this location index
351  std::set<CellPtr> cells = mLocationCellMap[index];
352 
353  // Return whether there is a cell attached to the location index
354  return !(cells.empty());
355 }
356 
357 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
359 {
360 }
361 
362 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
364 {
365  // Clear the maps
366  mLocationCellMap[index].clear();
367  mCellLocationMap.erase(pCell.get());
368 
369  // Replace with new cell
370  mLocationCellMap[index].insert(pCell);
371 
372  // Do other half of the map
373  mCellLocationMap[pCell.get()] = index;
374 }
375 
376 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
378 {
379  mLocationCellMap[index].insert(pCell);
380  mCellLocationMap[pCell.get()] = index;
381 }
382 
383 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
385 {
386  std::set<CellPtr>::iterator cell_iter = mLocationCellMap[index].find(pCell);
387 
388  if (cell_iter == mLocationCellMap[index].end())
389  {
390  EXCEPTION("Tried to remove a cell which is not attached to the given location index");
391  }
392  else
393  {
394  mLocationCellMap[index].erase(cell_iter);
395  mCellLocationMap.erase(pCell.get());
396  }
397 }
398 
399 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
400 void AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::MoveCellInLocationMap(CellPtr pCell, unsigned old_index, unsigned new_index)
401 {
402  // Remove the cell from its current location
403  RemoveCellUsingLocationIndex(old_index, pCell);
404 
405  // Add it to the new location
406  AddCellUsingLocationIndex(new_index, pCell);
407 }
408 
409 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
411 {
412  // Check the cell is in the map
413  assert(this->mCellLocationMap.find(pCell.get()) != this->mCellLocationMap.end());
414 
415  return mCellLocationMap[pCell.get()];
416 }
417 
418 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
420 {
421  return mpCellPropertyRegistry;
422 }
423 
424 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
426 {
427  boost::shared_ptr<CellPropertyRegistry> p_registry = GetCellPropertyRegistry();
428  if (!p_registry->HasOrderingBeenSpecified())
429  {
430  std::vector<boost::shared_ptr<AbstractCellProperty> > mutations_and_proliferative_types;
431  mutations_and_proliferative_types.push_back(p_registry->Get<WildTypeCellMutationState>());
432  mutations_and_proliferative_types.push_back(p_registry->Get<ApcOneHitCellMutationState>());
433  mutations_and_proliferative_types.push_back(p_registry->Get<ApcTwoHitCellMutationState>());
434  mutations_and_proliferative_types.push_back(p_registry->Get<BetaCateninOneHitCellMutationState>());
435  mutations_and_proliferative_types.push_back(p_registry->Get<StemCellProliferativeType>());
436  mutations_and_proliferative_types.push_back(p_registry->Get<TransitCellProliferativeType>());
437  mutations_and_proliferative_types.push_back(p_registry->Get<DifferentiatedCellProliferativeType>());
438 
439  // Parallel process with no cells won't have the default property, so add it in
440  mutations_and_proliferative_types.push_back(p_registry->Get<DefaultCellProliferativeType>());
441  p_registry->SpecifyOrdering(mutations_and_proliferative_types);
442  }
443 }
444 
445 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
447 {
448  mCentroid = zero_vector<double>(SPACE_DIM);
449  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter = this->Begin();
450  cell_iter != this->End();
451  ++cell_iter)
452  {
453  mCentroid += GetLocationOfCellCentre(*cell_iter);
454  }
455  mCentroid /= this->GetNumRealCells();
456 
457  return mCentroid;
458 }
459 
460 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
462 {
463 }
464 
465 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
467 {
468  typedef AbstractCellWriter<ELEMENT_DIM, SPACE_DIM> cell_writer_t;
469  BOOST_FOREACH(boost::shared_ptr<cell_writer_t> p_cell_writer, mCellWriters)
470  {
471  p_cell_writer->CloseFile();
472  }
473 
475  BOOST_FOREACH(boost::shared_ptr<pop_writer_t> p_pop_writer, mCellPopulationWriters)
476  {
477  p_pop_writer->CloseFile();
478  }
479 }
480 
481 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
483 {
484  typedef AbstractCellWriter<ELEMENT_DIM, SPACE_DIM> cell_writer_t;
485  BOOST_FOREACH(boost::shared_ptr<cell_writer_t> p_cell_writer, mCellWriters)
486  {
487  p_cell_writer->CloseFile();
488  }
489 
491  BOOST_FOREACH(boost::shared_ptr<pop_writer_t> p_pop_writer, mCellPopulationWriters)
492  {
493  p_pop_writer->CloseFile();
494  }
495 
496 #ifdef CHASTE_VTK
497  *mpVtkMetaFile << " </Collection>\n";
498  *mpVtkMetaFile << "</VTKFile>\n";
499  mpVtkMetaFile->close();
500 #endif //CHASTE_VTK
501 }
502 
503 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
505 {
506 #ifdef CHASTE_VTK
507  mpVtkMetaFile = rOutputFileHandler.OpenOutputFile("results.pvd");
508  *mpVtkMetaFile << "<?xml version=\"1.0\"?>\n";
509  *mpVtkMetaFile << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n";
510  *mpVtkMetaFile << " <Collection>\n";
511 #endif //CHASTE_VTK
512 
513  if (mOutputResultsForChasteVisualizer)
514  {
515  if (!HasWriter<NodeLocationWriter>())
516  {
517  AddPopulationWriter<NodeLocationWriter>();
518  }
519  if (!HasWriter<BoundaryNodeWriter>())
520  {
521  AddPopulationWriter<BoundaryNodeWriter>();
522  }
523  if (!HasWriter<CellProliferativeTypesWriter>())
524  {
525  AddCellWriter<CellProliferativeTypesWriter>();
526  }
527  }
528 
529  // Open output files for any cell writers
530  typedef AbstractCellWriter<ELEMENT_DIM, SPACE_DIM> cell_writer_t;
531  BOOST_FOREACH(boost::shared_ptr<cell_writer_t> p_cell_writer, mCellWriters)
532  {
533  p_cell_writer->OpenOutputFile(rOutputFileHandler);
534  }
535 
536  // Open output files and write headers for any population writers
538  BOOST_FOREACH(boost::shared_ptr<pop_writer_t> p_pop_writer, mCellPopulationWriters)
539  {
540  p_pop_writer->OpenOutputFile(rOutputFileHandler);
541  p_pop_writer->WriteHeader(this);
542  }
543 
544  // Open output files and write headers for any population count writers
546  BOOST_FOREACH(boost::shared_ptr<count_writer_t> p_count_writer, mCellPopulationCountWriters)
547  {
548  p_count_writer->OpenOutputFile(rOutputFileHandler);
549  p_count_writer->WriteHeader(this);
550  }
551 }
552 
553 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
555 {
556  typedef AbstractCellWriter<ELEMENT_DIM, SPACE_DIM> cell_writer_t;
558  BOOST_FOREACH(boost::shared_ptr<cell_writer_t> p_cell_writer, mCellWriters)
559  {
560  p_cell_writer->OpenOutputFileForAppend(rOutputFileHandler);
561  }
562  BOOST_FOREACH(boost::shared_ptr<pop_writer_t> p_pop_writer, mCellPopulationWriters)
563  {
564  p_pop_writer->OpenOutputFileForAppend(rOutputFileHandler);
565  }
566 }
567 
568 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
570 {
571  typedef AbstractCellWriter<ELEMENT_DIM, SPACE_DIM> cell_writer_t;
573  OutputFileHandler output_file_handler(rDirectory, false);
574 
575  if (!(mCellWriters.empty() && mCellPopulationWriters.empty() && mCellPopulationCountWriters.empty()))
576  {
577  // An ordering must be specified for cell mutation states and cell proliferative types
578  SetDefaultCellMutationStateAndProliferativeTypeOrdering();
579 
581  {
582  OpenRoundRobinWritersFilesForAppend(output_file_handler);
583 
584  // The master process writes time stamps
585  if (PetscTools::AmMaster())
586  {
587  BOOST_FOREACH(boost::shared_ptr<cell_writer_t> p_cell_writer, mCellWriters)
588  {
589  p_cell_writer->WriteTimeStamp();
590  }
591  BOOST_FOREACH(boost::shared_ptr<pop_writer_t> p_pop_writer, mCellPopulationWriters)
592  {
593  p_pop_writer->WriteTimeStamp();
594  }
595  }
596 
597  for (typename std::vector<boost::shared_ptr<AbstractCellPopulationWriter<ELEMENT_DIM, SPACE_DIM> > >::iterator pop_writer_iter = mCellPopulationWriters.begin();
598  pop_writer_iter != mCellPopulationWriters.end();
599  ++pop_writer_iter)
600  {
601  AcceptPopulationWriter(*pop_writer_iter);
602  }
603 
604  AcceptCellWritersAcrossPopulation();
605 
606  // The top-most process adds a newline
607  if (PetscTools::AmTopMost())
608  {
609  BOOST_FOREACH(boost::shared_ptr<cell_writer_t> p_cell_writer, mCellWriters)
610  {
611  p_cell_writer->WriteNewline();
612  }
613  BOOST_FOREACH(boost::shared_ptr<pop_writer_t> p_pop_writer, mCellPopulationWriters)
614  {
615  p_pop_writer->WriteNewline();
616  }
617  }
618  CloseRoundRobinWritersFiles();
619  }
621 
622  // Outside the round robin, deal with population count writers
624 
625  if (PetscTools::AmMaster())
626  {
627  // Open mCellPopulationCountWriters in append mode for writing, and write time stamps
628  BOOST_FOREACH(boost::shared_ptr<count_writer_t> p_count_writer, mCellPopulationCountWriters)
629  {
630  p_count_writer->OpenOutputFileForAppend(output_file_handler);
631  p_count_writer->WriteTimeStamp();
632  }
633  }
634  for (typename std::vector<boost::shared_ptr<AbstractCellPopulationCountWriter<ELEMENT_DIM, SPACE_DIM> > >::iterator count_writer_iter = mCellPopulationCountWriters.begin();
635  count_writer_iter != mCellPopulationCountWriters.end();
636  ++count_writer_iter)
637  {
638  AcceptPopulationCountWriter(*count_writer_iter);
639  }
640 
641  if (PetscTools::AmMaster())
642  {
643  // Add a newline and close any output files
644  BOOST_FOREACH(boost::shared_ptr<count_writer_t> p_count_writer, mCellPopulationCountWriters)
645  {
646  p_count_writer->WriteNewline();
647  p_count_writer->CloseFile();
648  }
649  }
650  }
651 
652  // VTK can only be written in 2 or 3 dimensions
653  if (SPACE_DIM > 1)
654  {
655  WriteVtkResultsToFile(rDirectory);
656  }
657 }
658 
659 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
661 {
662  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter = this->Begin();
663  cell_iter != this->End();
664  ++cell_iter)
665  {
666  for (typename std::vector<boost::shared_ptr<AbstractCellWriter<ELEMENT_DIM, SPACE_DIM> > >::iterator cell_writer_iter = mCellWriters.begin();
667  cell_writer_iter != mCellWriters.end();
668  ++cell_writer_iter)
669  {
670  AcceptCellWriter(*cell_writer_iter, *cell_iter);
671  }
672  }
673 }
674 
675 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
677 {
678  std::string cell_population_type = GetIdentifier();
679 
680  *rParamsFile << "\t<" << cell_population_type << ">\n";
681  OutputCellPopulationParameters(rParamsFile);
682  *rParamsFile << "\t</" << cell_population_type << ">\n";
683  *rParamsFile << "\n";
684  *rParamsFile << "\t<CellCycleModels>\n";
685 
692  std::set<std::string> unique_cell_cycle_models;
693  std::vector<CellPtr> first_cell_with_unique_CCM;
694  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter = this->Begin();
695  cell_iter != this->End();
696  ++cell_iter)
697  {
698  std::string identifier = cell_iter->GetCellCycleModel()->GetIdentifier();
699  if (unique_cell_cycle_models.count(identifier) == 0)
700  {
701  unique_cell_cycle_models.insert(identifier);
702  first_cell_with_unique_CCM.push_back((*cell_iter));
703  }
704  }
705 
706  // Loop over unique cell-cycle models
707  for (unsigned i=0; i<first_cell_with_unique_CCM.size(); i++)
708  {
709  // Output cell-cycle model details
710  first_cell_with_unique_CCM[i]->GetCellCycleModel()->OutputCellCycleModelInfo(rParamsFile);
711  }
712  *rParamsFile << "\t</CellCycleModels>\n";
713 
714  *rParamsFile << "\n";
715  *rParamsFile << "\t<SrnModels>\n";
716 
723  std::set<std::string> unique_srn_models;
724  std::vector<CellPtr> first_cell_with_unique_SRN;
725  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter = this->Begin();
726  cell_iter != this->End();
727  ++cell_iter)
728  {
729  std::string identifier = cell_iter->GetSrnModel()->GetIdentifier();
730  if (unique_srn_models.count(identifier) == 0)
731  {
732  unique_srn_models.insert(identifier);
733  first_cell_with_unique_SRN.push_back((*cell_iter));
734  }
735  }
736 
737  // Loop over unique SRN models
738  for (unsigned i=0; i<first_cell_with_unique_SRN.size(); i++)
739  {
740  // Output SRN model details
741  first_cell_with_unique_SRN[i]->GetSrnModel()->OutputSrnModelInfo(rParamsFile);
742  }
743 
744  *rParamsFile << "\t</SrnModels>\n";
745 }
746 
747 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
749 {
750  *rParamsFile << "\t\t<OutputResultsForChasteVisualizer>" << mOutputResultsForChasteVisualizer << "</OutputResultsForChasteVisualizer>\n";
751 }
752 
753 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
755 {
756 }
757 
758 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
760 {
761  return mOutputResultsForChasteVisualizer;
762 }
763 
764 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
766 {
767  mOutputResultsForChasteVisualizer = outputResultsForChasteVisualizer;
768 }
769 
770 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
772 {
773  return true;
774 }
775 
776 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
778 {
779  // Compute the centre of mass of the cell population
780  c_vector<double,SPACE_DIM> centre = GetCentroidOfCellPopulation();
781 
782  // Loop over cells and find the maximum distance from the centre of mass in each dimension
783  c_vector<double,SPACE_DIM> max_distance_from_centre = zero_vector<double>(SPACE_DIM);
784  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter = this->Begin();
785  cell_iter != this->End();
786  ++cell_iter)
787  {
788  c_vector<double,SPACE_DIM> cell_location = GetLocationOfCellCentre(*cell_iter);
789 
790  // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
791  c_vector<double,SPACE_DIM> displacement;
792  displacement = centre - cell_location;
793 
794  for (unsigned i=0; i<SPACE_DIM; i++)
795  {
796  if (displacement[i] > max_distance_from_centre[i])
797  {
798  max_distance_from_centre[i] = displacement[i];
799  }
800  }
801  }
802 
803  return max_distance_from_centre;
804 }
805 
806 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
807 std::pair<unsigned,unsigned> AbstractCellPopulation<ELEMENT_DIM,SPACE_DIM>::CreateOrderedPair(unsigned index1, unsigned index2)
808 {
809  assert(index1 != index2);
810 
811  std::pair<unsigned, unsigned> ordered_pair;
812  if (index1 < index2)
813  {
814  ordered_pair.first = index1;
815  ordered_pair.second = index2;
816  }
817  else
818  {
819  ordered_pair.first = index2;
820  ordered_pair.second = index1;
821  }
822  return ordered_pair;
823 }
824 
825 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
827 {
828  bool non_apoptotic_cell_present = false;
829 
830  if (IsCellAttachedToLocationIndex(pdeNodeIndex))
831  {
832  non_apoptotic_cell_present = !(GetCellUsingLocationIndex(pdeNodeIndex)->template HasCellProperty<ApoptoticCellProperty>());
833  }
834 
835  return non_apoptotic_cell_present;
836 }
837 
838 // Explicit instantiation
839 template class AbstractCellPopulation<1,1>;
840 template class AbstractCellPopulation<1,2>;
841 template class AbstractCellPopulation<2,2>;
842 template class AbstractCellPopulation<1,3>;
843 template class AbstractCellPopulation<2,3>;
844 template class AbstractCellPopulation<3,3>;
virtual void SimulationSetupHook(AbstractCellBasedSimulation< ELEMENT_DIM, SPACE_DIM > *pSimulation)
c_vector< double, SPACE_DIM > GetCentroidOfCellPopulation()
void OpenOutputFileForAppend(OutputFileHandler &rOutputFileHandler)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
void SetDefaultCellMutationStateAndProliferativeTypeOrdering()
unsigned GetLocationIndexUsingCell(CellPtr pCell)
void SetDataOnAllCells(const std::string &rDataName, double dataValue)
boost::shared_ptr< CellPropertyRegistry > mpCellPropertyRegistry
virtual void RemoveCellUsingLocationIndex(unsigned index, CellPtr pCell)
virtual void OutputCellPopulationParameters(out_stream &rParamsFile)=0
#define EXCEPTION(message)
Definition: Exception.hpp:143
static bool AmMaster()
Definition: PetscTools.cpp:120
static MPI_Comm GetWorld()
Definition: PetscTools.cpp:174
std::set< CellPtr > GetCellsUsingLocationIndex(unsigned index)
virtual bool IsCellAttachedToLocationIndex(unsigned index)
#define NEVER_REACHED
Definition: Exception.hpp:206
virtual void WriteResultsToFiles(const std::string &rDirectory)
void SetOutputResultsForChasteVisualizer(bool outputResultsForChasteVisualizer)
void MoveCellInLocationMap(CellPtr pCell, unsigned old_index, unsigned new_index)
std::list< CellPtr > & rGetCells()
std::vector< unsigned > GetCellCyclePhaseCount()
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
c_vector< double, SPACE_DIM > GetSizeOfCellPopulation()
virtual void AcceptCellWritersAcrossPopulation()
std::pair< unsigned, unsigned > CreateOrderedPair(unsigned index1, unsigned index2)
std::map< unsigned, std::set< CellPtr > > mLocationCellMap
void SetCellUsingLocationIndex(unsigned index, CellPtr pCell)
std::map< Cell *, unsigned > mCellLocationMap
static void EndRoundRobin()
Definition: PetscTools.cpp:160
boost::shared_ptr< CellPropertyRegistry > GetCellPropertyRegistry()
static void BeginRoundRobin()
Definition: PetscTools.cpp:150
void OpenRoundRobinWritersFilesForAppend(OutputFileHandler &rOutputFileHandler)
virtual bool IsPdeNodeAssociatedWithNonApoptoticCell(unsigned pdeNodeIndex)
virtual void OpenWritersFiles(OutputFileHandler &rOutputFileHandler)
static bool IsParallel()
Definition: PetscTools.cpp:97
void OutputCellPopulationInfo(out_stream &rParamsFile)
virtual bool IsRoomToDivide(CellPtr pCell)
AbstractMesh< ELEMENT_DIM, SPACE_DIM > & rGetMesh()
virtual void WriteDataToVisualizerSetupFile(out_stream &pVizSetupFile)
virtual void OpenOutputFile(OutputFileHandler &rOutputFileHandler)
static bool AmTopMost()
Definition: PetscTools.cpp:126
AbstractCellPopulation(AbstractMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
virtual void AddCellUsingLocationIndex(unsigned index, CellPtr pCell)
std::vector< unsigned > GetCellProliferativeTypeCount()
#define MAKE_PTR_ARGS(TYPE, NAME, ARGS)
std::vector< unsigned > GetCellMutationStateCount()
std::set< unsigned > GetCellAncestors()