Chaste  Release::2018.1
HeterotypicBoundaryLengthWriter.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 "HeterotypicBoundaryLengthWriter.hpp"
37 
38 #include "AbstractCellPopulation.hpp"
39 #include "MeshBasedCellPopulation.hpp"
40 #include "CaBasedCellPopulation.hpp"
41 #include "NodeBasedCellPopulation.hpp"
42 #include "PottsBasedCellPopulation.hpp"
43 #include "VertexBasedCellPopulation.hpp"
44 
45 #include "CellLabel.hpp"
46 
47 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
49  : AbstractCellPopulationWriter<ELEMENT_DIM, SPACE_DIM>("heterotypicboundary.dat")
50 {
51 }
52 
53 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
55 {
56  // Initialise helper variables
57  double heterotypic_boundary_length = 0.0;
58  double total_shared_edges_length = 0.0;
59  double num_heterotypic_pairs = 0.0;
60  double total_num_pairs = 0.0;
61 
62  // Make sure the Voronoi tessellation has been created
64  pCellPopulation->CreateVoronoiTessellation();
65 
66  // Iterate over cells
67  for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
68  cell_iter != pCellPopulation->End();
69  ++cell_iter)
70  {
71  // Get the location index corresponding to this cell
72  unsigned index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
73 
74  // Store whether this cell is labelled
75  bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
76 
77  // Get the set of neighbouring location indices
78  std::set<unsigned> neighbour_indices = pCellPopulation->GetNeighbouringNodeIndices(index);
79 
80  // Iterate over these neighbours
81  for (std::set<unsigned>::iterator neighbour_iter = neighbour_indices.begin();
82  neighbour_iter != neighbour_indices.end();
83  ++neighbour_iter)
84  {
85  // Get the length of the edge shared with this neighbour
86  unsigned neighbour_index = *neighbour_iter;
87  double edge_length = pCellPopulation->GetVoronoiEdgeLength(index,neighbour_index);
88 
89  // Ignore ghost nodes (in the case of a MeshBasedCellPopulationWithGhostNodes)
91  if (!pCellPopulation->IsGhostNode(neighbour_index))
92  {
93  total_shared_edges_length += edge_length;
94  total_num_pairs += 1.0;
95 
96  // Store whether this neighbour is labelled
97  CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(neighbour_index);
98  bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
99 
100  // If this cell is labelled and its neighbour is not, or vice versa...
101  if (cell_is_labelled != neighbour_is_labelled)
102  {
103  // ... then increment the fractional boundary length
104  heterotypic_boundary_length += edge_length;
105  num_heterotypic_pairs += 1.0;
106  }
107  }
108  }
109  }
110 
111  // We have counted each cell-cell edge twice
112  heterotypic_boundary_length *= 0.5;
113  total_shared_edges_length *= 0.5;
114 
115  // We have counted each pair of neighbouring cells twice
116  num_heterotypic_pairs *= 0.5;
117  total_num_pairs *= 0.5;
118 
119  *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
120 }
121 
122 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
124 {
125  // Initialise helper variables
126  double heterotypic_boundary_length = 0.0;
127  double total_shared_edges_length = 0.0;
128  double num_heterotypic_pairs = 0.0;
129  double total_num_pairs = 0.0;
130 
131  // Iterate over cells
132  for (typename AbstractCellPopulation<SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
133  cell_iter != pCellPopulation->End();
134  ++cell_iter)
135  {
136  // Get the location index corresponding to this cell
137  unsigned index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
138 
139  // Store whether this cell is labelled
140  bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
141 
142  // Get this node's von Neumann neighbours (not Moore neighbours, since they must share an edge)
143  std::set<unsigned> neighbour_node_indices = pCellPopulation->rGetMesh().GetVonNeumannNeighbouringNodeIndices(index);
144 
145  // Iterate over these neighbours
146  for (std::set<unsigned>::iterator neighbour_iter = neighbour_node_indices.begin();
147  neighbour_iter != neighbour_node_indices.end();
148  ++neighbour_iter)
149  {
150  // Assume the lattice is comprised of uniform unit lattice sites
151  double edge_length = 1.0;
152 
153  unsigned neighbour_index = *neighbour_iter;
154 
155  if (pCellPopulation->IsCellAttachedToLocationIndex(neighbour_index))
156  {
157  CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(neighbour_index);
158 
159  total_shared_edges_length += edge_length;
160  total_num_pairs += 1.0;
161 
162  // Store whether this neighbour is labelled
163  bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
164 
165  // If this cell is labelled and its neighbour is not, or vice versa...
166  if (cell_is_labelled != neighbour_is_labelled)
167  {
168  // ... then increment the fractional boundary length
169  heterotypic_boundary_length += edge_length;
170  num_heterotypic_pairs += 1.0;
171  }
172  }
173  }
174  }
175 
176  // We have counted each cell-cell edge twice
177  heterotypic_boundary_length *= 0.5;
178  total_shared_edges_length *= 0.5;
179 
180  // We have counted each pair of neighbouring cells twice
181  num_heterotypic_pairs *= 0.5;
182  total_num_pairs *= 0.5;
183 
184  *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
185 }
186 
187 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
189 {
190  // Make sure the cell population is updated so that mNodeNeighbours is set up
192  pCellPopulation->Update();
193 
194  // Initialise helper variables
195  double heterotypic_boundary_length = 0.0;
196  double total_shared_edges_length = 0.0;
197  double num_heterotypic_pairs = 0.0;
198  double total_num_pairs = 0.0;
199 
200  // Loop over cells
201  for (typename AbstractCellPopulation<SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
202  cell_iter != pCellPopulation->End();
203  ++cell_iter)
204  {
205  // Store whether this cell is labelled
206  bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
207 
208  // Store the radius of the node corresponding to this cell
209  unsigned node_index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
210  double node_radius = pCellPopulation->GetNode(node_index)->GetRadius();
211 
212  // Get the set of neighbouring node indices
213  std::set<unsigned> neighbour_indices = pCellPopulation->GetNeighbouringNodeIndices(node_index);
214 
215  if (!neighbour_indices.empty())
216  {
217  // Iterate over these neighbours
218  for (std::set<unsigned>::iterator neighbour_iter = neighbour_indices.begin();
219  neighbour_iter != neighbour_indices.end();
220  ++neighbour_iter)
221  {
222  // Store the radius of the node corresponding to this neighbour
223  double neighbour_radius = pCellPopulation->GetNode(*neighbour_iter)->GetRadius();
224 
225  // Get the (approximate) length of the edge shared with this neighbour
226  double separation = pCellPopulation->rGetMesh().GetDistanceBetweenNodes(node_index, *neighbour_iter);
227  double sum_of_radii = node_radius + neighbour_radius;
228 
229  // If the neighbours are close enough, then approximate their 'edge length'
230  if (separation < sum_of_radii)
231  {
232  // Use Heron's formula to compute the edge length
233  double a = node_radius;
234  double b = neighbour_radius;
235  double c = separation;
236  double s = 0.5*(a + b + c);
237  double A = sqrt(s*(s-a)*(s-b)*(s-c));
238  double edge_length = 4.0*A/c;
239 
240  total_shared_edges_length += edge_length;
241  total_num_pairs += 1.0;
242 
243  // Store whether this neighbour is labelled
244  CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(*neighbour_iter);
245  bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
246 
247  // If this cell is labelled and its neighbour is not, or vice versa...
248  if (cell_is_labelled != neighbour_is_labelled)
249  {
250  // ... then increment the fractional boundary length
251  heterotypic_boundary_length += edge_length;
252  num_heterotypic_pairs += 1.0;
253  }
254  }
255  }
256  }
257  }
258 
259  // We have counted each cell-cell edge twice
260  heterotypic_boundary_length *= 0.5;
261  total_shared_edges_length *= 0.5;
262 
263  // We have counted each pair of neighbouring cells twice
264  num_heterotypic_pairs *= 0.5;
265  total_num_pairs *= 0.5;
266 
267  *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
268 }
269 
270 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
272 {
274 
275  // Initialise helper variables
276  double heterotypic_boundary_length = 0.0;
277  double total_shared_edges_length = 0.0;
278  double num_heterotypic_pairs = 0.0;
279  double total_num_pairs = 0.0;
280 
281  // Iterate over cells
282  for (typename AbstractCellPopulation<SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
283  cell_iter != pCellPopulation->End();
284  ++cell_iter)
285  {
286  // Store whether this cell is labelled
287  bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
288 
289  unsigned elem_index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
290  unsigned num_nodes_in_elem = pCellPopulation->rGetMesh().GetElement(elem_index)->GetNumNodes();
291 
292  // Iterate over nodes contained in the element corresponding to this cell
293  for (unsigned local_index=0; local_index<num_nodes_in_elem; local_index++)
294  {
295  // Get this node's von Neumann neighbours (not Moore neighbours, since they must share an edge)
296  unsigned global_index = pCellPopulation->rGetMesh().GetElement(elem_index)->GetNodeGlobalIndex(local_index);
297  std::set<unsigned> neighbour_node_indices = pCellPopulation->rGetMesh().GetVonNeumannNeighbouringNodeIndices(global_index);
298 
299  // Iterate over these neighbours
300  for (std::set<unsigned>::iterator neighbour_iter = neighbour_node_indices.begin();
301  neighbour_iter != neighbour_node_indices.end();
302  ++neighbour_iter)
303  {
304  // Get the elements containing this neighbour
305  std::set<unsigned> neighbour_elem_indices = pCellPopulation->GetNode(*neighbour_iter)->rGetContainingElementIndices();
306 
307  if (neighbour_elem_indices.size() == 1)
308  {
309  unsigned neigbour_elem_index = *(neighbour_elem_indices.begin());
310 
311  if (neigbour_elem_index != elem_index)
312  {
313  // Edge is between two different elements
314  total_shared_edges_length += 1.0;
315 
316  // Store whether the cell corresponding to this element index is labelled
317  CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(neigbour_elem_index);
318  bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
319 
320  // If this cell is labelled and its neighbour is not, or vice versa...
321  if (cell_is_labelled != neighbour_is_labelled)
322  {
323  // ... then increment the fractional boundary length
324  heterotypic_boundary_length += 1.0;
325  }
326  }
327  }
328  else
329  {
330  // Original node is on boundary of mesh so don't include in the edge calculation.
331  }
332  }
333  }
334 
335  std::set<unsigned> neighbour_node_indices = pCellPopulation->GetNeighbouringLocationIndices(*cell_iter);
336 
337  // Iterate over these neighbours
338  for (std::set<unsigned>::iterator neighbour_iter = neighbour_node_indices.begin();
339  neighbour_iter != neighbour_node_indices.end();
340  ++neighbour_iter)
341  {
342  total_num_pairs += 1.0;
343 
344  CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(*neighbour_iter);
345  bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
346 
347  // If this cell is labelled and its neighbour is not, or vice versa...
348  if (cell_is_labelled != neighbour_is_labelled)
349  {
350  // ... then increment the fractional boundary length
351  num_heterotypic_pairs += 1.0;
352  }
353  }
354  }
355 
356  // We have counted each cell-cell edge twice
357  heterotypic_boundary_length *= 0.5;
358  total_shared_edges_length *= 0.5;
359 
360  // We have counted each pair of neighbouring cells twice
361  num_heterotypic_pairs *= 0.5;
362  total_num_pairs *= 0.5;
363 
364  *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
365 }
366 
367 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
369 {
370  // Make sure the cell population is updated
372  pCellPopulation->Update();
373 
374  // Initialise helper variables
375  double heterotypic_boundary_length = 0.0;
376  double total_shared_edges_length = 0.0;
377  double num_heterotypic_pairs = 0.0;
378  double total_num_pairs = 0.0;
379 
380  // Iterate over cells
381  for (typename AbstractCellPopulation<SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
382  cell_iter != pCellPopulation->End();
383  ++cell_iter)
384  {
385  // Store whether this cell is labelled
386  bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
387 
388  // Get the set of neighbouring element indices
389  unsigned elem_index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
390  std::set<unsigned> neighbour_elem_indices = pCellPopulation->rGetMesh().GetNeighbouringElementIndices(elem_index);
391 
392  // Iterate over these neighbours
393  for (std::set<unsigned>::iterator neighbour_iter = neighbour_elem_indices.begin();
394  neighbour_iter != neighbour_elem_indices.end();
395  ++neighbour_iter)
396  {
397  // Get the length of the edge shared with this neighbour
398  unsigned neighbour_index = *neighbour_iter;
399  double edge_length = pCellPopulation->rGetMesh().GetEdgeLength(elem_index, neighbour_index);
400 
401  total_shared_edges_length += edge_length;
402  total_num_pairs += 1.0;
403 
404  // Store whether this neighbour is labelled
405  CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(*neighbour_iter);
406  bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
407 
408  // If this cell is labelled and its neighbour is not, or vice versa...
409  if (cell_is_labelled != neighbour_is_labelled)
410  {
411  // ... then increment the fractional boundary length
412  heterotypic_boundary_length += edge_length;
413  num_heterotypic_pairs += 1.0;
414  }
415  }
416  }
417 
418  // We have counted each cell-cell edge twice
419  heterotypic_boundary_length *= 0.5;
420  total_shared_edges_length *= 0.5;
421 
422  // We have counted each pair of neighbouring cells twice
423  num_heterotypic_pairs *= 0.5;
424  total_num_pairs *= 0.5;
425 
426  *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
427 }
428 
429 // Explicit instantiation
436 
438 // Declare identifier for the serializer
double GetVoronoiEdgeLength(unsigned index1, unsigned index2)
NodesOnlyMesh< DIM > & rGetMesh()
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
PottsMesh< DIM > & rGetMesh()
unsigned GetLocationIndexUsingCell(CellPtr pCell)
Node< DIM > * GetNode(unsigned index)
Node< DIM > * GetNode(unsigned index)
std::set< unsigned > GetNeighbouringLocationIndices(CellPtr pCell)
virtual bool IsCellAttachedToLocationIndex(unsigned index)
double GetDistanceBetweenNodes(unsigned indexA, unsigned indexB)
MutableVertexMesh< DIM, DIM > & rGetMesh()
std::set< unsigned > GetNeighbouringNodeIndices(unsigned index)
void Update(bool hasHadBirthsOrDeaths=true)
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
std::set< unsigned > GetNeighbouringNodeIndices(unsigned index)
virtual void Visit(MeshBasedCellPopulation< ELEMENT_DIM, SPACE_DIM > *pCellPopulation)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
std::set< unsigned > GetNeighbouringElementIndices(unsigned elementIndex)
Definition: VertexMesh.cpp:783
void Update(bool hasHadBirthsOrDeaths=true)
double GetEdgeLength(unsigned elementIndex1, unsigned elementIndex2)
Definition: VertexMesh.cpp:415