Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
CaBasedCellPopulation.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include <boost/scoped_array.hpp>
37
38#include "CaBasedCellPopulation.hpp"
39#include "MutableMesh.hpp"
40#include "AbstractCaUpdateRule.hpp"
41#include "AbstractCaSwitchingUpdateRule.hpp"
42#include "RandomNumberGenerator.hpp"
43#include "CellLocationIndexWriter.hpp"
44#include "ExclusionCaBasedDivisionRule.hpp"
45#include "NodesOnlyMesh.hpp"
46#include "ApoptoticCellProperty.hpp"
47
48// Needed to convert mesh in order to write nodes to VTK (visualize as glyphs)
49#include "VtkMeshWriter.hpp"
50
51// LCOV_EXCL_START
52template<unsigned DIM>
57// LCOV_EXCL_STOP
58
59template<unsigned DIM>
61 std::vector<CellPtr>& rCells,
62 const std::vector<unsigned> locationIndices,
63 unsigned latticeCarryingCapacity,
64 bool deleteMesh,
65 bool validate)
66 : AbstractOnLatticeCellPopulation<DIM>(rMesh, rCells, locationIndices, deleteMesh),
67 mLatticeCarryingCapacity(latticeCarryingCapacity)
68{
69 mAvailableSpaces = std::vector<unsigned>(this->GetNumNodes(), latticeCarryingCapacity);
71
72 // This must always be true
73 assert(this->mCells.size() <= this->mrMesh.GetNumNodes()*latticeCarryingCapacity);
74
75 if (locationIndices.empty())
76 {
77 EXCEPTION("No location indices being passed. Specify where cells lie before creating the cell population.");
78 }
79 else
80 {
81 // Create a set of node indices corresponding to empty sites.
82 // Note iterating over mCells is OK as it has the same order as location indices at this point (its just coppied from rCells)
83 std::list<CellPtr>::iterator it = this->mCells.begin();
84 for (unsigned i=0; it != this->mCells.end(); ++it, ++i)
85 {
86 assert(i < locationIndices.size());
87 if (!IsSiteAvailable(locationIndices[i],*it))
88 {
89 EXCEPTION("One of the lattice sites has more cells than the carrying capacity. Check the initial cell locations.");
90 }
91 mAvailableSpaces[locationIndices[i]]--;
92 }
93 }
94 if (validate)
95 {
96 EXCEPTION("There is no validation for CaBasedCellPopulation.");
97 }
98}
99
100template<unsigned DIM>
105
106template<unsigned DIM>
108{
109 if (this->mDeleteMesh)
110 {
111 delete &this->mrMesh;
112 }
113}
114
115template<unsigned DIM>
117{
118 return mAvailableSpaces;
119}
120
121template<unsigned DIM>
122bool CaBasedCellPopulation<DIM>::IsSiteAvailable(unsigned index, CellPtr pCell)
123{
125 return (mAvailableSpaces[index] != 0);
126}
127
128template<unsigned DIM>
130{
131 return static_cast<PottsMesh<DIM>& >((this->mrMesh));
132}
133
134template<unsigned DIM>
136{
137 return static_cast<PottsMesh<DIM>& >((this->mrMesh));
138}
139
140template<unsigned DIM>
142{
143 std::vector<Node<DIM>*> temp_nodes;
144
145 // Create nodes at the centre of the cells
146 unsigned cell_index = 0;
147 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
148 cell_iter != this->End();
149 ++cell_iter)
150 {
151 temp_nodes.push_back(new Node<DIM>(cell_index, this->GetLocationOfCellCentre(*cell_iter)));
152 cell_index++;
153 }
154
155 return new MutableMesh<DIM,DIM>(temp_nodes);
156}
157
158template<unsigned DIM>
160{
161 return this->mrMesh.GetNode(index);
162}
163
164template<unsigned DIM>
166{
167 return this->mrMesh.GetNumNodes();
168}
169
170template<unsigned DIM>
172{
173 unsigned index = this->GetLocationIndexUsingCell(pCell);
174 std::set<unsigned> candidates = static_cast<PottsMesh<DIM>& >((this->mrMesh)).GetMooreNeighbouringNodeIndices(index);
175
176 std::set<unsigned> neighbour_indices;
177 for (std::set<unsigned>::iterator iter = candidates.begin();
178 iter != candidates.end();
179 ++iter)
180 {
181 if (!IsSiteAvailable(*iter, pCell))
182 {
183 neighbour_indices.insert(*iter);
184 }
185 }
186
187 return neighbour_indices;
188}
189
190template<unsigned DIM>
191c_vector<double, DIM> CaBasedCellPopulation<DIM>::GetLocationOfCellCentre(CellPtr pCell)
192{
193 return this->mrMesh.GetNode(this->GetLocationIndexUsingCell(pCell))->rGetLocation();
194}
195
196template<unsigned DIM>
198{
199 return this->mrMesh.GetNode(this->GetLocationIndexUsingCell(pCell));
200}
201
202template<unsigned DIM>
204{
205 if (!IsSiteAvailable(index, pCell))
206 {
207 EXCEPTION("No available spaces at location index " << index << ".");
208 }
209
210 mAvailableSpaces[index]--;
212}
213
214template<unsigned DIM>
216{
218
219 mAvailableSpaces[index]++;
220
221 assert(mAvailableSpaces[index] <= mLatticeCarryingCapacity);
222}
223
224template<unsigned DIM>
226{
227 return mpCaBasedDivisionRule->IsRoomToDivide(pCell, *this);
228}
229
230template<unsigned DIM>
231CellPtr CaBasedCellPopulation<DIM>::AddCell(CellPtr pNewCell, CellPtr pParentCell)
232{
233 unsigned daughter_node_index = mpCaBasedDivisionRule->CalculateDaughterNodeIndex(pNewCell,pParentCell,*this);
234
235 // Associate the new cell with the neighbouring node
236 this->mCells.push_back(pNewCell);
237
238 // Update location cell map
239 CellPtr p_created_cell = this->mCells.back();
240 AddCellUsingLocationIndex(daughter_node_index,p_created_cell);
241
242 return p_created_cell;
243}
244
245template<unsigned DIM>
247 unsigned targetNodeIndex,
248 CellPtr pCell)
249{
250 return 1.0;
251}
252
253template<unsigned DIM>
255{
256 unsigned num_removed = 0;
257
258 for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
259 cell_iter != this->mCells.end();
260 )
261 {
262 if ((*cell_iter)->IsDead())
263 {
264 // Get the location index corresponding to this cell
265 unsigned location_index = this->GetLocationIndexUsingCell(*cell_iter);
266
267 // Use this to remove the cell from the population
268 RemoveCellUsingLocationIndex(location_index, (*cell_iter));
269
270 // Erase cell and update counter
271 cell_iter = this->mCells.erase(cell_iter);
272 num_removed++;
273 }
274 else
275 {
276 ++cell_iter;
277 }
278 }
279 return num_removed;
280}
281
282template<unsigned DIM>
284{
285 /*
286 * Here we loop over the nodes and calculate the probability of moving
287 * and then select the node to move to.
288 */
289 if (!(this->mUpdateRuleCollection.empty()))
290 {
291 // Iterate over cells
293 for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
294 cell_iter != this->mCells.end();
295 ++cell_iter)
296 {
297 // Loop over neighbours and calculate probability of moving (make sure all probabilities are <1)
298 unsigned node_index = this->GetLocationIndexUsingCell(*cell_iter);
299
300 // Find a random available neighbouring node to overwrite current site
301 std::set<unsigned> neighbouring_node_indices = static_cast<PottsMesh<DIM>& >((this->mrMesh)).GetMooreNeighbouringNodeIndices(node_index);
302 std::vector<double> neighbouring_node_propensities;
303 std::vector<unsigned> neighbouring_node_indices_vector;
304
305 if (!neighbouring_node_indices.empty())
306 {
307 unsigned num_neighbours = neighbouring_node_indices.size();
308 double probability_of_not_moving = 1.0;
309
310 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
311 iter != neighbouring_node_indices.end();
312 ++iter)
313 {
314 double probability_of_moving = 0.0;
315
316 neighbouring_node_indices_vector.push_back(*iter);
317
318 if (IsSiteAvailable(*iter, *cell_iter))
319 {
320 // Iterating over the update rule
321 for (typename std::vector<boost::shared_ptr<AbstractUpdateRule<DIM> > >::iterator iter_rule = this->mUpdateRuleCollection.begin();
322 iter_rule != this->mUpdateRuleCollection.end();
323 ++iter_rule)
324 {
325 // This static cast is fine, since we assert the update rule must be a CA update rule in AddUpdateRule()
326 double p = (boost::static_pointer_cast<AbstractCaUpdateRule<DIM> >(*iter_rule))->EvaluateProbability(node_index, *iter, *this, dt, 1, *cell_iter);
327 probability_of_moving += p;
328 if (probability_of_moving < 0)
329 {
330 EXCEPTION("The probability of cellular movement is smaller than zero. In order to prevent it from happening you should change your time step and parameters");
331 }
332
333 if (probability_of_moving > 1)
334 {
335 EXCEPTION("The probability of the cellular movement is bigger than one. In order to prevent it from happening you should change your time step and parameters");
336 }
337 }
338
339 probability_of_not_moving -= probability_of_moving;
340 neighbouring_node_propensities.push_back(probability_of_moving);
341 }
342 else
343 {
344 neighbouring_node_propensities.push_back(0.0);
345 }
346 }
347 if (probability_of_not_moving < 0)
348 {
349 EXCEPTION("The probability of the cell not moving is smaller than zero. In order to prevent it from happening you should change your time step and parameters");
350 }
351
352 // Sample random number to specify which move to make
354 double random_number = p_gen->ranf();
355
356 double total_probability = 0.0;
357 for (unsigned counter=0; counter<num_neighbours; counter++)
358 {
359 total_probability += neighbouring_node_propensities[counter];
360 if (total_probability >= random_number)
361 {
362 // Move the cell to this neighbour location
363 unsigned chosen_neighbour_location_index = neighbouring_node_indices_vector[counter];
364 this->MoveCellInLocationMap((*cell_iter), node_index, chosen_neighbour_location_index);
365 break;
366 }
367 }
368 // If loop completes with total_probability < random_number then stay in the same location
369 }
370 else
371 {
372 // Each node in the mesh must have at least one neighbour
374 }
375 }
376 }
377
378 /*
379 * Here we loop over the nodes and select a neighbour to test if
380 * the cells (associated with the nodes) should swap locations
381 * or if a cell should move to an empty node
382 * Note this currently only works for latticeCarryingCapacity = 1
383 */
384 if (!(mSwitchingUpdateRuleCollection.empty()))
385 {
386 assert(mLatticeCarryingCapacity == 1);
387
389 unsigned num_nodes = this->mrMesh.GetNumNodes();
390
391 // Randomly permute mUpdateRuleCollection if specified
392 if (this->mIterateRandomlyOverUpdateRuleCollection)
393 {
394 // Randomly permute mUpdateRuleCollection
395 p_gen->Shuffle(mSwitchingUpdateRuleCollection);
396 }
397
398 for (unsigned i=0; i<num_nodes; i++)
399 {
400 unsigned node_index;
401
402 if (this->mUpdateNodesInRandomOrder)
403 {
404 node_index = p_gen->randMod(num_nodes);
405 }
406 else
407 {
408 // Loop over nodes in index order
409 node_index = i%num_nodes;
410 }
411
412 // Find a random available neighbouring node to switch cells with the current site
413 std::set<unsigned> neighbouring_node_indices = static_cast<PottsMesh<DIM>& >((this->mrMesh)).GetMooreNeighbouringNodeIndices(node_index);
414
415 unsigned neighbour_location_index;
416
417 if (!neighbouring_node_indices.empty())
418 {
419 unsigned num_neighbours = neighbouring_node_indices.size();
420 unsigned chosen_neighbour = p_gen->randMod(num_neighbours);
421
422 std::set<unsigned>::iterator neighbour_iter = neighbouring_node_indices.begin();
423 for (unsigned j=0; j<chosen_neighbour; j++)
424 {
425 neighbour_iter++;
426 }
427 neighbour_location_index = *neighbour_iter;
428
429 bool is_cell_on_node_index = mAvailableSpaces[node_index] == 0 ? true : false;
430 bool is_cell_on_neighbour_location_index = mAvailableSpaces[neighbour_location_index] == 0 ? true : false;
431
432 if (is_cell_on_node_index || is_cell_on_neighbour_location_index)
433 {
434 double probability_of_switch = 0.0;
435
436 // Now add contributions to the probability from each CA switching update rule
437 for (typename std::vector<boost::shared_ptr<AbstractUpdateRule<DIM> > >::iterator iter_rule = mSwitchingUpdateRuleCollection.begin();
438 iter_rule != mSwitchingUpdateRuleCollection.end();
439 ++iter_rule)
440 {
441 // This static cast is fine, since we assert the update rule must be a CA switching update rule in AddUpdateRule()
442 double p = (boost::static_pointer_cast<AbstractCaSwitchingUpdateRule<DIM> >(*iter_rule))->EvaluateSwitchingProbability(node_index, neighbour_location_index, *this, dt, 1);
443 probability_of_switch += p;
444 }
445
446 assert(probability_of_switch >= 0);
447 assert(probability_of_switch <= 1);
448
449 // Generate a uniform random number to do the random switch
450 double random_number = p_gen->ranf();
451
452 if (random_number < probability_of_switch)
453 {
454 if (is_cell_on_node_index && is_cell_on_neighbour_location_index)
455 {
456 // Swap the cells associated with the node and the neighbour node
457 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
458 CellPtr p_neighbour_cell = this->GetCellUsingLocationIndex(neighbour_location_index);
459
460 // Remove the cells from their current location
461 RemoveCellUsingLocationIndex(node_index, p_cell);
462 RemoveCellUsingLocationIndex(neighbour_location_index, p_neighbour_cell);
463
464 // Add cells to their new locations
465 AddCellUsingLocationIndex(node_index, p_neighbour_cell);
466 AddCellUsingLocationIndex(neighbour_location_index, p_cell);
467 }
468 else if (is_cell_on_node_index && !is_cell_on_neighbour_location_index)
469 {
470 // Move the cells associated with the node to the neighbour node
471 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
472 RemoveCellUsingLocationIndex(node_index, p_cell);
473 AddCellUsingLocationIndex(neighbour_location_index, p_cell);
474 }
475 else if (!is_cell_on_node_index && is_cell_on_neighbour_location_index)
476 {
477 // Move the cell associated with the neighbour node onto the node
478 CellPtr p_neighbour_cell = this->GetCellUsingLocationIndex(neighbour_location_index);
479 RemoveCellUsingLocationIndex(neighbour_location_index, p_neighbour_cell);
480 AddCellUsingLocationIndex(node_index, p_neighbour_cell);
481 }
482 else
483 {
485 }
486 }
487 }
488 }
489 }
490 }
491}
492
493template<unsigned DIM>
495{
496 return false;
497}
498
499template<unsigned DIM>
500void CaBasedCellPopulation<DIM>::Update(bool hasHadBirthsOrDeaths)
501{
502}
503
504template<unsigned DIM>
506{
507 pPopulationWriter->Visit(this);
508}
509
510template<unsigned DIM>
512{
513 pPopulationCountWriter->Visit(this);
514}
515
516template<unsigned DIM>
518{
519 pPopulationEventWriter->Visit(this);
520}
521
522template<unsigned DIM>
523void CaBasedCellPopulation<DIM>::AcceptCellWriter(boost::shared_ptr<AbstractCellWriter<DIM, DIM> > pCellWriter, CellPtr pCell)
524{
525 pCellWriter->VisitCell(pCell, this);
526}
527
528template<unsigned DIM>
530{
531 if (this->mOutputResultsForChasteVisualizer)
532 {
533 if (!this-> template HasWriter<CellLocationIndexWriter>())
534 {
535 this-> template AddCellWriter<CellLocationIndexWriter>();
536 }
537 }
538
540}
541
542template<unsigned DIM>
544{
545 double cell_volume = 1.0;
546 return cell_volume;
547}
548
549template<unsigned DIM>
550double CaBasedCellPopulation<DIM>::GetWidth(const unsigned& rDimension)
551{
552 double width = this->mrMesh.GetWidth(rDimension);
553 return width;
554}
555
556template<unsigned DIM>
558{
559 // The update rule must be derived from AbstractCaUpdateRule or AbstractCaSwitchingUpdateRule
560 assert(bool(dynamic_cast<AbstractCaUpdateRule<DIM>*>(pUpdateRule.get())) ||
561 bool(dynamic_cast<AbstractCaSwitchingUpdateRule<DIM>*>(pUpdateRule.get())));
562
563 if (bool(dynamic_cast<AbstractCaUpdateRule<DIM>*>(pUpdateRule.get())))
564 {
565 this->mUpdateRuleCollection.push_back(pUpdateRule);
566 }
567 else
568 {
569 mSwitchingUpdateRuleCollection.push_back(pUpdateRule);
570 }
571}
572
573template<unsigned DIM>
575{
576 // Clear mSwitchingUpdateRuleCollection
577 mSwitchingUpdateRuleCollection.clear();
578
579 // Clear mUpdateRuleCollection
581}
582
583template<unsigned DIM>
584const std::vector<boost::shared_ptr<AbstractUpdateRule<DIM> > > CaBasedCellPopulation<DIM>::GetUpdateRuleCollection() const
585{
586 std::vector<boost::shared_ptr<AbstractUpdateRule<DIM> > > update_rules;
587
588 for (unsigned i=0; i<this->mUpdateRuleCollection.size(); i++)
589 {
590 update_rules.push_back(this->mUpdateRuleCollection[i]);
591 }
592 for (unsigned i=0; i<mSwitchingUpdateRuleCollection.size(); i++)
593 {
594 update_rules.push_back(mSwitchingUpdateRuleCollection[i]);
595 }
596
597 return update_rules;
598}
599
600template<unsigned DIM>
601boost::shared_ptr<AbstractCaBasedDivisionRule<DIM> > CaBasedCellPopulation<DIM>::GetCaBasedDivisionRule()
602{
603 return mpCaBasedDivisionRule;
604}
605
606template<unsigned DIM>
608{
609 mpCaBasedDivisionRule = pCaBasedDivisionRule;
610}
611
612template<unsigned DIM>
614{
615 // Add the division rule parameters
616 *rParamsFile << "\t\t<CaBasedDivisionRule>\n";
617 mpCaBasedDivisionRule->OutputCellCaBasedDivisionRuleInfo(rParamsFile);
618 *rParamsFile << "\t\t</CaBasedDivisionRule>\n";
619
620 // Call method on direct parent class
622}
623
624template<unsigned DIM>
625void CaBasedCellPopulation<DIM>::WriteVtkResultsToFile(const std::string& rDirectory)
626{
627#ifdef CHASTE_VTK
628 // Store the present time as a string
629 unsigned num_timesteps = SimulationTime::Instance()->GetTimeStepsElapsed();
630 std::stringstream time;
631 time << num_timesteps;
632
633 // Store the number of cells for which to output data to VTK
634 unsigned num_cells = this->GetNumRealCells();
635
636 // When outputting any CellData, we assume that the first cell is representative of all cells
637 unsigned num_cell_data_items = 0u;
638 std::vector<std::string> cell_data_names;
639 if (num_cells > 0u)
640 {
641 num_cell_data_items = this->Begin()->GetCellData()->GetNumItems();
642 cell_data_names = this->Begin()->GetCellData()->GetKeys();
643 }
644 std::vector<std::vector<double> > cell_data;
645 for (unsigned var=0; var<num_cell_data_items; var++)
646 {
647 std::vector<double> cell_data_var(num_cells);
648 cell_data.push_back(cell_data_var);
649 }
650
651 // Create mesh writer for VTK output
652 VtkMeshWriter<DIM, DIM> mesh_writer(rDirectory, "results_"+time.str(), false);
653
654 // Create a counter to keep track of how many cells are at a lattice site
655 unsigned num_sites = this->mrMesh.GetNumNodes();
656 boost::scoped_array<unsigned> number_of_cells_at_site(new unsigned[num_sites]);
657 for (unsigned i=0; i<num_sites; i++)
658 {
659 number_of_cells_at_site[i] = 0;
660 }
661
662 // Populate a vector of nodes associated with cell locations, by iterating through the list of cells
663 std::vector<Node<DIM>*> nodes;
664 unsigned node_index = 0;
665 for (typename AbstractCellPopulation<DIM,DIM>::Iterator cell_iter = this->Begin();
666 cell_iter != this->End();
667 ++cell_iter)
668 {
669 // Get the location index of this cell and update the counter number_of_cells_at_site
670 unsigned location_index = this->GetLocationIndexUsingCell(*cell_iter);
671 number_of_cells_at_site[location_index]++;
672 assert(number_of_cells_at_site[location_index] <= mLatticeCarryingCapacity);
673
674 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
675 c_vector<double, DIM> coords;
676 coords = this->mrMesh.GetNode(location_index)->rGetLocation();
677
678 // Move the coordinates slightly so that we can visualise all cells in a lattice site if there is more than one per site
679 if (mLatticeCarryingCapacity > 1)
680 {
681 c_vector<double, DIM> offset;
682
683 if (DIM == 2)
684 {
685 double angle = (double)number_of_cells_at_site[location_index]*2.0*M_PI/(double)mLatticeCarryingCapacity;
686 offset[0] = 0.2*sin(angle);
687 offset[1] = 0.2*cos(angle);
688 }
689 else
690 {
692 for (unsigned i=0; i<DIM; i++)
693 {
694 offset[i] = p_gen->ranf(); // This assumes that all sites are 1 unit apart
695 }
696 }
697
698 for (unsigned i=0; i<DIM; i++)
699 {
700 coords[i] += offset[i];
701 }
702 }
703
704 nodes.push_back(new Node<DIM>(node_index, coords, false));
705 node_index++;
706 }
707
708 // Iterate over any cell writers that are present
709 for (typename std::vector<boost::shared_ptr<AbstractCellWriter<DIM, DIM> > >::iterator cell_writer_iter = this->mCellWriters.begin();
710 cell_writer_iter != this->mCellWriters.end();
711 ++cell_writer_iter)
712 {
713 // Create vector to store VTK cell data
714 // (using a default value of -1 to correspond to an empty lattice site)
715 std::vector<double> vtk_cell_data(num_cells, -1.0);
716
717 // Loop over cells
718 unsigned cell_index = 0;
719 for (typename AbstractCellPopulation<DIM,DIM>::Iterator cell_iter = this->Begin();
720 cell_iter != this->End();
721 ++cell_iter)
722 {
723 // Populate the vector of VTK cell data
724 vtk_cell_data[cell_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(*cell_iter, this);
725 cell_index++;
726 }
727
728 mesh_writer.AddPointData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
729 }
730
731 // Loop over cells to output the cell data
732 unsigned cell_index = 0;
733 for (typename AbstractCellPopulation<DIM,DIM>::Iterator cell_iter = this->Begin();
734 cell_iter != this->End();
735 ++cell_iter)
736 {
737 for (unsigned var=0; var<num_cell_data_items; var++)
738 {
739 cell_data[var][cell_index] = cell_iter->GetCellData()->GetItem(cell_data_names[var]);
740 }
741 cell_index++;
742 }
743 for (unsigned var=0; var<num_cell_data_items; var++)
744 {
745 mesh_writer.AddPointData(cell_data_names[var], cell_data[var]);
746 }
747
748 /*
749 * At present, the VTK writer can only write things which inherit from AbstractTetrahedralMeshWriter.
750 * For now, we do an explicit conversion to NodesOnlyMesh. This can be written to VTK, then visualized
751 * as glyphs in Paraview.
752 */
753 NodesOnlyMesh<DIM> temp_mesh;
754
755 // Use an approximation of the node spacing as the interaction distance for the nodes only mesh. This is to
756 // avoid rounding errors in distributed box collection.
757 double volume = this->mrMesh.GetWidth(0);
758 for (unsigned idx=1; idx<DIM; idx++)
759 {
760 volume *= this->mrMesh.GetWidth(idx);
761 }
762
763 double spacing;
764 if (this->mrMesh.GetNumNodes() >0 && volume > 0.0)
765 {
766 spacing = std::pow(volume / double(this->mrMesh.GetNumNodes()), 1.0/double(DIM));
767 }
768 else
769 {
770 spacing = 1.0;
771 }
772
773 temp_mesh.ConstructNodesWithoutMesh(nodes, spacing * 1.2);
774 mesh_writer.WriteFilesUsingMesh(temp_mesh);
775
776 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
777 *(this->mpVtkMetaFile) << num_timesteps;
778 *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
779 *(this->mpVtkMetaFile) << num_timesteps;
780 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
781
782 // Tidy up
783 for (unsigned i=0; i<nodes.size(); i++)
784 {
785 delete nodes[i];
786 }
787#endif //CHASTE_VTK
788}
789
790template<unsigned DIM>
792 unsigned pdeNodeIndex,
793 std::string& rVariableName,
794 bool dirichletBoundaryConditionApplies,
795 double dirichletBoundaryValue)
796{
797 unsigned counter = 0;
798 typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
799 while (counter != pdeNodeIndex)
800 {
801 ++cell_iter;
802 counter++;
803 }
804
805 double value = cell_iter->GetCellData()->GetItem(rVariableName);
806
807 return value;
808}
809
810template<unsigned DIM>
812{
813 // pdeNodeIndex corresponds to the 'position' of the cell to interrogate in the vector of cells
814 typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
815
816 assert(pdeNodeIndex < this->GetNumRealCells());
817 for (unsigned i=0; i<pdeNodeIndex; i++)
818 {
819 ++cell_iter;
820 }
821 bool is_cell_apoptotic = cell_iter->template HasCellProperty<ApoptoticCellProperty>();
822
823 return !is_cell_apoptotic;
824}
825
826// Explicit instantiation
827template class CaBasedCellPopulation<1>;
828template class CaBasedCellPopulation<2>;
829template class CaBasedCellPopulation<3>;
830
831// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define NEVER_REACHED
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
virtual void AddCellUsingLocationIndex(unsigned index, CellPtr pCell)
virtual void OpenWritersFiles(OutputFileHandler &rOutputFileHandler)
virtual void RemoveCellUsingLocationIndex(unsigned index, CellPtr pCell)
virtual void OutputCellPopulationParameters(out_stream &rParamsFile)
void SetCaBasedDivisionRule(boost::shared_ptr< AbstractCaBasedDivisionRule< DIM > > pCaBasedDivisionRule)
Node< DIM > * GetNode(unsigned index)
virtual void AddUpdateRule(boost::shared_ptr< AbstractUpdateRule< DIM > > pUpdateRule)
CaBasedCellPopulation(PottsMesh< DIM > &rMesh, std::vector< CellPtr > &rCells, const std::vector< unsigned > locationIndices, unsigned latticeCarryingCapacity=1u, bool deleteMesh=false, bool validate=false)
std::vector< unsigned > mAvailableSpaces
boost::shared_ptr< AbstractCaBasedDivisionRule< DIM > > GetCaBasedDivisionRule()
virtual double EvaluateDivisionPropensity(unsigned currentNodeIndex, unsigned targetNodeIndex, CellPtr pCell)
void AddCellUsingLocationIndex(unsigned index, CellPtr pCell)
virtual const std::vector< boost::shared_ptr< AbstractUpdateRule< DIM > > > GetUpdateRuleCollection() const
CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell=CellPtr())
c_vector< double, DIM > GetLocationOfCellCentre(CellPtr pCell)
virtual void WriteVtkResultsToFile(const std::string &rDirectory)
double GetWidth(const unsigned &rDimension)
virtual void AcceptCellWriter(boost::shared_ptr< AbstractCellWriter< DIM, DIM > > pCellWriter, CellPtr pCell)
std::set< unsigned > GetNeighbouringLocationIndices(CellPtr pCell)
void OutputCellPopulationParameters(out_stream &rParamsFile)
void RemoveCellUsingLocationIndex(unsigned index, CellPtr pCell)
virtual bool IsPdeNodeAssociatedWithNonApoptoticCell(unsigned pdeNodeIndex)
std::vector< unsigned > & rGetAvailableSpaces()
virtual void AcceptPopulationCountWriter(boost::shared_ptr< AbstractCellPopulationCountWriter< DIM, DIM > > pPopulationCountWriter)
boost::shared_ptr< AbstractCaBasedDivisionRule< DIM > > mpCaBasedDivisionRule
PottsMesh< DIM > & rGetMesh()
virtual void AcceptPopulationWriter(boost::shared_ptr< AbstractCellPopulationWriter< DIM, DIM > > pPopulationWriter)
void Update(bool hasHadBirthsOrDeaths=true)
virtual double GetCellDataItemAtPdeNode(unsigned pdeNodeIndex, std::string &rVariableName, bool dirichletBoundaryConditionApplies=false, double dirichletBoundaryValue=0.0)
bool IsCellAssociatedWithADeletedLocation(CellPtr pCell)
virtual void AcceptPopulationEventWriter(boost::shared_ptr< AbstractCellPopulationEventWriter< DIM, DIM > > pPopulationEventWriter)
Node< DIM > * GetNodeCorrespondingToCell(CellPtr pCell)
virtual TetrahedralMesh< DIM, DIM > * GetTetrahedralMeshForPdeModifier()
bool IsRoomToDivide(CellPtr pCell)
virtual void OpenWritersFiles(OutputFileHandler &rOutputFileHandler)
virtual bool IsSiteAvailable(unsigned index, CellPtr pCell)
double GetVolumeOfCell(CellPtr pCell)
Definition Node.hpp:59
void ConstructNodesWithoutMesh(const std::vector< Node< SPACE_DIM > * > &rNodes, double maxInteractionDistance)
double GetWidth(const unsigned &rDimension) const
std::set< unsigned > GetMooreNeighbouringNodeIndices(unsigned nodeIndex)
static RandomNumberGenerator * Instance()
void Shuffle(std::vector< boost::shared_ptr< T > > &rValues)
unsigned randMod(unsigned base)
static SimulationTime * Instance()
unsigned GetTimeStepsElapsed() const
void AddPointData(std::string name, std::vector< double > data)
void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)