Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
NodeBasedCellPopulation.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 "NodeBasedCellPopulation.hpp"
38#include "VtkMeshWriter.hpp"
39
40template<unsigned DIM>
42 std::vector<CellPtr>& rCells,
43 const std::vector<unsigned> locationIndices,
44 bool deleteMesh,
45 bool validate)
46 : AbstractCentreBasedCellPopulation<DIM>(rMesh, rCells, locationIndices),
47 mDeleteMesh(deleteMesh),
48 mUseVariableRadii(false),
49 mLoadBalanceMesh(false),
50 mLoadBalanceFrequency(100)
51{
52 mpNodesOnlyMesh = static_cast<NodesOnlyMesh<DIM>* >(&(this->mrMesh));
53
54 if (validate)
55 {
56 Validate();
57 }
58}
59
60template<unsigned DIM>
63 mDeleteMesh(true),
64 mUseVariableRadii(false), // will be set by serialize() method
65 mLoadBalanceMesh(false),
66 mLoadBalanceFrequency(100)
67{
68 mpNodesOnlyMesh = static_cast<NodesOnlyMesh<DIM>* >(&(this->mrMesh));
69}
70
71template<unsigned DIM>
73{
74 Clear();
75 if (mDeleteMesh)
76 {
77 delete &this->mrMesh;
78 }
79}
80
81template<unsigned DIM>
83{
84 return *mpNodesOnlyMesh;
85}
86
87template<unsigned DIM>
89{
90 return *mpNodesOnlyMesh;
91}
92
93template<unsigned DIM>
95{
96 // Note that this code does not yet work in parallel.
98
99 std::vector<Node<DIM>*> temp_nodes;
100
101 // Get the nodes of mpNodesOnlyMesh
102 for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = mpNodesOnlyMesh->GetNodeIteratorBegin();
103 node_iter != mpNodesOnlyMesh->GetNodeIteratorEnd();
104 ++node_iter)
105 {
106 temp_nodes.push_back(new Node<DIM>(node_iter->GetIndex(), node_iter->rGetLocation(), node_iter->IsBoundaryNode()));
107 }
108
109 return new MutableMesh<DIM,DIM>(temp_nodes);
110}
111
112template<unsigned DIM>
114{
115 mNodePairs.clear();
116}
117
118template<unsigned DIM>
120{
121 for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
122 node_iter != this->mrMesh.GetNodeIteratorEnd();
123 ++node_iter)
124 {
125 try
126 {
127 this->GetCellUsingLocationIndex(node_iter->GetIndex());
128 }
129 catch (Exception&)
130 {
131 EXCEPTION("At time " << SimulationTime::Instance()->GetTime() << ", Node " << node_iter->GetIndex() << " does not appear to have a cell associated with it");
132 }
133 }
134}
135
136template<unsigned DIM>
138{
139 return mpNodesOnlyMesh->GetNodeOrHaloNode(index);
140}
141
142template<unsigned DIM>
143void NodeBasedCellPopulation<DIM>::SetNode(unsigned nodeIndex, ChastePoint<DIM>& rNewLocation)
144{
145 mpNodesOnlyMesh->SetNode(nodeIndex, rNewLocation, false);
146}
147
148template<unsigned DIM>
149void NodeBasedCellPopulation<DIM>::Update(bool hasHadBirthsOrDeaths)
150{
151 UpdateCellProcessLocation();
152
153 mpNodesOnlyMesh->UpdateBoxCollection();
154
155 if (mLoadBalanceMesh)
156 {
157 if ((SimulationTime::Instance()->GetTimeStepsElapsed() % mLoadBalanceFrequency) == 0)
158 {
159 mpNodesOnlyMesh->LoadBalanceMesh();
160
161 UpdateCellProcessLocation();
162
163 mpNodesOnlyMesh->UpdateBoxCollection();
164 }
165 }
166
167 RefreshHaloCells();
168
169 mpNodesOnlyMesh->CalculateInteriorNodePairs(mNodePairs);
170
171 AddReceivedHaloCells();
172
173 mpNodesOnlyMesh->CalculateBoundaryNodePairs(mNodePairs);
174
175 /*
176 * Update cell radii based on CellData
177 */
178 if (mUseVariableRadii)
179 {
180 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
181 cell_iter != this->End();
182 ++cell_iter)
183 {
184 double cell_radius = cell_iter->GetCellData()->GetItem("Radius");
185 unsigned node_index = this->GetLocationIndexUsingCell(*cell_iter);
186 this->GetNode(node_index)->SetRadius(cell_radius);
187 }
188 }
189
190 // Make sure that everyone exits update together so that all asynchronous communications are complete.
191 PetscTools::Barrier("Update");
192}
193
194template<unsigned DIM>
196{
197 if (!map.IsIdentityMap())
198 {
199 UpdateParticlesAfterReMesh(map);
200
201 // Update the mappings between cells and location indices
203 std::map<Cell*, unsigned> old_map = this->mCellLocationMap;
204
205 // Remove any dead pointers from the maps (needed to avoid archiving errors)
206 this->mLocationCellMap.clear();
207 this->mCellLocationMap.clear();
208
209 for (std::list<CellPtr>::iterator it = this->mCells.begin();
210 it != this->mCells.end();
211 ++it)
212 {
213 unsigned old_node_index = old_map[(*it).get()];
214
215 // This shouldn't ever happen, as the cell vector only contains living cells
216 assert(!map.IsDeleted(old_node_index));
217
218 unsigned new_node_index = map.GetNewIndex(old_node_index);
219 this->SetCellUsingLocationIndex(new_node_index,*it);
220 }
221
222 this->Validate();
223 }
224}
225
226template<unsigned DIM>
228{
229 unsigned num_removed = 0;
230 for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
231 cell_iter != this->mCells.end();
232 )
233 {
234 if ((*cell_iter)->IsDead())
235 {
236 // Remove the node from the mesh
237 num_removed++;
238 unsigned location_index = mpNodesOnlyMesh->SolveNodeMapping(this->GetLocationIndexUsingCell((*cell_iter)));
239 mpNodesOnlyMesh->DeleteNodePriorToReMesh(location_index);
240
241 // Update mappings between cells and location indices
242 unsigned location_index_of_removed_node = this->GetLocationIndexUsingCell((*cell_iter));
243 this->RemoveCellUsingLocationIndex(location_index_of_removed_node, (*cell_iter));
244
245 // Update vector of cells
246 cell_iter = this->mCells.erase(cell_iter);
247 }
248 else
249 {
250 ++cell_iter;
251 }
252 }
253
254 return num_removed;
255}
256
257template<unsigned DIM>
259{
260 return mpNodesOnlyMesh->AddNode(pNewNode);
261}
262
263template<unsigned DIM>
265{
266 return mpNodesOnlyMesh->GetNumNodes();
267}
268
269template<unsigned DIM>
271{
272 std::map<unsigned, CellPtr>::iterator iter = mLocationHaloCellMap.find(index);
273 if (iter != mLocationHaloCellMap.end())
274 {
275 return iter->second;
276 }
277 else
278 {
280 }
281}
282
283template<unsigned DIM>
287
288template<unsigned DIM>
289std::vector< std::pair<Node<DIM>*, Node<DIM>* > >& NodeBasedCellPopulation<DIM>::rGetNodePairs()
290{
291 return mNodePairs;
292}
293
294template<unsigned DIM>
296{
297 *rParamsFile << "\t\t<MechanicsCutOffLength>" << mpNodesOnlyMesh->GetMaximumInteractionDistance() << "</MechanicsCutOffLength>\n";
298 *rParamsFile << "\t\t<UseVariableRadii>" << mUseVariableRadii << "</UseVariableRadii>\n";
299
300 // Call method on direct parent class
302}
303
304template<unsigned DIM>
306{
307 pPopulationWriter->Visit(this);
308}
309
310template<unsigned DIM>
312{
313 pPopulationCountWriter->Visit(this);
314}
315
316template<unsigned DIM>
318{
319 pPopulationEventWriter->Visit(this);
320}
321
322template<unsigned DIM>
323void NodeBasedCellPopulation<DIM>::AcceptCellWriter(boost::shared_ptr<AbstractCellWriter<DIM, DIM> > pCellWriter, CellPtr pCell)
324{
325 pCellWriter->VisitCell(pCell, this);
326}
327
328template<unsigned DIM>
330{
331 return mpNodesOnlyMesh->GetMaximumInteractionDistance();
332}
333
334template<unsigned DIM>
336{
337 return mUseVariableRadii;
338}
339
340template<unsigned DIM>
342{
343 mUseVariableRadii = useVariableRadii;
344}
345
346template<unsigned DIM>
348{
349 mLoadBalanceMesh = loadBalanceMesh;
350}
351
352template<unsigned DIM>
354{
355 mLoadBalanceFrequency = loadBalanceFrequency;
356}
357
358template<unsigned DIM>
359double NodeBasedCellPopulation<DIM>::GetWidth(const unsigned& rDimension)
360{
361 return mpNodesOnlyMesh->GetWidth(rDimension);
362}
363
364template<unsigned DIM>
366{
367 c_vector<double, DIM> local_size = AbstractCellPopulation<DIM, DIM>::GetSizeOfCellPopulation();
368 c_vector<double, DIM> global_size;
369
370 for (unsigned i=0; i<DIM; i++)
371 {
372 MPI_Allreduce(&local_size[i], &global_size[i], 1, MPI_DOUBLE, MPI_MAX, PetscTools::GetWorld());
373 }
374
375 return global_size;
376}
377
378template<unsigned DIM>
379std::set<unsigned> NodeBasedCellPopulation<DIM>::GetNodesWithinNeighbourhoodRadius(unsigned index, double neighbourhoodRadius)
380{
381 // Check neighbourhoodRadius is less than the interaction radius. If not you wont return all the correct nodes
382 if (neighbourhoodRadius > mpNodesOnlyMesh->GetMaximumInteractionDistance())
383 {
384 EXCEPTION("neighbourhoodRadius should be less than or equal to the the maximum interaction radius defined on the NodesOnlyMesh");
385 }
386
387 std::set<unsigned> neighbouring_node_indices;
388
389 // Get location
390 Node<DIM>* p_node_i = this->GetNode(index);
391 const c_vector<double, DIM>& r_node_i_location = p_node_i->rGetLocation();
392
393 // Check the mNodeNeighbours has been set up correctly
394 if (!(p_node_i->GetNeighboursSetUp()))
395 {
396 EXCEPTION("mNodeNeighbours not set up. Call Update() before GetNodesWithinNeighbourhoodRadius()");
397 }
398
399 // Get set of 'candidate' neighbours.
400 std::vector<unsigned>& near_nodes = p_node_i->rGetNeighbours();
401
402 // Find which ones are actually close
403 for (std::vector<unsigned>::iterator iter = near_nodes.begin();
404 iter != near_nodes.end();
405 ++iter)
406 {
407 // Be sure not to return the index itself.
408 if ((*iter) != index)
409 {
410 Node<DIM>* p_node_j = this->GetNode((*iter));
411
412 // Get the location of this node
413 const c_vector<double, DIM>& r_node_j_location = p_node_j->rGetLocation();
414
415 // Get the vector the two nodes (using GetVectorFromAtoB to catch periodicities etc.)
416 c_vector<double, DIM> node_to_node_vector = mpNodesOnlyMesh->GetVectorFromAtoB(r_node_j_location, r_node_i_location);
417
418 // Calculate the distance between the two nodes
419 double distance_between_nodes = norm_2(node_to_node_vector);
420
421 // If the cell j is within the neighbourhood of radius neighbourhoodRadius
422 // of cell i
423 if (distance_between_nodes <= neighbourhoodRadius)// + DBL_EPSILSON)
424 {
425 // ...then add this node index to the set of neighbouring node indices
426 neighbouring_node_indices.insert((*iter));
427 }
428 }
429 }
430
431 return neighbouring_node_indices;
432}
433
434template<unsigned DIM>
436{
437 std::set<unsigned> neighbouring_node_indices;
438
439 // Get location and radius of node
440 Node<DIM>* p_node_i = this->GetNode(index);
441 const c_vector<double, DIM>& r_node_i_location = p_node_i->rGetLocation();
442 double radius_of_cell_i = p_node_i->GetRadius();
443
444 // Check the mNodeNeighbours has been set up correctly
445 if (!(p_node_i->GetNeighboursSetUp()))
446 {
447 EXCEPTION("mNodeNeighbours not set up. Call Update() before GetNeighbouringNodeIndices()");
448 }
449
450 // Make sure that the max_interaction distance is smaller than or equal to the box collection size
451 if (!(radius_of_cell_i * 2.0 <= mpNodesOnlyMesh->GetMaximumInteractionDistance()))
452 {
453 EXCEPTION("mpNodesOnlyMesh::mMaxInteractionDistance is smaller than twice the radius of cell " << index << " (" << radius_of_cell_i << ") so interactions may be missed. Make the cut-off larger to avoid errors.");
454 }
455
456 // Get set of 'candidate' neighbours
457 std::vector<unsigned>& near_nodes = p_node_i->rGetNeighbours();
458
459 // Find which ones are actually close
460 for (std::vector<unsigned>::iterator iter = near_nodes.begin();
461 iter != near_nodes.end();
462 ++iter)
463 {
464 // Be sure not to return the index itself
465 if ((*iter) != index)
466 {
467 Node<DIM>* p_node_j = this->GetNode((*iter));
468
469 // Get the location of this node
470 const c_vector<double, DIM>& r_node_j_location = p_node_j->rGetLocation();
471
472 // Get the vector the two nodes (using GetVectorFromAtoB to catch periodicities etc.)
473 c_vector<double, DIM> node_to_node_vector = mpNodesOnlyMesh->GetVectorFromAtoB(r_node_j_location, r_node_i_location);
474
475 // Calculate the distance between the two nodes
476 double distance_between_nodes = norm_2(node_to_node_vector);
477
478 // Get the radius of the cell corresponding to this node
479 double radius_of_cell_j = p_node_j->GetRadius();
480
481 // If the cells are close enough to exert a force on each other...
482 double max_interaction_distance = radius_of_cell_i + radius_of_cell_j;
483
484 // Make sure that the max_interaction distance is smaller than or equal to the box collection size
485 if (!(max_interaction_distance <= mpNodesOnlyMesh->GetMaximumInteractionDistance()))
486 {
487 EXCEPTION("mpNodesOnlyMesh::mMaxInteractionDistance is smaller than the sum of radius of cell " << index << " (" << radius_of_cell_i << ") and cell " << (*iter) << " (" << radius_of_cell_j <<"). Make the cut-off larger to avoid errors.");
488 }
489 if (distance_between_nodes <= max_interaction_distance)// + DBL_EPSILSON) //Assumes that max_interaction_distance is of order 1
490 {
491 // ...then add this node index to the set of neighbouring node indices
492 neighbouring_node_indices.insert((*iter));
493 }
494 }
495 }
496
497 return neighbouring_node_indices;
498}
499
500template <unsigned DIM>
501double NodeBasedCellPopulation<DIM>::GetVolumeOfCell([[maybe_unused]] CellPtr pCell)
502{
503 // Not implemented or tested in 1D
504 if constexpr ((DIM == 2) || (DIM == 3))
505 {
506 // Get node index corresponding to this cell
507 unsigned node_index = this->GetLocationIndexUsingCell(pCell);
508 Node<DIM>* p_node = this->GetNode(node_index);
509
510 // Get cell radius
511 double cell_radius = p_node->GetRadius();
512
513 // Begin code to approximate cell volume
514 double averaged_cell_radius = 0.0;
515 unsigned num_cells = 0;
516
517 // Get the location of this node
518 const c_vector<double, DIM>& r_node_i_location = GetNode(node_index)->rGetLocation();
519
520 // Get the set of node indices corresponding to this cell's neighbours
521 std::set<unsigned> neighbouring_node_indices = GetNeighbouringNodeIndices(node_index);
522
523 // THe number of neighbours in equilibrium configuration, from sphere packing problem
524 unsigned num_neighbours_equil;
525 if (DIM == 2)
526 {
527 num_neighbours_equil = 6;
528 }
529 else // DIM == 3
530 {
531 num_neighbours_equil = 12;
532 }
533
534 // Loop over this set
535 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
536 iter != neighbouring_node_indices.end();
537 ++iter)
538 {
539 Node<DIM>* p_node_j = this->GetNode(*iter);
540
541 // Get the location of the neighbouring node
542 const c_vector<double, DIM>& r_node_j_location = p_node_j->rGetLocation();
543
544 double neighbouring_cell_radius = p_node_j->GetRadius();
545
546 // If this throws then you may not be considering all cell interactions use a larger cut off length
547 assert(cell_radius+neighbouring_cell_radius<mpNodesOnlyMesh->GetMaximumInteractionDistance());
548
549 // Calculate the distance between the two nodes and add to cell radius
550 double separation = norm_2(mpNodesOnlyMesh->GetVectorFromAtoB(r_node_j_location, r_node_i_location));
551
552 if (separation < cell_radius+neighbouring_cell_radius)
553 {
554 // The effective radius is the mid point of the overlap
555 averaged_cell_radius = averaged_cell_radius + cell_radius - (cell_radius+neighbouring_cell_radius-separation)/2.0;
556 num_cells++;
557 }
558 }
559 if (num_cells < num_neighbours_equil)
560 {
561 averaged_cell_radius += (num_neighbours_equil-num_cells)*cell_radius;
562
563 averaged_cell_radius /= num_neighbours_equil;
564 }
565 else
566 {
567 averaged_cell_radius /= num_cells;
568 }
569 assert(averaged_cell_radius < mpNodesOnlyMesh->GetMaximumInteractionDistance()/2.0);
570
571 cell_radius = averaged_cell_radius;
572
573 // End code to approximate cell volume
574
575 // Calculate cell volume from radius of cell
576 double cell_volume = 0.0;
577 if (DIM == 2)
578 {
579 cell_volume = M_PI*cell_radius*cell_radius;
580 }
581 else // DIM == 3
582 {
583 cell_volume = (4.0/3.0)*M_PI*cell_radius*cell_radius*cell_radius;
584 }
585
586 return cell_volume;
587 }
588 else
589 {
591 }
592}
593
594template<unsigned DIM>
596{
597#ifdef CHASTE_VTK
598 // Store the present time as a string
599 std::stringstream time;
601
602 // Make sure the nodes are ordered contiguously in memory
603 NodeMap map(1 + this->mpNodesOnlyMesh->GetMaximumNodeIndex());
604 this->mpNodesOnlyMesh->ReMesh(map);
605
606 // Create mesh writer for VTK output
607 VtkMeshWriter<DIM, DIM> mesh_writer(rDirectory, "results_"+time.str(), false);
608 mesh_writer.SetParallelFiles(*mpNodesOnlyMesh);
609
610 auto num_nodes = GetNumNodes();
611
612 // For each cell that this process owns, find the corresponding node index, which we only want to calculate once
613 std::vector<unsigned> node_indices_in_cell_order;
614 for (auto cell_iter = this->Begin(); cell_iter != this->End(); ++cell_iter)
615 {
616 // Get the node index corresponding to this cell
617 unsigned global_index = this->GetLocationIndexUsingCell(*cell_iter);
618 unsigned node_index = this->rGetMesh().SolveNodeMapping(global_index);
619
620 node_indices_in_cell_order.emplace_back(node_index);
621 }
622
623 // Iterate over any cell writers that are present. This is in a separate loop to below, because the writer loop
624 // needs to the the outer loop.
625 for (auto&& p_cell_writer : this->mCellWriters)
626 {
627 // Add any scalar data
628 if (p_cell_writer->GetOutputScalarData())
629 {
630 std::vector<double> vtk_cell_data(num_nodes);
631
632 unsigned loop_it = 0;
633 for (auto cell_iter = this->Begin(); cell_iter != this->End(); ++cell_iter, ++loop_it)
634 {
635 unsigned node_idx = node_indices_in_cell_order[loop_it];
636 vtk_cell_data[node_idx] = p_cell_writer->GetCellDataForVtkOutput(*cell_iter, this);
637 }
638
639 mesh_writer.AddPointData(p_cell_writer->GetVtkCellDataName(), vtk_cell_data);
640 }
641
642 // Add any vector data
643 if (p_cell_writer->GetOutputVectorData())
644 {
645 std::vector<c_vector<double, DIM>> vtk_cell_data(num_nodes);
646
647 unsigned loop_it = 0;
648 for (auto cell_iter = this->Begin(); cell_iter != this->End(); ++cell_iter, ++loop_it)
649 {
650 unsigned node_idx = node_indices_in_cell_order[loop_it];
651 vtk_cell_data[node_idx] = p_cell_writer->GetVectorCellDataForVtkOutput(*cell_iter, this);
652 }
653
654 mesh_writer.AddPointData(p_cell_writer->GetVtkVectorCellDataName(), vtk_cell_data);
655 }
656 }
657
658 // Process rank and cell data can be collected on a cell-by-cell basis, and both occur in the following loop
659 // We assume the first cell is representative of all cells
660 if (this->Begin() != this->End()) // some processes may own no cells, so can't do this->Begin()->GetCellData()
661 {
662 auto num_cell_data_items = this->Begin()->GetCellData()->GetNumItems();
663 std::vector<std::string> cell_data_names = this->Begin()->GetCellData()->GetKeys();
664 std::vector<std::vector<double>> cell_data(num_cell_data_items, std::vector<double>(num_nodes));
665
666 std::vector<double> rank(num_nodes);
667
668 unsigned loop_it = 0;
669 for (auto cell_iter = this->Begin(); cell_iter != this->End(); ++cell_iter, ++loop_it)
670 {
671 unsigned node_idx = node_indices_in_cell_order[loop_it];
672
673 for (unsigned cell_data_idx = 0; cell_data_idx < num_cell_data_items; ++cell_data_idx)
674 {
675 cell_data[cell_data_idx][node_idx] = cell_iter->GetCellData()->GetItem(cell_data_names[cell_data_idx]);
676 }
677
678 rank[node_idx] = (PetscTools::GetMyRank());
679 }
680
681 // Add point data to writers
682 mesh_writer.AddPointData("Process rank", rank);
683 for (unsigned cell_data_idx = 0; cell_data_idx < num_cell_data_items; ++cell_data_idx)
684 {
685 mesh_writer.AddPointData(cell_data_names[cell_data_idx], cell_data[cell_data_idx]);
686 }
687 }
688
689 mesh_writer.WriteFilesUsingMesh(*mpNodesOnlyMesh);
690
691 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
692 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
693 *(this->mpVtkMetaFile) << R"(" group="" part="0" file="results_)";
694 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
696 {
697 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
698 }
699 else
700 {
701 // Parallel vtu files .vtu -> .pvtu
702 *(this->mpVtkMetaFile) << ".pvtu\"/>\n";
703 }
704#endif //CHASTE_VTK
705}
706
707template<unsigned DIM>
708CellPtr NodeBasedCellPopulation<DIM>::AddCell(CellPtr pNewCell, CellPtr pParentCell)
709{
710 assert(pNewCell);
711
712 // Add new cell to population
713 CellPtr p_created_cell = AbstractCentreBasedCellPopulation<DIM>::AddCell(pNewCell, pParentCell);
714 assert(p_created_cell == pNewCell);
715
716 // Then set the new cell radius in the NodesOnlyMesh
717 unsigned parent_node_index = this->GetLocationIndexUsingCell(pParentCell);
718 Node<DIM>* p_parent_node = this->GetNode(parent_node_index);
719
720 unsigned new_node_index = this->GetLocationIndexUsingCell(p_created_cell);
721 Node<DIM>* p_new_node = this->GetNode(new_node_index);
722
723 p_new_node->SetRadius(p_parent_node->GetRadius());
724
725 // Return pointer to new cell
726 return p_created_cell;
727}
728
729template<unsigned DIM>
730void NodeBasedCellPopulation<DIM>::AddMovedCell(CellPtr pCell, boost::shared_ptr<Node<DIM> > pNode)
731{
732 // Create a new node
733 mpNodesOnlyMesh->AddMovedNode(pNode);
734
735 // Update cells vector
736 this->mCells.push_back(pCell);
737
738 // Update mappings between cells and location indices
739 this->AddCellUsingLocationIndex(pNode->GetIndex(), pCell);
740}
741
742template<unsigned DIM>
744{
745 mpNodesOnlyMesh->DeleteMovedNode(index);
746
747 // Update vector of cells
748 for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
749 cell_iter != this->mCells.end();
750 ++cell_iter)
751 {
752 if (this->GetLocationIndexUsingCell(*cell_iter) == index)
753 {
754 // Update mappings between cells and location indices
755 this->RemoveCellUsingLocationIndex(index, (*cell_iter));
756 cell_iter = this->mCells.erase(cell_iter);
757
758 break;
759 }
760 }
761}
762
763template<unsigned DIM>
765{
766 MPI_Status status;
767
769 {
770 boost::shared_ptr<std::vector<std::pair<CellPtr, Node<DIM>* > > > p_cells_right(&mCellsToSendRight, null_deleter());
771 mpCellsRecvRight = mRightCommunicator.SendRecvObject(p_cells_right, PetscTools::GetMyRank() + 1, mCellCommunicationTag, PetscTools::GetMyRank() + 1, mCellCommunicationTag, status);
772 }
774 {
775 boost::shared_ptr<std::vector<std::pair<CellPtr, Node<DIM>* > > > p_cells_left(&mCellsToSendLeft, null_deleter());
776 mpCellsRecvLeft = mLeftCommunicator.SendRecvObject(p_cells_left, PetscTools::GetMyRank() - 1, mCellCommunicationTag, PetscTools::GetMyRank() - 1, mCellCommunicationTag, status);
777 }
778 else if ( mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
779 {
780 boost::shared_ptr<std::vector<std::pair<CellPtr, Node<DIM>* > > > p_cells_left(&mCellsToSendLeft, null_deleter());
781 mpCellsRecvLeft = mLeftCommunicator.SendRecvObject(p_cells_left, PetscTools::GetNumProcs() - 1, mCellCommunicationTag, PetscTools::GetNumProcs() - 1, mCellCommunicationTag, status);
782}
783
784 // We need to leave this to the end (rather than as an else if for AmTopMost())
785 // otherwise there will be a cyclic send-receive and it will stall
786 if ( PetscTools::AmTopMost() && mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
787 {
788 boost::shared_ptr<std::vector<std::pair<CellPtr, Node<DIM>* > > > p_cells_right(&mCellsToSendRight, null_deleter());
789 mpCellsRecvRight = mRightCommunicator.SendRecvObject(p_cells_right, 0, mCellCommunicationTag, 0, mCellCommunicationTag, status);
790 }
791}
792
793// helper function for NonBlockingSendCellsToNeighbourProcess() to calculate the tag
794template<unsigned DIM>
795unsigned NodeBasedCellPopulation<DIM>::CalculateMessageTag(unsigned senderI, unsigned receiverJ)
796{
797 /* Old function was overloading integer type for nProcs > 12
798 unsigned tag = SmallPow(2u, 1+ PetscTools::GetMyRank() ) * SmallPow (3u, 1 + PetscTools::GetMyRank() + 1);
799 Instead we use a Cantor pairing function which produces lower paired values
800 See: https://en.wikipedia.org/wiki/Pairing_function */
801 unsigned tag = 0.5*((senderI+receiverJ)*(senderI+receiverJ+1) + 2*receiverJ);
802 assert(tag < UINT_MAX); //Just make sure doesnt hit UINT_MAX as old method did.
803 return tag;
804}
805
806template<unsigned DIM>
808{
810 {
811 boost::shared_ptr<std::vector<std::pair<CellPtr, Node<DIM>* > > > p_cells_right(&mCellsToSendRight, null_deleter());
812 unsigned tag = CalculateMessageTag( PetscTools::GetMyRank(), PetscTools::GetMyRank()+1 );
813 mRightCommunicator.ISendObject(p_cells_right, PetscTools::GetMyRank() + 1, tag);
814 }
816 {
817 unsigned tag = CalculateMessageTag( PetscTools::GetMyRank(), PetscTools::GetMyRank()-1 );
818 boost::shared_ptr<std::vector<std::pair<CellPtr, Node<DIM>* > > > p_cells_left(&mCellsToSendLeft, null_deleter());
819 mLeftCommunicator.ISendObject(p_cells_left, PetscTools::GetMyRank() - 1, tag);
820 }
821 else if ( mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
822 {
823 unsigned tag = CalculateMessageTag( PetscTools::GetMyRank(), PetscTools::GetNumProcs()-1 );
824 boost::shared_ptr<std::vector<std::pair<CellPtr, Node<DIM>* > > > p_cells_left(&mCellsToSendLeft, null_deleter());
825 mLeftCommunicator.ISendObject(p_cells_left, PetscTools::GetNumProcs()-1, tag);
826 }
827 if ( PetscTools::AmTopMost() && mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
828 {
829 unsigned tag = CalculateMessageTag( PetscTools::GetMyRank(), 0 );
830 boost::shared_ptr<std::vector<std::pair<CellPtr, Node<DIM>* > > > p_cells_right(&mCellsToSendRight, null_deleter());
831 mRightCommunicator.ISendObject(p_cells_right, 0, tag);
832 }
833
834 // Now post receives to start receiving data before returning.
836 {
837 unsigned tag = CalculateMessageTag( PetscTools::GetMyRank()+1, PetscTools::GetMyRank() );
838 mRightCommunicator.IRecvObject(PetscTools::GetMyRank() + 1, tag);
839 }
841 {
842 unsigned tag = CalculateMessageTag( PetscTools::GetMyRank()-1, PetscTools::GetMyRank() );
843 mLeftCommunicator.IRecvObject(PetscTools::GetMyRank() - 1, tag);
844 }
845 else if ( mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
846 {
847 unsigned tag = CalculateMessageTag( PetscTools::GetNumProcs()-1, PetscTools::GetMyRank() );
848 mLeftCommunicator.IRecvObject(PetscTools::GetNumProcs() - 1, tag);
849}
850 if ( PetscTools::AmTopMost() && mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
851 {
852 unsigned tag = CalculateMessageTag( 0, PetscTools::GetMyRank() );
853 mRightCommunicator.IRecvObject(0, tag);
854}
855}
856
857template<unsigned DIM>
859{
861 {
862 mpCellsRecvRight = mRightCommunicator.GetRecvObject();
863 }
864 if (!PetscTools::AmMaster() || mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
865 {
866 mpCellsRecvLeft = mLeftCommunicator.GetRecvObject();
867 }
868 // Periodicity across processors has to be in this set order
869 if ( PetscTools::AmTopMost() && mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
870 {
871 mpCellsRecvRight = mRightCommunicator.GetRecvObject();
872}
873}
874
875template<unsigned DIM>
876std::pair<CellPtr, Node<DIM>* > NodeBasedCellPopulation<DIM>::GetCellNodePair(unsigned nodeIndex)
877{
878 Node<DIM>* p_node = this->GetNode(nodeIndex);
879
880 CellPtr p_cell = this->GetCellUsingLocationIndex(nodeIndex);
881
882 std::pair<CellPtr, Node<DIM>* > new_pair(p_cell, p_node);
883
884 return new_pair;
885}
886
887template<unsigned DIM>
889{
890 std::pair<CellPtr, Node<DIM>* > pair = GetCellNodePair(nodeIndex);
891
892 mCellsToSendRight.push_back(pair);
893}
894
895template<unsigned DIM>
897{
898 std::pair<CellPtr, Node<DIM>* > pair = GetCellNodePair(nodeIndex);
899
900 mCellsToSendLeft.push_back(pair);
901}
902
903template<unsigned DIM>
905{
906 if (!PetscTools::AmMaster() || mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection() )
907 {
908 for (typename std::vector<std::pair<CellPtr, Node<DIM>* > >::iterator iter = mpCellsRecvLeft->begin();
909 iter != mpCellsRecvLeft->end();
910 ++iter)
911 {
912 // Make a shared pointer to the node to make sure it is correctly deleted.
913 boost::shared_ptr<Node<DIM> > p_node(iter->second);
914 AddMovedCell(iter->first, p_node);
915 }
916 }
917 if (!PetscTools::AmTopMost() || mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection())
918 {
919 for (typename std::vector<std::pair<CellPtr, Node<DIM>* > >::iterator iter = mpCellsRecvRight->begin();
920 iter != mpCellsRecvRight->end();
921 ++iter)
922 {
923 boost::shared_ptr<Node<DIM> > p_node(iter->second);
924 AddMovedCell(iter->first, p_node);
925 }
926 }
927}
928
929template<unsigned DIM>
931{
932 mpNodesOnlyMesh->ResizeBoxCollection();
933
934 mpNodesOnlyMesh->CalculateNodesOutsideLocalDomain();
935
936 std::vector<unsigned> nodes_to_send_right = mpNodesOnlyMesh->rGetNodesToSendRight();
937 AddCellsToSendRight(nodes_to_send_right);
938
939 std::vector<unsigned> nodes_to_send_left = mpNodesOnlyMesh->rGetNodesToSendLeft();
940 AddCellsToSendLeft(nodes_to_send_left);
941
942 // Post non-blocking send / receives so communication on both sides can start.
943 SendCellsToNeighbourProcesses();
944
945 // Post blocking receive calls that wait until communication complete.
946 //GetReceivedCells();
947
948 for (std::vector<unsigned>::iterator iter = nodes_to_send_right.begin();
949 iter != nodes_to_send_right.end();
950 ++iter)
951 {
952 DeleteMovedCell(*iter);
953 }
954
955 for (std::vector<unsigned>::iterator iter = nodes_to_send_left.begin();
956 iter != nodes_to_send_left.end();
957 ++iter)
958 {
959 DeleteMovedCell(*iter);
960 }
961
962 AddReceivedCells();
963
964 NodeMap map(1 + mpNodesOnlyMesh->GetMaximumNodeIndex());
965 mpNodesOnlyMesh->ReMesh(map);
966 UpdateMapsAfterRemesh(map);
967}
968
969template<unsigned DIM>
971{
972 mpNodesOnlyMesh->ClearHaloNodes();
973
974 mHaloCells.clear();
975 mHaloCellLocationMap.clear();
976 mLocationHaloCellMap.clear();
977
978 std::vector<unsigned> halos_to_send_right = mpNodesOnlyMesh->rGetHaloNodesToSendRight();
979 AddCellsToSendRight(halos_to_send_right);
980
981 std::vector<unsigned> halos_to_send_left = mpNodesOnlyMesh->rGetHaloNodesToSendLeft();
982 AddCellsToSendLeft(halos_to_send_left);
983
984 NonBlockingSendCellsToNeighbourProcesses();
985}
986
987template<unsigned DIM>
988void NodeBasedCellPopulation<DIM>::AddCellsToSendRight(std::vector<unsigned>& cellLocationIndices)
989{
990 mCellsToSendRight.clear();
991
992 for (unsigned i=0; i < cellLocationIndices.size(); i++)
993 {
994 AddNodeAndCellToSendRight(cellLocationIndices[i]);
995 }
996}
997
998template<unsigned DIM>
999void NodeBasedCellPopulation<DIM>::AddCellsToSendLeft(std::vector<unsigned>& cellLocationIndices)
1000{
1001 mCellsToSendLeft.clear();
1002
1003 for (unsigned i=0; i < cellLocationIndices.size(); i++)
1004 {
1005 AddNodeAndCellToSendLeft(cellLocationIndices[i]);
1006 }
1007}
1008
1009template<unsigned DIM>
1011{
1012 GetReceivedCells();
1013
1014 if (!PetscTools::AmMaster() || mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection())
1015 {
1016 for (typename std::vector<std::pair<CellPtr, Node<DIM>* > >::iterator iter = mpCellsRecvLeft->begin();
1017 iter != mpCellsRecvLeft->end();
1018 ++iter)
1019 {
1020 boost::shared_ptr<Node<DIM> > p_node(iter->second);
1021 AddHaloCell(iter->first, p_node);
1022
1023 }
1024 }
1025 if (!PetscTools::AmTopMost() || mpNodesOnlyMesh->GetIsPeriodicAcrossProcsFromBoxCollection())
1026 {
1027 for (typename std::vector<std::pair<CellPtr, Node<DIM>* > >::iterator iter = mpCellsRecvRight->begin();
1028 iter != mpCellsRecvRight->end();
1029 ++iter)
1030 {
1031 boost::shared_ptr<Node<DIM> > p_node(iter->second);
1032 AddHaloCell(iter->first, p_node);
1033 }
1034 }
1035
1036 mpNodesOnlyMesh->AddHaloNodesToBoxes();
1037}
1038
1039template<unsigned DIM>
1040void NodeBasedCellPopulation<DIM>::AddHaloCell(CellPtr pCell, boost::shared_ptr<Node<DIM> > pNode)
1041{
1042 mHaloCells.push_back(pCell);
1043
1044 mpNodesOnlyMesh->AddHaloNode(pNode);
1045
1046 mHaloCellLocationMap[pCell] = pNode->GetIndex();
1047
1048 mLocationHaloCellMap[pNode->GetIndex()] = pCell;
1049}
1050
1051// Explicit instantiation
1052template class NodeBasedCellPopulation<1>;
1053template class NodeBasedCellPopulation<2>;
1054template class NodeBasedCellPopulation<3>;
1055
1056// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define NEVER_REACHED
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
AbstractMesh< ELEMENT_DIM, SPACE_DIM > & mrMesh
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
c_vector< double, SPACE_DIM > GetSizeOfCellPopulation()
virtual void OutputCellPopulationParameters(out_stream &rParamsFile)
CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell=CellPtr())
c_vector< double, DIM > GetSizeOfCellPopulation()
unsigned CalculateMessageTag(unsigned senderI, unsigned receiverJ)
void SetNode(unsigned nodeIndex, ChastePoint< DIM > &rNewLocation)
void AddCellsToSendRight(std::vector< unsigned > &cellLocationIndices)
void AddHaloCell(CellPtr pCell, boost::shared_ptr< Node< DIM > > pNode)
Node< DIM > * GetNode(unsigned index)
virtual void AcceptPopulationCountWriter(boost::shared_ptr< AbstractCellPopulationCountWriter< DIM, DIM > > pPopulationCountWriter)
void OutputCellPopulationParameters(out_stream &rParamsFile)
void AddNodeAndCellToSendLeft(unsigned nodeIndex)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
NodesOnlyMesh< DIM > & rGetMesh()
NodeBasedCellPopulation(NodesOnlyMesh< DIM > &rMesh, std::vector< CellPtr > &rCells, const std::vector< unsigned > locationIndices=std::vector< unsigned >(), bool deleteMesh=false, bool validate=true)
void SetUseVariableRadii(bool useVariableRadii=true)
std::pair< CellPtr, Node< DIM > * > GetCellNodePair(unsigned nodeIndex)
void SetLoadBalanceMesh(bool loadBalanceMesh)
virtual void WriteVtkResultsToFile(const std::string &rDirectory)
void AddCellsToSendLeft(std::vector< unsigned > &cellLocationIndices)
virtual void AcceptCellWriter(boost::shared_ptr< AbstractCellWriter< DIM, DIM > > pCellWriter, CellPtr pCell)
unsigned AddNode(Node< DIM > *pNewNode)
std::set< unsigned > GetNodesWithinNeighbourhoodRadius(unsigned index, double neighbourhoodRadius)
double GetWidth(const unsigned &rDimension)
virtual TetrahedralMesh< DIM, DIM > * GetTetrahedralMeshForPdeModifier()
NodesOnlyMesh< DIM > * mpNodesOnlyMesh
void SetLoadBalanceFrequency(unsigned loadBalanceFrequency)
double GetVolumeOfCell(CellPtr pCell)
void Update(bool hasHadBirthsOrDeaths=true)
virtual CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell)
void AddNodeAndCellToSendRight(unsigned nodeIndex)
void AddMovedCell(CellPtr pCell, boost::shared_ptr< Node< DIM > > pNode)
std::set< unsigned > GetNeighbouringNodeIndices(unsigned index)
virtual void UpdateParticlesAfterReMesh(NodeMap &rMap)
virtual void AcceptPopulationEventWriter(boost::shared_ptr< AbstractCellPopulationEventWriter< DIM, DIM > > pPopulationEventWriter)
std::vector< std::pair< Node< DIM > *, Node< DIM > * > > & rGetNodePairs()
virtual void AcceptPopulationWriter(boost::shared_ptr< AbstractCellPopulationWriter< DIM, DIM > > pPopulationWriter)
unsigned GetNewIndex(unsigned oldIndex) const
Definition NodeMap.cpp:87
bool IsIdentityMap()
Definition NodeMap.cpp:100
bool IsDeleted(unsigned index)
Definition NodeMap.cpp:82
Definition Node.hpp:59
void SetRadius(double radius)
Definition Node.cpp:256
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition Node.cpp:139
bool GetNeighboursSetUp()
Definition Node.cpp:368
std::vector< unsigned > & rGetNeighbours()
Definition Node.cpp:376
double GetRadius()
Definition Node.cpp:248
static MPI_Comm GetWorld()
static bool AmMaster()
static bool AmTopMost()
static void Barrier(const std::string callerId="")
static bool IsSequential()
static unsigned GetMyRank()
static unsigned GetNumProcs()
static SimulationTime * Instance()
unsigned GetTimeStepsElapsed() const
void AddPointData(std::string name, std::vector< double > data)
void SetParallelFiles(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)