Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
HeterotypicBoundaryLengthWriter.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 "HeterotypicBoundaryLengthWriter.hpp"
37
38#include "AbstractCellPopulation.hpp"
39#include "CaBasedCellPopulation.hpp"
40#include "ImmersedBoundaryCellPopulation.hpp"
41#include "MeshBasedCellPopulation.hpp"
42#include "NodeBasedCellPopulation.hpp"
43#include "PottsBasedCellPopulation.hpp"
44#include "VertexBasedCellPopulation.hpp"
45
46#include "CellLabel.hpp"
47
48#include <boost/polygon/voronoi.hpp>
49#include <set>
50#include <utility>
51
52template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
57
58template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
60{
61 // Initialise helper variables
62 double heterotypic_boundary_length = 0.0;
63 double total_shared_edges_length = 0.0;
64 double num_heterotypic_pairs = 0.0;
65 double total_num_pairs = 0.0;
66
67 // Make sure the Voronoi tessellation has been created
69 pCellPopulation->CreateVoronoiTessellation();
70
71 // Iterate over cells
72 for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
73 cell_iter != pCellPopulation->End();
74 ++cell_iter)
75 {
76 // Get the location index corresponding to this cell
77 unsigned index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
78
79 // Store whether this cell is labelled
80 bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
81
82 // Get the set of neighbouring location indices
83 std::set<unsigned> neighbour_indices = pCellPopulation->GetNeighbouringNodeIndices(index);
84
85 // Iterate over these neighbours
86 for (std::set<unsigned>::iterator neighbour_iter = neighbour_indices.begin();
87 neighbour_iter != neighbour_indices.end();
88 ++neighbour_iter)
89 {
90 // Get the length of the edge shared with this neighbour
91 unsigned neighbour_index = *neighbour_iter;
92 double edge_length = pCellPopulation->GetVoronoiEdgeLength(index,neighbour_index);
93
94 // Ignore ghost nodes (in the case of a MeshBasedCellPopulationWithGhostNodes)
96 if (!pCellPopulation->IsGhostNode(neighbour_index))
97 {
98 total_shared_edges_length += edge_length;
99 total_num_pairs += 1.0;
100
101 // Store whether this neighbour is labelled
102 CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(neighbour_index);
103 bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
104
105 // If this cell is labelled and its neighbour is not, or vice versa...
106 if (cell_is_labelled != neighbour_is_labelled)
107 {
108 // ... then increment the fractional boundary length
109 heterotypic_boundary_length += edge_length;
110 num_heterotypic_pairs += 1.0;
111 }
112 }
113 }
114 }
115
116 // We have counted each cell-cell edge twice
117 heterotypic_boundary_length *= 0.5;
118 total_shared_edges_length *= 0.5;
119
120 // We have counted each pair of neighbouring cells twice
121 num_heterotypic_pairs *= 0.5;
122 total_num_pairs *= 0.5;
123
124 *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
125}
126
127template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
129{
130 // Initialise helper variables
131 double heterotypic_boundary_length = 0.0;
132 double total_shared_edges_length = 0.0;
133 double num_heterotypic_pairs = 0.0;
134 double total_num_pairs = 0.0;
135
136 // Iterate over cells
137 for (typename AbstractCellPopulation<SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
138 cell_iter != pCellPopulation->End();
139 ++cell_iter)
140 {
141 // Get the location index corresponding to this cell
142 unsigned index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
143
144 // Store whether this cell is labelled
145 bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
146
147 // Get this node's von Neumann neighbours (not Moore neighbours, since they must share an edge)
148 std::set<unsigned> neighbour_node_indices = pCellPopulation->rGetMesh().GetVonNeumannNeighbouringNodeIndices(index);
149
150 // Iterate over these neighbours
151 for (std::set<unsigned>::iterator neighbour_iter = neighbour_node_indices.begin();
152 neighbour_iter != neighbour_node_indices.end();
153 ++neighbour_iter)
154 {
155 // Assume the lattice is comprised of uniform unit lattice sites
156 double edge_length = 1.0;
157
158 unsigned neighbour_index = *neighbour_iter;
159
160 if (pCellPopulation->IsCellAttachedToLocationIndex(neighbour_index))
161 {
162 CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(neighbour_index);
163
164 total_shared_edges_length += edge_length;
165 total_num_pairs += 1.0;
166
167 // Store whether this neighbour is labelled
168 bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
169
170 // If this cell is labelled and its neighbour is not, or vice versa...
171 if (cell_is_labelled != neighbour_is_labelled)
172 {
173 // ... then increment the fractional boundary length
174 heterotypic_boundary_length += edge_length;
175 num_heterotypic_pairs += 1.0;
176 }
177 }
178 }
179 }
180
181 // We have counted each cell-cell edge twice
182 heterotypic_boundary_length *= 0.5;
183 total_shared_edges_length *= 0.5;
184
185 // We have counted each pair of neighbouring cells twice
186 num_heterotypic_pairs *= 0.5;
187 total_num_pairs *= 0.5;
188
189 *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
190}
191
192template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
194{
195 // Make sure the cell population is updated so that mNodeNeighbours is set up
197 pCellPopulation->Update();
198
199 // Initialise helper variables
200 double heterotypic_boundary_length = 0.0;
201 double total_shared_edges_length = 0.0;
202 double num_heterotypic_pairs = 0.0;
203 double total_num_pairs = 0.0;
204
205 // Loop over cells
206 for (typename AbstractCellPopulation<SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
207 cell_iter != pCellPopulation->End();
208 ++cell_iter)
209 {
210 // Store whether this cell is labelled
211 bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
212
213 // Store the radius of the node corresponding to this cell
214 unsigned node_index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
215 double node_radius = pCellPopulation->GetNode(node_index)->GetRadius();
216
217 // Get the set of neighbouring node indices
218 std::set<unsigned> neighbour_indices = pCellPopulation->GetNeighbouringNodeIndices(node_index);
219
220 if (!neighbour_indices.empty())
221 {
222 // Iterate over these neighbours
223 for (std::set<unsigned>::iterator neighbour_iter = neighbour_indices.begin();
224 neighbour_iter != neighbour_indices.end();
225 ++neighbour_iter)
226 {
227 // Store the radius of the node corresponding to this neighbour
228 double neighbour_radius = pCellPopulation->GetNode(*neighbour_iter)->GetRadius();
229
230 // Get the (approximate) length of the edge shared with this neighbour
231 double separation = pCellPopulation->rGetMesh().GetDistanceBetweenNodes(node_index, *neighbour_iter);
232 double sum_of_radii = node_radius + neighbour_radius;
233
234 // If the neighbours are close enough, then approximate their 'edge length'
235 if (separation < sum_of_radii)
236 {
237 // Use Heron's formula to compute the edge length
238 double a = node_radius;
239 double b = neighbour_radius;
240 double c = separation;
241 double s = 0.5*(a + b + c);
242 double A = sqrt(s*(s-a)*(s-b)*(s-c));
243 double edge_length = 4.0*A/c;
244
245 total_shared_edges_length += edge_length;
246 total_num_pairs += 1.0;
247
248 // Store whether this neighbour is labelled
249 CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(*neighbour_iter);
250 bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
251
252 // If this cell is labelled and its neighbour is not, or vice versa...
253 if (cell_is_labelled != neighbour_is_labelled)
254 {
255 // ... then increment the fractional boundary length
256 heterotypic_boundary_length += edge_length;
257 num_heterotypic_pairs += 1.0;
258 }
259 }
260 }
261 }
262 }
263
264 // We have counted each cell-cell edge twice
265 heterotypic_boundary_length *= 0.5;
266 total_shared_edges_length *= 0.5;
267
268 // We have counted each pair of neighbouring cells twice
269 num_heterotypic_pairs *= 0.5;
270 total_num_pairs *= 0.5;
271
272 *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
273}
274
275template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
277{
279
280 // Initialise helper variables
281 double heterotypic_boundary_length = 0.0;
282 double total_shared_edges_length = 0.0;
283 double num_heterotypic_pairs = 0.0;
284 double total_num_pairs = 0.0;
285
286 // Iterate over cells
287 for (typename AbstractCellPopulation<SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
288 cell_iter != pCellPopulation->End();
289 ++cell_iter)
290 {
291 // Store whether this cell is labelled
292 bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
293
294 unsigned elem_index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
295 unsigned num_nodes_in_elem = pCellPopulation->rGetMesh().GetElement(elem_index)->GetNumNodes();
296
297 // Iterate over nodes contained in the element corresponding to this cell
298 for (unsigned local_index=0; local_index<num_nodes_in_elem; local_index++)
299 {
300 // Get this node's von Neumann neighbours (not Moore neighbours, since they must share an edge)
301 unsigned global_index = pCellPopulation->rGetMesh().GetElement(elem_index)->GetNodeGlobalIndex(local_index);
302 std::set<unsigned> neighbour_node_indices = pCellPopulation->rGetMesh().GetVonNeumannNeighbouringNodeIndices(global_index);
303
304 // Iterate over these neighbours
305 for (std::set<unsigned>::iterator neighbour_iter = neighbour_node_indices.begin();
306 neighbour_iter != neighbour_node_indices.end();
307 ++neighbour_iter)
308 {
309 // Get the elements containing this neighbour
310 std::set<unsigned> neighbour_elem_indices = pCellPopulation->GetNode(*neighbour_iter)->rGetContainingElementIndices();
311
312 if (neighbour_elem_indices.size() == 1)
313 {
314 unsigned neigbour_elem_index = *(neighbour_elem_indices.begin());
315
316 if (neigbour_elem_index != elem_index)
317 {
318 // Edge is between two different elements
319 total_shared_edges_length += 1.0;
320
321 // Store whether the cell corresponding to this element index is labelled
322 CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(neigbour_elem_index);
323 bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
324
325 // If this cell is labelled and its neighbour is not, or vice versa...
326 if (cell_is_labelled != neighbour_is_labelled)
327 {
328 // ... then increment the fractional boundary length
329 heterotypic_boundary_length += 1.0;
330 }
331 }
332 }
333 else
334 {
335 // Original node is on boundary of mesh so don't include in the edge calculation.
336 }
337 }
338 }
339
340 std::set<unsigned> neighbour_node_indices = pCellPopulation->GetNeighbouringLocationIndices(*cell_iter);
341
342 // Iterate over these neighbours
343 for (std::set<unsigned>::iterator neighbour_iter = neighbour_node_indices.begin();
344 neighbour_iter != neighbour_node_indices.end();
345 ++neighbour_iter)
346 {
347 total_num_pairs += 1.0;
348
349 CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(*neighbour_iter);
350 bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
351
352 // If this cell is labelled and its neighbour is not, or vice versa...
353 if (cell_is_labelled != neighbour_is_labelled)
354 {
355 // ... then increment the fractional boundary length
356 num_heterotypic_pairs += 1.0;
357 }
358 }
359 }
360
361 // We have counted each cell-cell edge twice
362 heterotypic_boundary_length *= 0.5;
363 total_shared_edges_length *= 0.5;
364
365 // We have counted each pair of neighbouring cells twice
366 num_heterotypic_pairs *= 0.5;
367 total_num_pairs *= 0.5;
368
369 *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
370}
371
372template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
374{
375 // Make sure the cell population is updated
377 pCellPopulation->Update();
378
379 // Initialise helper variables
380 double heterotypic_boundary_length = 0.0;
381 double total_shared_edges_length = 0.0;
382 double num_heterotypic_pairs = 0.0;
383 double total_num_pairs = 0.0;
384
385 // Iterate over cells
386 for (typename AbstractCellPopulation<SPACE_DIM>::Iterator cell_iter = pCellPopulation->Begin();
387 cell_iter != pCellPopulation->End();
388 ++cell_iter)
389 {
390 // Store whether this cell is labelled
391 bool cell_is_labelled = cell_iter->template HasCellProperty<CellLabel>();
392
393 // Get the set of neighbouring element indices
394 unsigned elem_index = pCellPopulation->GetLocationIndexUsingCell(*cell_iter);
395 std::set<unsigned> neighbour_elem_indices = pCellPopulation->rGetMesh().GetNeighbouringElementIndices(elem_index);
396
397 // Iterate over these neighbours
398 for (std::set<unsigned>::iterator neighbour_iter = neighbour_elem_indices.begin();
399 neighbour_iter != neighbour_elem_indices.end();
400 ++neighbour_iter)
401 {
402 // Get the length of the edge shared with this neighbour
403 unsigned neighbour_index = *neighbour_iter;
404 double edge_length = pCellPopulation->rGetMesh().GetEdgeLength(elem_index, neighbour_index);
405
406 total_shared_edges_length += edge_length;
407 total_num_pairs += 1.0;
408
409 // Store whether this neighbour is labelled
410 CellPtr p_neighbour_cell = pCellPopulation->GetCellUsingLocationIndex(*neighbour_iter);
411 bool neighbour_is_labelled = p_neighbour_cell->template HasCellProperty<CellLabel>();
412
413 // If this cell is labelled and its neighbour is not, or vice versa...
414 if (cell_is_labelled != neighbour_is_labelled)
415 {
416 // ... then increment the fractional boundary length
417 heterotypic_boundary_length += edge_length;
418 num_heterotypic_pairs += 1.0;
419 }
420 }
421 }
422
423 // We have counted each cell-cell edge twice
424 heterotypic_boundary_length *= 0.5;
425 total_shared_edges_length *= 0.5;
426
427 // We have counted each pair of neighbouring cells twice
428 num_heterotypic_pairs *= 0.5;
429 total_num_pairs *= 0.5;
430
431 *this->mpOutStream << heterotypic_boundary_length << "\t" << total_shared_edges_length << "\t" << num_heterotypic_pairs << "\t" << total_num_pairs;
432}
433
434template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
436{
437 auto& r_mesh = pCellPopulation->rGetMesh();
438
439 // Initialise helper variables
440 double heterotypic_boundary_length = 0.0;
441 double total_shared_edges_length = 0.0;
442
443 // Keep track of which elements are interacting with which others
444 std::set<std::pair<unsigned, unsigned>> total_elem_pairs;
445 std::set<std::pair<unsigned, unsigned>> heterotypic_elem_pairs;
446
447 // Get the updated voronoi diagram for nodes
448 const bool update_diagram = true;
449 const boost::polygon::voronoi_diagram<double>& vd = r_mesh.rGetNodeLocationsVoronoiDiagram(update_diagram);
450
451 // Get the vector that allows us to map between node index and voronoi cell ID
452 const std::vector<unsigned>& node_idx_to_voronoi_cell_id_map = r_mesh.GetVoronoiCellIdsIndexedByNodeIndex();
453
454 for (const auto& p_node : r_mesh.rGetNodes())
455 {
456 const unsigned this_node_idx = p_node->GetIndex();
457 const unsigned this_elem_idx = *p_node->ContainingElementsBegin();
458
459 // Get the voronoi cell that corresponds to this node
460 const unsigned voronoi_cell_id = node_idx_to_voronoi_cell_id_map[this_node_idx];
461 const auto& voronoi_cell = vd.cells()[voronoi_cell_id];
462
463 // Iterate over the edges of this cell. Note that due to the halo region none of these edges should be infinite.
464 auto p_edge = voronoi_cell.incident_edge();
465 do
466 {
467 if (p_edge->is_finite())
468 {
469 // The global node index corresponding to a voronoi cell cell is encoded in its 'color' variable
470 const unsigned twin_node_idx = p_edge->twin()->cell()->color();
471 const unsigned twin_elem_idx = *r_mesh.GetNode(twin_node_idx)->ContainingElementsBegin();
472
473 // Check if the nodes are in different elements
474 if (this_elem_idx != twin_elem_idx)
475 {
476 // Add a pair to the set, which will record the total number of unique elem-elem interactions
477 const unsigned lo_idx = std::min(this_elem_idx, twin_elem_idx);
478 const unsigned hi_idx = std::max(this_elem_idx, twin_elem_idx);
479 total_elem_pairs.insert(std::make_pair(lo_idx, hi_idx));
480
481 const bool this_is_labelled = pCellPopulation->GetCellUsingLocationIndex(this_elem_idx)->template HasCellProperty<CellLabel>();
482 const bool twin_is_labelled = pCellPopulation->GetCellUsingLocationIndex(twin_elem_idx)->template HasCellProperty<CellLabel>();
483
484 const double edge_length = r_mesh.CalculateLengthOfVoronoiEdge(*p_edge);
485
486 total_shared_edges_length += edge_length;
487
488 if (this_is_labelled != twin_is_labelled)
489 {
490 heterotypic_boundary_length += edge_length;
491 heterotypic_elem_pairs.insert(std::make_pair(lo_idx, hi_idx));
492 }
493 }
494 }
495
496 p_edge = p_edge->next();
497 } while (p_edge != voronoi_cell.incident_edge());
498 }
499
500 // Every edge has been counted twice
501 heterotypic_boundary_length *= 0.5;
502 total_shared_edges_length *= 0.5;
503
504 *this->mpOutStream << heterotypic_boundary_length << '\t'
505 << total_shared_edges_length << '\t'
506 << heterotypic_elem_pairs.size() << '\t'
507 << total_elem_pairs.size();
508}
509
510// Explicit instantiation
517
519// Declare identifier for the serializer
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
virtual bool IsCellAttachedToLocationIndex(unsigned index)
unsigned GetLocationIndexUsingCell(CellPtr pCell)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
double GetDistanceBetweenNodes(unsigned indexA, unsigned indexB)
PottsMesh< DIM > & rGetMesh()
virtual void Visit(MeshBasedCellPopulation< ELEMENT_DIM, SPACE_DIM > *pCellPopulation)
ImmersedBoundaryMesh< DIM, DIM > & rGetMesh()
const boost::polygon::voronoi_diagram< double > & rGetNodeLocationsVoronoiDiagram(bool update=true)
double GetVoronoiEdgeLength(unsigned index1, unsigned index2)
std::set< unsigned > GetNeighbouringNodeIndices(unsigned index)
Node< DIM > * GetNode(unsigned index)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
NodesOnlyMesh< DIM > & rGetMesh()
void Update(bool hasHadBirthsOrDeaths=true)
std::set< unsigned > GetNeighbouringNodeIndices(unsigned index)
std::set< unsigned > & rGetContainingElementIndices()
Definition Node.cpp:300
double GetRadius()
Definition Node.cpp:248
std::set< unsigned > GetNeighbouringLocationIndices(CellPtr pCell)
Node< DIM > * GetNode(unsigned index)
MutableVertexMesh< DIM, DIM > & rGetMesh()
void Update(bool hasHadBirthsOrDeaths=true)
std::set< unsigned > GetNeighbouringElementIndices(unsigned elementIndex)
double GetEdgeLength(unsigned elementIndex1, unsigned elementIndex2)