Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
MeshBasedCellPopulationWithGhostNodes.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 "MeshBasedCellPopulationWithGhostNodes.hpp"
37#include "CellLocationIndexWriter.hpp"
38#include "VtkMeshWriter.hpp"
39
40template<unsigned DIM>
43 std::vector<CellPtr>& rCells,
44 const std::vector<unsigned> locationIndices,
45 bool deleteMesh,
46 double ghostCellSpringStiffness,
47 double ghostGhostSpringStiffness,
48 double ghostSpringRestLength)
49 : MeshBasedCellPopulation<DIM,DIM>(rMesh, rCells, locationIndices, deleteMesh, false), // do not call the base class Validate()
50 mGhostCellSpringStiffness(ghostCellSpringStiffness),
51 mGhostGhostSpringStiffness(ghostGhostSpringStiffness),
52 mGhostSpringRestLength(ghostSpringRestLength)
53{
54 if (!locationIndices.empty())
55 {
56 // Create a set of node indices corresponding to ghost nodes
57 std::set<unsigned> node_indices;
58 std::set<unsigned> location_indices;
59 std::set<unsigned> ghost_node_indices;
60
61 for (unsigned i=0; i<this->GetNumNodes(); i++)
62 {
63 node_indices.insert(this->GetNode(i)->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(ghost_node_indices, ghost_node_indices.begin()));
73
74 // This method finishes and then calls Validate()
75 SetGhostNodes(ghost_node_indices);
76 }
77 else
78 {
79 this->mIsGhostNode = std::vector<bool>(this->GetNumNodes(), false);
80 Validate();
81 }
82}
83
84template<unsigned DIM>
86 double ghostCellSpringStiffness,
87 double ghostGhostSpringStiffness,
88 double ghostSpringRestLength)
89 : MeshBasedCellPopulation<DIM,DIM>(rMesh),
90 mGhostCellSpringStiffness(ghostCellSpringStiffness),
91 mGhostGhostSpringStiffness(ghostGhostSpringStiffness),
92 mGhostSpringRestLength(ghostSpringRestLength)
93{
94}
95
96template<unsigned DIM>
100
101template<unsigned DIM>
103{
104 EXCEPTION("Currently can't solve PDEs on meshes with ghost nodes");
105 return static_cast<TetrahedralMesh<DIM, DIM>*>(&(this->mrMesh));
106}
107
108template<unsigned DIM>
110{
111 return this->mIsGhostNode;
112}
113
114template<unsigned DIM>
116{
117 return this->mIsGhostNode[index];
118}
119
120template<unsigned DIM>
122{
123 std::set<unsigned> ghost_node_indices;
124 for (unsigned i=0; i<this->mIsGhostNode.size(); i++)
125 {
126 if (this->mIsGhostNode[i])
127 {
128 ghost_node_indices.insert(i);
129 }
130 }
131 return ghost_node_indices;
132}
133
134template<unsigned DIM>
135void MeshBasedCellPopulationWithGhostNodes<DIM>::SetGhostNodes(const std::set<unsigned>& rGhostNodeIndices)
136{
137 // Reinitialise all entries of mIsGhostNode to false
138 this->mIsGhostNode = std::vector<bool>(this->mrMesh.GetNumNodes(), false);
139
140 // Update mIsGhostNode
141 for (std::set<unsigned>::iterator iter=rGhostNodeIndices.begin(); iter!=rGhostNodeIndices.end(); ++iter)
142 {
143 this->mIsGhostNode[*iter] = true;
144 }
145
146 Validate();
147}
148
149template<unsigned DIM>
150c_vector<double, DIM> MeshBasedCellPopulationWithGhostNodes<DIM>::CalculateForceBetweenGhostNodes(const unsigned& rNodeAGlobalIndex, const unsigned& rNodeBGlobalIndex)
151{
152 assert(rNodeAGlobalIndex != rNodeBGlobalIndex);
153 c_vector<double, DIM> unit_difference;
154 const c_vector<double, DIM>& r_node_a_location = this->GetNode(rNodeAGlobalIndex)->rGetLocation();
155 const c_vector<double, DIM>& r_node_b_location = this->GetNode(rNodeBGlobalIndex)->rGetLocation();
156
157 // There is reason not to subtract one position from the other (cylindrical meshes)
158 unit_difference = this->mrMesh.GetVectorFromAtoB(r_node_a_location, r_node_b_location);
159
160 double distance_between_nodes = norm_2(unit_difference);
161 unit_difference /= distance_between_nodes;
162
163 double rest_length = 1.0; // TODO this could also be a parameter.
164 double spring_stiffness = mGhostCellSpringStiffness;
165 if (this->mIsGhostNode[rNodeAGlobalIndex] && this->mIsGhostNode[rNodeBGlobalIndex])
166 {
167 rest_length = mGhostSpringRestLength;
168 spring_stiffness = mGhostGhostSpringStiffness;
169 }
170
171 return spring_stiffness * unit_difference * (distance_between_nodes - rest_length);
172}
173
174template<unsigned DIM>
175CellPtr MeshBasedCellPopulationWithGhostNodes<DIM>::AddCell(CellPtr pNewCell, CellPtr pParentCell)
176{
177 // Add new cell to population
178 CellPtr p_created_cell = MeshBasedCellPopulation<DIM,DIM>::AddCell(pNewCell, pParentCell);
179 assert(p_created_cell == pNewCell);
180
181 // Update size of mIsGhostNode if necessary
182 unsigned new_node_index = this->GetLocationIndexUsingCell(p_created_cell);
183
184 if (this->GetNumNodes() > this->mIsGhostNode.size())
185 {
186 this->mIsGhostNode.resize(this->GetNumNodes());
187 this->mIsGhostNode[new_node_index] = false;
188 }
189
190 // Return pointer to new cell
191 return p_created_cell;
192}
193
194template<unsigned DIM>
196{
197 // Get a list of all the nodes that are ghosts
198 std::vector<bool> validated_node = mIsGhostNode;
199 assert(mIsGhostNode.size()==this->GetNumNodes());
200
201 // Look through all of the cells and record what node they are associated with.
202 for (typename AbstractCellPopulation<DIM,DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
203 {
204 unsigned node_index = this->GetLocationIndexUsingCell((*cell_iter));
205
206 // If the node attached to this cell is labelled as a ghost node, then throw an error
207 if (mIsGhostNode[node_index])
208 {
209 EXCEPTION("Node " << node_index << " is labelled as a ghost node and has a cell attached");
210 }
211 validated_node[node_index] = true;
212 }
213
214 for (unsigned i=0; i<validated_node.size(); i++)
215 {
216 if (!validated_node[i])
217 {
218 EXCEPTION("Node " << i << " does not appear to be a ghost node or have a cell associated with it");
219 }
220 }
221}
222
223template<unsigned DIM>
225{
226 assert(mIsGhostNode[nodeIndex]);
227
228 static_cast<MutableMesh<DIM,DIM>&>((this->mrMesh)).DeleteNodePriorToReMesh(nodeIndex);
229}
230
231template<unsigned DIM>
233{
234 // Copy mIsGhostNode to a temporary vector
235 std::vector<bool> ghost_nodes_before_remesh = mIsGhostNode;
236
237 // Reinitialise mIsGhostNode
238 mIsGhostNode.clear();
239 mIsGhostNode.resize(this->GetNumNodes());
240
241 // Update mIsGhostNode using the node map
242 for (unsigned old_index=0; old_index<rMap.GetSize(); old_index++)
243 {
244 if (!rMap.IsDeleted(old_index))
245 {
246 unsigned new_index = rMap.GetNewIndex(old_index);
247 mIsGhostNode[new_index] = ghost_nodes_before_remesh[old_index];
248 }
249 }
250}
251
252template<unsigned DIM>
254{
255 unsigned node_index = this->GetLocationIndexUsingCell(pCell);
256 std::set<unsigned> neighbour_indices = this->GetNeighbouringNodeIndices(node_index);
257
258 // Remove ghost nodes from the neighbour indices
259 for (std::set<unsigned>::iterator iter = neighbour_indices.begin();
260 iter != neighbour_indices.end();)
261 {
262 if (this->IsGhostNode(*iter))
263 {
264 neighbour_indices.erase(iter++);
265 }
266 else
267 {
268 ++iter;
269 }
270 }
271
272 return neighbour_indices;
273}
274
275template<unsigned DIM>
277{
278 for (typename AbstractMesh<DIM, DIM>::NodeIterator node_iter = this->rGetMesh().GetNodeIteratorBegin();
279 node_iter != this->rGetMesh().GetNodeIteratorEnd();
280 ++node_iter)
281 {
282 // If it isn't a ghost node then there might be cell writers attached
283 if (! this->IsGhostNode(node_iter->GetIndex()))
284 {
285 for (typename std::vector<boost::shared_ptr<AbstractCellWriter<DIM, DIM> > >::iterator cell_writer_iter = this->mCellWriters.begin();
286 cell_writer_iter != this->mCellWriters.end();
287 ++cell_writer_iter)
288 {
289 CellPtr cell_from_node = this->GetCellUsingLocationIndex(node_iter->GetIndex());
290 this->AcceptCellWriter(*cell_writer_iter, cell_from_node);
291 }
292 }
293 }
294}
295
296template<unsigned DIM>
298{
299 // Initialise vector of forces on ghost nodes
300 std::vector<c_vector<double, DIM> > drdt(this->GetNumNodes());
301 for (unsigned i=0; i<drdt.size(); i++)
302 {
303 drdt[i] = zero_vector<double>(DIM);
304 }
305
306 // Calculate forces on ghost nodes
307 for (typename MutableMesh<DIM, DIM>::EdgeIterator edge_iterator = static_cast<MutableMesh<DIM, DIM>&>((this->mrMesh)).EdgesBegin();
308 edge_iterator != static_cast<MutableMesh<DIM, DIM>&>((this->mrMesh)).EdgesEnd();
309 ++edge_iterator)
310 {
311 unsigned nodeA_global_index = edge_iterator.GetNodeA()->GetIndex();
312 unsigned nodeB_global_index = edge_iterator.GetNodeB()->GetIndex();
313
314 c_vector<double, DIM> force = CalculateForceBetweenGhostNodes(nodeA_global_index, nodeB_global_index);
315
316 if (!this->mIsGhostNode[nodeA_global_index])
317 {
318 drdt[nodeB_global_index] -= force;
319 }
320 else
321 {
322 drdt[nodeA_global_index] += force;
323
324 if (this->mIsGhostNode[nodeB_global_index])
325 {
326 drdt[nodeB_global_index] -= force;
327 }
328 }
329 }
330
331 for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
332 node_iter != this->mrMesh.GetNodeIteratorEnd();
333 ++node_iter)
334 {
335 unsigned node_index = node_iter->GetIndex();
336 if (this->mIsGhostNode[node_index])
337 {
338 node_iter->ClearAppliedForce();
339 node_iter->AddAppliedForceContribution(drdt[node_index]);
340 }
341 }
342}
343
344template<unsigned DIM>
346{
347 if (this->mOutputResultsForChasteVisualizer)
348 {
349 if (!this-> template HasWriter<CellLocationIndexWriter>())
350 {
351 this-> template AddCellWriter<CellLocationIndexWriter>();
352 }
353 }
354
356}
357
358template<unsigned DIM>
360{
361#ifdef CHASTE_VTK
362 // Store the present time as a string
363 unsigned num_timesteps = SimulationTime::Instance()->GetTimeStepsElapsed();
364 std::stringstream time;
365 time << num_timesteps;
366
367 if (this->mWriteVtkAsPoints)
368 {
369 // Create mesh writer for VTK output
370 VtkMeshWriter<DIM, DIM> mesh_writer(rDirectory, "mesh_results_"+time.str(), false);
371
372 // Iterate over any cell writers that are present
373 unsigned num_vtk_cells = this->rGetMesh().GetNumNodes();
374 for (typename std::vector<boost::shared_ptr<AbstractCellWriter<DIM, DIM> > >::iterator cell_writer_iter = this->mCellWriters.begin();
375 cell_writer_iter != this->mCellWriters.end();
376 ++cell_writer_iter)
377 {
378 // Create vector to store VTK cell data
379 std::vector<double> vtk_cell_data(num_vtk_cells);
380
381 // Loop over nodes of mesh
382 for (typename AbstractMesh<DIM, DIM>::NodeIterator node_iter = this->rGetMesh().GetNodeIteratorBegin();
383 node_iter != this->rGetMesh().GetNodeIteratorEnd();
384 ++node_iter)
385 {
386 // Get the indices of this node
387 unsigned node_index = node_iter->GetIndex();
388
389 // If this node corresponds to a ghost node, set any "cell" data to be -1.0
390 if (this->IsGhostNode(node_index))
391 {
392 // Populate the vector of VTK cell data
393 vtk_cell_data[node_index] = -1.0;
394 }
395 else
396 {
397 // Get the cell corresponding to this node
398 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
399
400 // Populate the vector of VTK cell data
401 vtk_cell_data[node_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(p_cell, this);
402 }
403 }
404
405 mesh_writer.AddPointData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
406 }
407
408 // Next, record which nodes are ghost nodes
409 // Note that the cell writer hierarchy can not be used to do this as ghost nodes don't have corresponding cells.
410 std::vector<double> ghosts(num_vtk_cells);
411 for (typename AbstractMesh<DIM, DIM>::NodeIterator node_iter = this->rGetMesh().GetNodeIteratorBegin();
412 node_iter != this->rGetMesh().GetNodeIteratorEnd();
413 ++node_iter)
414 {
415 unsigned node_index = node_iter->GetIndex();
416 ghosts[node_index] = (double) (this->IsGhostNode(node_index));
417 }
418 mesh_writer.AddPointData("Non-ghosts", ghosts);
419
421
422 mesh_writer.WriteFilesUsingMesh(this->rGetMesh());
423 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
424 *(this->mpVtkMetaFile) << num_timesteps;
425 *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"mesh_results_";
426 *(this->mpVtkMetaFile) << num_timesteps;
427 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
428 }
429 if (this->mpVoronoiTessellation != nullptr)
430 {
431 // Create mesh writer for VTK output
432 VertexMeshWriter<DIM, DIM> mesh_writer(rDirectory, "voronoi_results", false);
433
434 // Iterate over any cell writers that are present
435 unsigned num_vtk_cells = this->mpVoronoiTessellation->GetNumElements();
436 for (typename std::vector<boost::shared_ptr<AbstractCellWriter<DIM, DIM> > >::iterator cell_writer_iter = this->mCellWriters.begin();
437 cell_writer_iter != this->mCellWriters.end();
438 ++cell_writer_iter)
439 {
440 // Create vector to store VTK cell data
441 std::vector<double> vtk_cell_data(num_vtk_cells);
442
443 // Loop over elements of mpVoronoiTessellation
444 for (typename VertexMesh<DIM, DIM>::VertexElementIterator elem_iter = this->mpVoronoiTessellation->GetElementIteratorBegin();
445 elem_iter != this->mpVoronoiTessellation->GetElementIteratorEnd();
446 ++elem_iter)
447 {
448 // Get the indices of this element and the corresponding node in mrMesh
449 unsigned elem_index = elem_iter->GetIndex();
450 unsigned node_index = this->mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
451
452 // If this node corresponds to a ghost node, set any "cell" data to be -1.0
453 if (this->IsGhostNode(node_index))
454 {
455 // Populate the vector of VTK cell data
456 vtk_cell_data[elem_index] = -1.0;
457 }
458 else
459 {
460 // Get the cell corresponding to this node
461 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
462
463 // Populate the vector of VTK cell data
464 vtk_cell_data[elem_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(p_cell, this);
465 }
466 }
467
468 mesh_writer.AddCellData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
469 }
470
471 // Next, record which nodes are ghost nodes
472 // Note that the cell writer hierarchy can not be used to do this as ghost nodes don't have corresponding cells.
473 std::vector<double> ghosts(num_vtk_cells);
474 for (typename VertexMesh<DIM, DIM>::VertexElementIterator elem_iter = this->mpVoronoiTessellation->GetElementIteratorBegin();
475 elem_iter != this->mpVoronoiTessellation->GetElementIteratorEnd();
476 ++elem_iter)
477 {
478 unsigned elem_index = elem_iter->GetIndex();
479 unsigned node_index = this->mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
480 ghosts[elem_index] = (double) (this->IsGhostNode(node_index));
481 }
482 mesh_writer.AddCellData("Non-ghosts", ghosts);
483
485
486 mesh_writer.WriteVtkUsingMesh(*(this->mpVoronoiTessellation), time.str());
487 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
488 *(this->mpVtkMetaFile) << num_timesteps;
489 *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"voronoi_results_";
490 *(this->mpVtkMetaFile) << num_timesteps;
491 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
492 }
493#endif //CHASTE_VTK
494}
495
496template<unsigned DIM>
498{
499 *rParamsFile << "\t\t<GhostCellSpringStiffness>" << mGhostCellSpringStiffness << "</GhostCellSpringStiffness>\n";
500 *rParamsFile << "\t\t<GhostGhostSpringStiffness>" << mGhostGhostSpringStiffness << "</GhostGhostSpringStiffness>\n";
501 *rParamsFile << "\t\t<GhostSpringRestLength>" << mGhostSpringRestLength << "</GhostSpringRestLength>\n";
502
503 // Call method on direct parent class
505}
506
507// Explicit instantiation
511
512// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
std::set< unsigned > GetNeighbouringLocationIndices(CellPtr pCell)
virtual void OpenWritersFiles(OutputFileHandler &rOutputFileHandler)
MeshBasedCellPopulationWithGhostNodes(MutableMesh< DIM, DIM > &rMesh, std::vector< CellPtr > &rCells, const std::vector< unsigned > locationIndices=std::vector< unsigned >(), bool deleteMesh=false, double ghostCellSpringStiffness=15.0, double ghostGhostSpringStiffness=15.0, double ghostSpringRestLength=1.0)
c_vector< double, DIM > CalculateForceBetweenGhostNodes(const unsigned &rNodeAGlobalIndex, const unsigned &rNodeBGlobalIndex)
virtual void WriteVtkResultsToFile(const std::string &rDirectory)
CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell)
void SetGhostNodes(const std::set< unsigned > &rGhostNodeIndices)
virtual TetrahedralMesh< DIM, DIM > * GetTetrahedralMeshForPdeModifier()
virtual void OpenWritersFiles(OutputFileHandler &rOutputFileHandler)
void OutputCellPopulationParameters(out_stream &rParamsFile)
virtual CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell)
Node< SPACE_DIM > * GetNode(unsigned index)
void DeleteNodePriorToReMesh(unsigned index)
unsigned GetNewIndex(unsigned oldIndex) const
Definition NodeMap.cpp:87
unsigned GetSize()
Definition NodeMap.cpp:105
bool IsDeleted(unsigned index)
Definition NodeMap.cpp:82
static SimulationTime * Instance()
unsigned GetTimeStepsElapsed() const
EdgeIterator EdgesEnd()
void WriteVtkUsingMesh(VertexMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, std::string stamp="")
void AddCellData(std::string dataName, std::vector< double > dataPayload)
void AddPointData(std::string name, std::vector< double > data)
void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)