Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
NodeBasedCellPopulationWithParticles.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 "NodeBasedCellPopulationWithParticles.hpp"
37#include "VtkMeshWriter.hpp"
38#include "CellAgesWriter.hpp"
39#include "CellAncestorWriter.hpp"
40#include "CellProliferativePhasesWriter.hpp"
41#include "CellVolumesWriter.hpp"
42#include "CellMutationStatesCountWriter.hpp"
43
44template<unsigned DIM>
46 std::vector<CellPtr>& rCells,
47 const std::vector<unsigned> locationIndices,
48 bool deleteMesh)
49 : NodeBasedCellPopulation<DIM>(rMesh, rCells, locationIndices, deleteMesh, false)
50{
52 if (!locationIndices.empty())
53 {
54 // Create a set of node indices corresponding to particles
55 std::set<unsigned> node_indices;
56 std::set<unsigned> location_indices;
57 std::set<unsigned> particle_indices;
58
59 for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = rMesh.GetNodeIteratorBegin();
60 node_iter != rMesh.GetNodeIteratorEnd();
61 ++node_iter)
62 {
63 node_indices.insert(node_iter->GetIndex());
64 }
65 for (unsigned i=0; i<locationIndices.size(); i++)
66 {
67 location_indices.insert(locationIndices[i]);
68 }
69
70 std::set_difference(node_indices.begin(), node_indices.end(),
71 location_indices.begin(), location_indices.end(),
72 std::inserter(particle_indices, particle_indices.begin()));
73
74 // This method finishes and then calls Validate()
75 SetParticles(particle_indices);
76 }
77 else
78 {
79 for (typename NodesOnlyMesh<DIM>::NodeIterator node_iter = rMesh.GetNodeIteratorBegin();
80 node_iter != rMesh.GetNodeIteratorEnd();
81 ++node_iter)
82 {
83 (*node_iter).SetIsParticle(false);
84 }
86 }
87}
88
89template<unsigned DIM>
94
95template<unsigned DIM>
97{
98 return this->GetNode(index)->IsParticle();
99}
100
101template<unsigned DIM>
103{
104 std::set<unsigned> particle_indices;
105
106 for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
107 node_iter != this->mrMesh.GetNodeIteratorEnd();
108 ++node_iter)
109 {
110 if (node_iter->IsParticle())
111 {
112 particle_indices.insert(node_iter->GetIndex());
113 }
114 }
115
116 return particle_indices;
117}
118
119template<unsigned DIM>
120void NodeBasedCellPopulationWithParticles<DIM>::SetParticles(const std::set<unsigned>& rParticleIndices)
121{
122 for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
123 node_iter != this->mrMesh.GetNodeIteratorEnd();
124 ++node_iter)
125 {
126 if (rParticleIndices.find(node_iter->GetIndex()) != rParticleIndices.end())
127 {
128 node_iter->SetIsParticle(true);
129 }
130 }
132}
133
134template<unsigned DIM>
138
139template<unsigned DIM>
140CellPtr NodeBasedCellPopulationWithParticles<DIM>::AddCell(CellPtr pNewCell, CellPtr pParentCell)
141{
142 assert(pNewCell);
143
144 // Add new cell to population
145 CellPtr p_created_cell = AbstractCentreBasedCellPopulation<DIM>::AddCell(pNewCell, pParentCell);
146 assert(p_created_cell == pNewCell);
147
148 // Then set the new cell radius in the NodesOnlyMesh
149 unsigned node_index = this->GetLocationIndexUsingCell(p_created_cell);
150 this->GetNode(node_index)->SetRadius(0.5);
151 this->GetNode(node_index)->SetIsParticle(false);
152
153 // Return pointer to new cell
154 return p_created_cell;
155}
156
157template<unsigned DIM>
159{
160 std::map<unsigned, bool> validated_nodes;
161 for (typename AbstractMesh<DIM, DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
162 node_iter != this->mrMesh.GetNodeIteratorEnd();
163 ++node_iter)
164 {
165 validated_nodes[node_iter->GetIndex()] = node_iter->IsParticle();
166 }
167
168 // Look through all of the cells and record what node they are associated with.
169 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
170 {
171 unsigned node_index = this->GetLocationIndexUsingCell((*cell_iter));
172
173 // If the node attached to this cell is labelled as a particle, then throw an error
174 if (this->GetNode(node_index)->IsParticle())
175 {
176 EXCEPTION("Node " << node_index << " is labelled as a particle and has a cell attached");
177 }
178 validated_nodes[node_index] = true;
179 }
180
181 for (std::map<unsigned, bool>::iterator map_iter = validated_nodes.begin();
182 map_iter != validated_nodes.end();
183 map_iter++)
184 {
185 if (!map_iter->second)
186 {
187 EXCEPTION("Node " << map_iter->first << " does not appear to be a particle or has a cell associated with it");
188 }
189 }
190}
191
192template<unsigned DIM>
194{
195 for (typename AbstractMesh<DIM, DIM>::NodeIterator node_iter = this->rGetMesh().GetNodeIteratorBegin();
196 node_iter != this->rGetMesh().GetNodeIteratorEnd();
197 ++node_iter)
198 {
199 // If it isn't a particle then there might be cell writers attached
200 if (! this->IsParticle(node_iter->GetIndex()))
201 {
202 for (typename std::vector<boost::shared_ptr<AbstractCellWriter<DIM, DIM> > >::iterator cell_writer_iter = this->mCellWriters.begin();
203 cell_writer_iter != this->mCellWriters.end();
204 ++cell_writer_iter)
205 {
206 CellPtr cell_from_node = this->GetCellUsingLocationIndex(node_iter->GetIndex());
207 this->AcceptCellWriter(*cell_writer_iter, cell_from_node);
208 }
209 }
210 }
211}
212
213template<unsigned DIM>
215{
216#ifdef CHASTE_VTK
217 // Store the present time as a string
218 std::stringstream time;
220
221 // Make sure the nodes are ordered contiguously in memory
222 NodeMap map(1 + this->mpNodesOnlyMesh->GetMaximumNodeIndex());
223 this->mpNodesOnlyMesh->ReMesh(map);
224
225 // Store the number of cells for which to output data to VTK
226 unsigned num_nodes = this->GetNumNodes();
227 std::vector<double> rank(num_nodes);
228 std::vector<double> particles(num_nodes);
229
230 unsigned num_cell_data_items = 0;
231 std::vector<std::string> cell_data_names;
232
233 // We assume that the first cell is representative of all cells
234 if (num_nodes > 0)
235 {
236 num_cell_data_items = this->Begin()->GetCellData()->GetNumItems();
237 cell_data_names = this->Begin()->GetCellData()->GetKeys();
238 }
239
240 std::vector<std::vector<double> > cell_data;
241 for (unsigned var=0; var<num_cell_data_items; var++)
242 {
243 std::vector<double> cell_data_var(num_nodes);
244 cell_data.push_back(cell_data_var);
245 }
246
247 // Create mesh writer for VTK output
248 VtkMeshWriter<DIM, DIM> mesh_writer(rDirectory, "results_"+time.str(), false);
249 mesh_writer.SetParallelFiles(*(this->mpNodesOnlyMesh));
250
251 // Iterate over any cell writers that are present
252 for (typename std::vector<boost::shared_ptr<AbstractCellWriter<DIM, DIM> > >::iterator cell_writer_iter = this->mCellWriters.begin();
253 cell_writer_iter != this->mCellWriters.end();
254 ++cell_writer_iter)
255 {
256 // Create vector to store VTK cell data
257 std::vector<double> vtk_cell_data(num_nodes);
258
259 // Loop over nodes
260 for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
261 node_iter != this->mrMesh.GetNodeIteratorEnd();
262 ++node_iter)
263 {
264 unsigned node_index = node_iter->GetIndex();
265
266 // If this node is a particle (not a cell), then we set the 'dummy' VTK cell data for this to be -2.0...
267 if (this->IsParticle(node_index))
268 {
269 vtk_cell_data[node_index] = -2.0;
270 }
271 else
272 {
273 // ...otherwise we populate the vector of VTK cell data as usual
274 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
275 vtk_cell_data[node_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(p_cell, this);
276 }
277 }
278
279 mesh_writer.AddPointData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
280 }
281
282 // Loop over cells
283 for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
284 cell_iter != this->End();
285 ++cell_iter)
286 {
287 // Get the node index corresponding to this cell
288 unsigned global_index = this->GetLocationIndexUsingCell(*cell_iter);
289 unsigned node_index = this->rGetMesh().SolveNodeMapping(global_index);
290
291 for (unsigned var=0; var<num_cell_data_items; var++)
292 {
293 cell_data[var][node_index] = cell_iter->GetCellData()->GetItem(cell_data_names[var]);
294 }
295
296 rank[node_index] = (PetscTools::GetMyRank());
297 }
298
299 mesh_writer.AddPointData("Process rank", rank);
300
301 // Loop over nodes
302 for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
303 node_iter != this->mrMesh.GetNodeIteratorEnd();
304 ++node_iter)
305 {
306 unsigned node_index = node_iter->GetIndex();
307 particles[node_index] = (double) (this->IsParticle(node_index));
308 }
309
310 mesh_writer.AddPointData("Non-particles", particles);
311
312 if (num_cell_data_items > 0)
313 {
314 for (unsigned var=0; var<cell_data.size(); var++)
315 {
316 mesh_writer.AddPointData(cell_data_names[var], cell_data[var]);
317 }
318 }
319
320 mesh_writer.WriteFilesUsingMesh(*(this->mpNodesOnlyMesh));
321
322 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
323 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
324 *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
325 *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
327 {
328 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
329 }
330/* {
331 // Parallel vtu files .vtu -> .pvtu
332 *(this->mpVtkMetaFile) << ".pvtu\"/>\n";
333 }*/
334#endif //CHASTE_VTK
335}
336
337template<unsigned DIM>
339{
340 // Call method on direct parent class
342}
343
344// Explicit instantiation
348
349// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define EXCEPT_IF_NOT(test)
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell=CellPtr())
NodeIterator GetNodeIteratorEnd()
NodeIterator GetNodeIteratorBegin(bool skipDeletedNodes=true)
CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell)
virtual void WriteVtkResultsToFile(const std::string &rDirectory)
void SetParticles(const std::set< unsigned > &rParticleIndices)
NodeBasedCellPopulationWithParticles(NodesOnlyMesh< DIM > &rMesh, std::vector< CellPtr > &rCells, const std::vector< unsigned > locationIndices=std::vector< unsigned >(), bool deleteMesh=false)
void OutputCellPopulationParameters(out_stream &rParamsFile)
static bool IsSequential()
static unsigned GetMyRank()
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)