Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ShovingCaBasedDivisionRule.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 "ShovingCaBasedDivisionRule.hpp"
37#include "RandomNumberGenerator.hpp"
38
39template<unsigned SPACE_DIM>
41{
42 // This logic currently only works for Moore neighbourhood in 2D
43 bool is_central_node = false;
44 if (SPACE_DIM == 2)
45 {
46 if (numNeighbours == 8)
47 {
48 is_central_node = true;
49 }
50 }
51 else // Currently not tested in 3D
52 {
54 }
55 if (!is_central_node)
56 {
57 EXCEPTION("Cells reaching the boundary of the domain. Make the Potts mesh larger.");
58 }
59}
60
61template<unsigned SPACE_DIM>
63{
64 return true;
65}
66
67template<unsigned SPACE_DIM>
69 CellPtr pParentCell,
70 CaBasedCellPopulation<SPACE_DIM>& rCellPopulation)
71{
72 // Get node index corresponding to the parent cell
73 unsigned parent_node_index = rCellPopulation.GetLocationIndexUsingCell(pParentCell);
74
75 PottsMesh<SPACE_DIM>* p_static_cast_mesh = static_cast<PottsMesh<SPACE_DIM>*>(&(rCellPopulation.rGetMesh()));
76
77 // Get the set of neighbouring node indices
78 std::set<unsigned> neighbouring_node_indices = p_static_cast_mesh->GetMooreNeighbouringNodeIndices(parent_node_index);
79 unsigned num_neighbours = neighbouring_node_indices.size();
80
81 // Check cell is not on the boundary
82 IsNodeOnBoundary(num_neighbours);
83
84 std::vector<double> neighbouring_node_propensities;
85 std::vector<unsigned> neighbouring_node_indices_vector;
86
87 double total_propensity = 0.0;
88
89 // Select neighbour at random
90 for (std::set<unsigned>::iterator neighbour_iter = neighbouring_node_indices.begin();
91 neighbour_iter != neighbouring_node_indices.end();
92 ++neighbour_iter)
93 {
94 neighbouring_node_indices_vector.push_back(*neighbour_iter);
95
96 double propensity_dividing_into_neighbour = rCellPopulation.EvaluateDivisionPropensity(parent_node_index,*neighbour_iter,pParentCell);
97
98 neighbouring_node_propensities.push_back(propensity_dividing_into_neighbour);
99 total_propensity += propensity_dividing_into_neighbour;
100 }
101 assert(total_propensity>0); // if this trips the cell can't divided so need to include this in the IsSiteAvailable method
102
103 for (unsigned i=0; i<num_neighbours; i++)
104 {
105 neighbouring_node_propensities[i] /= total_propensity;
106 }
107
108 // Sample random number to specify which move to make
110 double random_number = p_gen->ranf();
111
112 double total_probability = 0.0;
113 unsigned daughter_node_index = UNSIGNED_UNSET;
114
115 unsigned counter;
116 for (counter=0; counter < num_neighbours; counter++)
117 {
118 total_probability += neighbouring_node_propensities[counter];
119 if (total_probability >= random_number)
120 {
121 // Divide the parent cell to this neighbour location
122 daughter_node_index = neighbouring_node_indices_vector[counter];
123 break;
124 }
125 }
126 // This loop should always break as sum(neighbouring_node_propensities) = 1
127
128 assert(daughter_node_index != UNSIGNED_UNSET);
129 assert(daughter_node_index < p_static_cast_mesh->GetNumNodes());
130
131 // If daughter node is occupied then move the cell in the direction of counter
132 if (!(rCellPopulation.IsSiteAvailable(daughter_node_index, pNewCell)))
133 {
134 std::list<std::pair<unsigned,unsigned> > cell_moves;
135
136 bool is_neighbour_occupied = true;
137 int max_moves = 1000;
138 int move_number = 0;
139 unsigned current_node_index = parent_node_index;
140 unsigned target_node_index = daughter_node_index;
141 while (is_neighbour_occupied && move_number < max_moves)
142 {
143 move_number++;
144 current_node_index = target_node_index;
145
146 std::set<unsigned> neighbouring_node_indices = p_static_cast_mesh->GetMooreNeighbouringNodeIndices(current_node_index);
147
148 // Check cell is not on the boundary
149 IsNodeOnBoundary(neighbouring_node_indices.size());
150
151 // Select the appropriate neighbour
152 std::set<unsigned>::iterator neighbour_iter = neighbouring_node_indices.begin();
153 for (unsigned i=0; i<counter; i++)
154 {
155 ++neighbour_iter;
156 }
157 assert(neighbour_iter != neighbouring_node_indices.end());
158
159 target_node_index = *neighbour_iter;
160
161 std::pair<unsigned, unsigned> new_move(current_node_index, target_node_index);
162
163 cell_moves.push_back(new_move);
164
165 // If target node is unoccupied move the cell on the current node to the target node and stop shoving cells
166 if (rCellPopulation.IsSiteAvailable(target_node_index, pNewCell))
167 {
168 is_neighbour_occupied = false;
169 }
170
171 // If target node is occupied then keep shoving the cells out of the way
172 current_node_index = target_node_index;
173
174 }
175 assert(move_number<max_moves);
176
177 // Do moves to free up the daughter node index
178 for (std::list<std::pair<unsigned, unsigned> >::reverse_iterator reverse_iter = cell_moves.rbegin();
179 reverse_iter != cell_moves.rend();
180 ++reverse_iter)
181 {
182 assert(rCellPopulation.IsSiteAvailable(reverse_iter->second, pNewCell));
183 assert(!(rCellPopulation.IsSiteAvailable(reverse_iter->first, pNewCell)));
184
185 // Move cell from first() to second()
186 rCellPopulation.MoveCellInLocationMap(rCellPopulation.GetCellUsingLocationIndex(reverse_iter->first), reverse_iter->first, reverse_iter->second);
187 }
188
189 // Check daughter site is now free
190 assert(rCellPopulation.IsSiteAvailable(daughter_node_index, pNewCell));
191
192 }
193 return daughter_node_index;
194}
195
196// Explicit instantiation
197template class ShovingCaBasedDivisionRule<1>;
198template class ShovingCaBasedDivisionRule<2>;
199template class ShovingCaBasedDivisionRule<3>;
200
201// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define NEVER_REACHED
const unsigned UNSIGNED_UNSET
Definition Exception.hpp:53
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
unsigned GetLocationIndexUsingCell(CellPtr pCell)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
void MoveCellInLocationMap(CellPtr pCell, unsigned old_index, unsigned new_index)
virtual double EvaluateDivisionPropensity(unsigned currentNodeIndex, unsigned targetNodeIndex, CellPtr pCell)
PottsMesh< DIM > & rGetMesh()
virtual bool IsSiteAvailable(unsigned index, CellPtr pCell)
std::set< unsigned > GetMooreNeighbouringNodeIndices(unsigned nodeIndex)
static RandomNumberGenerator * Instance()
virtual unsigned CalculateDaughterNodeIndex(CellPtr pNewCell, CellPtr pParentCell, CaBasedCellPopulation< SPACE_DIM > &rCellPopulation)
virtual bool IsRoomToDivide(CellPtr pParentCell, CaBasedCellPopulation< SPACE_DIM > &rCellPopulation)
void IsNodeOnBoundary(unsigned numNeighbours)