Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ImmersedBoundaryCellPopulation.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 "ImmersedBoundaryCellPopulation.hpp"
37
38#include <iomanip>
39
40#include "ApoptoticCellProperty.hpp"
41#include "CellPopulationElementWriter.hpp"
42#include "ImmersedBoundaryMeshWriter.hpp"
43#include "ImmersedBoundaryBoundaryCellWriter.hpp"
44#include "ShortAxisImmersedBoundaryDivisionRule.hpp"
45#include "StepSizeException.hpp"
46#include "Warnings.hpp"
47
48template <unsigned DIM>
51 std::vector<CellPtr>& rCells,
52 bool deleteMesh,
53 bool validate,
54 const std::vector<unsigned> locationIndices)
55 : AbstractOffLatticeCellPopulation<DIM>(rMesh, rCells, locationIndices),
56 mDeleteMesh(deleteMesh),
57 mIntrinsicSpacing(0.01),
58 mPopulationHasActiveSources(false),
59 mOutputNodeRegionToVtk(false),
60 mReMeshFrequency(UINT_MAX),
61 mThrowStepSizeException(true),
62 mCellRearrangementThreshold(0.5)
63{
66
67 /*
68 * If no location indices are specified, associate with elements from the
69 * mesh (assumed to be sequentially ordered).
70 */
71 std::list<CellPtr>::iterator it = this->mCells.begin();
72 for (unsigned i = 0; it != this->mCells.end(); ++it, ++i)
73 {
74 // Assume that the ordering matches
75 unsigned index = locationIndices.empty() ? i : locationIndices[i];
77 }
78
79 // Check each element has only one cell attached
80 if (validate)
81 {
82 Validate();
83 }
84
86
87 // Set the mesh division spacing distance
90}
91
92template <unsigned DIM>
96 mDeleteMesh(true),
97 mIntrinsicSpacing(0.01),
98 mPopulationHasActiveSources(false),
99 mOutputNodeRegionToVtk(false),
100 mReMeshFrequency(UINT_MAX),
101 mThrowStepSizeException(true),
102 mCellRearrangementThreshold(0.5)
103{
105}
106
107template <unsigned DIM>
109{
110 if (mDeleteMesh)
111 {
112 delete &this->mrMesh;
113 }
114}
115
116template <unsigned DIM>
118{
119 double average_intrinsic_size = 0.0;
120
121 for (auto elem_iter = mpImmersedBoundaryMesh->GetElementIteratorBegin();
122 elem_iter != mpImmersedBoundaryMesh->GetElementIteratorEnd();
123 ++elem_iter)
124 {
125 average_intrinsic_size += mpImmersedBoundaryMesh->GetVolumeOfElement(elem_iter->GetIndex());
126 }
127
128 return sqrt(average_intrinsic_size / mpImmersedBoundaryMesh->GetNumElements());
129}
130
131template <unsigned DIM>
133{
134 return 0.0;
135}
136
137template <unsigned DIM>
142
143template <unsigned DIM>
145{
146 return *mpImmersedBoundaryMesh;
147}
148
149template <unsigned DIM>
151{
152 return mpImmersedBoundaryMesh->GetElement(elementIndex);
153}
154
155template <unsigned DIM>
157{
158 return mpImmersedBoundaryMesh->GetLamina(laminaIndex);
159}
160
161template <unsigned DIM>
163{
164 return this->mrMesh.GetNumNodes();
165}
166
167template <unsigned DIM>
169{
170 return mpImmersedBoundaryMesh->GetCentroidOfElement(this->mCellLocationMap[pCell.get()]);
171}
172
173template <unsigned DIM>
175{
176 return this->mrMesh.GetNode(index);
177}
178
179template <unsigned DIM>
181{
182 assert(newDistance >= 0.0);
183 mInteractionDistance = newDistance;
184}
185
186template <unsigned DIM>
188{
189 return mInteractionDistance;
190}
191
192template <unsigned DIM>
194{
195 mReMeshFrequency = newFrequency;
196}
197
198template <unsigned DIM>
200{
201 return mReMeshFrequency;
202}
203
204template <unsigned DIM>
206{
207 mCellRearrangementThreshold = newThreshold;
208}
209
210template <unsigned DIM>
212{
213 return mCellRearrangementThreshold;
214}
215
216template <unsigned DIM>
218{
219 mThrowStepSizeException = throws;
220}
221
222template <unsigned DIM>
224{
225 return mThrowStepSizeException;
226}
227
228template <unsigned DIM>
230{
231 return mIntrinsicSpacing;
232}
233
234template <unsigned DIM>
236 CellPtr pCell)
237{
238 return mpImmersedBoundaryMesh->GetNeighbouringElementIndices(this->GetLocationIndexUsingCell(pCell));
239}
240
241template <unsigned DIM>
243{
244 if (pNewNode == nullptr)
245 {
246 return 0;
247 }
248 return mpImmersedBoundaryMesh->AddNode(pNewNode);
249}
250
251template <unsigned DIM>
253{
254 mpImmersedBoundaryMesh->SetNode(nodeIndex, rNewLocation);
255}
256
257template <unsigned DIM>
259 CellPtr pCell)
260{
261 return mpImmersedBoundaryMesh->GetElement(this->GetLocationIndexUsingCell(pCell));
262}
263
264template <unsigned DIM>
266{
267 return mpImmersedBoundaryMesh->GetNumElements();
268}
269
270template <unsigned DIM>
272{
273 return mpImmersedBoundaryMesh->GetNumLaminas();
274}
275
276template <unsigned DIM>
278 CellPtr pNewCell,
279 CellPtr pParentCell)
280{
281 // Get the element associated with this cell
282 ImmersedBoundaryElement<DIM, DIM>* p_element = GetElementCorrespondingToCell(pParentCell);
283
284 // Get the orientation of division
285 c_vector<double, DIM> division_vector = mpImmersedBoundaryDivisionRule->CalculateCellDivisionVector(pParentCell, *this);
286
287 // Divide the element
288 unsigned new_elem_idx = mpImmersedBoundaryMesh->DivideElementAlongGivenAxis(p_element, division_vector, true);
289
290 // Associate the new cell with the element
291 this->mCells.push_back(pNewCell);
292
293 // Update location cell map
294 CellPtr p_created_cell = this->mCells.back();
295 this->SetCellUsingLocationIndex(new_elem_idx, p_created_cell);
296 this->mCellLocationMap[p_created_cell.get()] = new_elem_idx;
297
298 return p_created_cell;
299}
300
301template <unsigned DIM>
303{
304 unsigned num_removed = 0;
305
306 for (std::list<CellPtr>::iterator it = this->mCells.begin();
307 it != this->mCells.end();
308 )
309 {
310 if ((*it)->IsDead())
311 {
312 // Count the cell as dead
313 num_removed++;
314
315 if (!(this->GetElement(this->GetLocationIndexUsingCell((*it)))->IsDeleted()))
316 {
317 this->GetElement(this->GetLocationIndexUsingCell(*it))->MarkAsDeleted();
318 }
319
320 // Delete the cell
321 it = this->mCells.erase(it);
322 }
323 else
324 {
325 ++it;
326 }
327 }
328 return num_removed;
329}
330
331template <unsigned DIM>
333 [[maybe_unused]] double dt) // [[maybe_unused]] due to unused-but-set-parameter warning in GCC 7,8,9
334{
335 if constexpr (DIM == 2)
336 {
337 // Helper variables, pre-declared for efficiency
338 unsigned num_grid_pts_x = this->rGetMesh().GetNumGridPtsX();
339 unsigned num_grid_pts_y = this->rGetMesh().GetNumGridPtsY();
340
341 double characteristic_spacing = this->rGetMesh().GetCharacteristicNodeSpacing();
342 double grid_spacing_x = 1.0 / (double)num_grid_pts_x;
343 double grid_spacing_y = 1.0 / (double)num_grid_pts_y;
344
345 unsigned first_idx_x;
346 unsigned first_idx_y;
347
348 std::vector<unsigned> x_indices(4);
349 std::vector<unsigned> y_indices(4);
350
351 std::vector<double> x_deltas(4);
352 std::vector<double> y_deltas(4);
353
354 double delta;
355
356 c_vector<double, DIM> displacement = zero_vector<double>(DIM);
357
358 // Get references to the fluid velocity grid
359 const multi_array<double, 3>& vel_grids = this->rGetMesh().rGet2dVelocityGrids();
360
361 // Iterate over all nodes
362 c_vector<double, DIM> node_location;
363 for (auto node_iter = this->rGetMesh().GetNodeIteratorBegin(false);
364 node_iter != this->rGetMesh().GetNodeIteratorEnd();
365 ++node_iter)
366 {
367 // Get location of current node
368 node_location = node_iter->rGetLocation();
369
370 // Get first grid index in each dimension, taking account of possible wrap-around
371 first_idx_x = unsigned(floor(node_location[0] / grid_spacing_x)) + num_grid_pts_x - 1;
372 first_idx_y = unsigned(floor(node_location[1] / grid_spacing_y)) + num_grid_pts_y - 1;
373
374 // Calculate all four indices and deltas in each dimension
375 for (unsigned i = 0; i < 4; i++)
376 {
377 x_indices[i] = (first_idx_x + i) % num_grid_pts_x;
378 y_indices[i] = (first_idx_y + i) % num_grid_pts_y;
379
380 x_deltas[i] = Delta1D(fabs(x_indices[i] * grid_spacing_x - node_location[0]), grid_spacing_x);
381 y_deltas[i] = Delta1D(fabs(y_indices[i] * grid_spacing_x - node_location[1]), grid_spacing_y);
382 }
383
384 // Loop over the 4x4 grid which will influence the displacement of the current node
385 for (unsigned x_idx = 0; x_idx < 4; ++x_idx)
386 {
387 for (unsigned y_idx = 0; y_idx < 4; ++y_idx)
388 {
389 // The applied velocity is weighted by the delta function
390 delta = x_deltas[x_idx] * y_deltas[y_idx];
391 displacement[0] += vel_grids[0][x_indices[x_idx]][y_indices[y_idx]] * delta;
392 displacement[1] += vel_grids[1][x_indices[x_idx]][y_indices[y_idx]] * delta;
393 }
394 }
395
396 // Normalise by timestep
397 displacement *= dt;
398
399 // If the displacement is too big, warn the user once and scale it back
400 if (norm_2(displacement) > characteristic_spacing)
401 {
402 if (norm_2(displacement) > 10.0 * characteristic_spacing)
403 {
404 EXCEPTION("Nodes are moving more than 10x CharacteristicNodeSpacing. Aborting.");
405 }
406
407 WARN_ONCE_ONLY("Nodes are moving more than the CharacteristicNodeSpacing. This could cause elements to become inverted so the motion has been restricted. Use a smaller timestep to avoid these warnings.");
408 displacement *= characteristic_spacing / norm_2(displacement);
409 }
410
411 CheckForStepSizeException(node_iter->GetIndex(), displacement, dt);
412
413 // Get new node location
414 node_location += displacement;
415
416 // Account for periodic boundary
417 for (unsigned i = 0; i < DIM; ++i)
418 {
419 node_location[i] = fmod(node_location[i] + 1.0, 1.0);
420 }
421
422 // Create ChastePoint for new node location
423 ChastePoint<DIM> new_point(node_location);
424
425 // Move the node
426 this->SetNode(node_iter->GetIndex(), new_point);
427 }
428
429 // If active sources, we need to update those location as well
430 if (this->DoesPopulationHaveActiveSources())
431 {
432 std::vector<std::shared_ptr<FluidSource<DIM>>>& r_element_sources = this->rGetMesh().rGetElementFluidSources();
433 std::vector<std::shared_ptr<FluidSource<DIM>>>& r_balance_sources = this->rGetMesh().rGetBalancingFluidSources();
434
435 // Construct a vector of all sources combined
436 std::vector<std::shared_ptr<FluidSource<DIM>>> combined_sources;
437 combined_sources.insert(combined_sources.end(), r_element_sources.begin(), r_element_sources.end());
438 combined_sources.insert(combined_sources.end(), r_balance_sources.begin(), r_balance_sources.end());
439
440 c_vector<double, DIM> source_location;
441
442 // Iterate over all sources and update their locations
443 for (unsigned source_idx = 0; source_idx < combined_sources.size(); source_idx++)
444 {
445 // Get location of current node
446 source_location = combined_sources[source_idx]->rGetLocation();
447
448 // Get first grid index in each dimension, taking account of possible wrap-around
449 first_idx_x = unsigned(floor(source_location[0] / grid_spacing_x)) + num_grid_pts_x - 1;
450 first_idx_y = unsigned(floor(source_location[1] / grid_spacing_y)) + num_grid_pts_y - 1;
451
452 // Calculate all four indices and deltas in each dimension
453 for (unsigned i = 0; i < 4; ++i)
454 {
455 x_indices[i] = (first_idx_x + i) % num_grid_pts_x;
456 y_indices[i] = (first_idx_y + i) % num_grid_pts_y;
457
458 x_deltas[i] = Delta1D(fabs(x_indices[i] * grid_spacing_x - source_location[0]), grid_spacing_x);
459 y_deltas[i] = Delta1D(fabs(y_indices[i] * grid_spacing_x - source_location[1]), grid_spacing_y);
460 }
461
462 // Loop over the 4x4 grid which will influence the displacement of the current node
463 for (unsigned x_idx = 0; x_idx < 4; ++x_idx)
464 {
465 for (unsigned y_idx = 0; y_idx < 4; ++y_idx)
466 {
467 // The applied velocity is weighted by the delta function
468 delta = x_deltas[x_idx] * y_deltas[y_idx];
469 displacement[0] += vel_grids[0][x_indices[x_idx]][y_indices[y_idx]] * delta;
470 displacement[1] += vel_grids[1][x_indices[x_idx]][y_indices[y_idx]] * delta;
471 }
472 }
473
474 // Normalise by timestep
475 displacement *= dt;
476
477 // If the displacement is too big, warn the user once and scale it back
478 if (norm_2(displacement) > characteristic_spacing)
479 {
480 if (norm_2(displacement) > 10.0 * characteristic_spacing)
481 {
482 EXCEPTION("Sources are moving more than 10x CharacteristicNodeSpacing. Aborting.");
483 }
484
485 WARN_ONCE_ONLY("Sources are moving more than the CharacteristicNodeSpacing. This could cause elements to become inverted so the motion has been restricted. Use a smaller timestep to avoid these warnings.");
486 displacement *= characteristic_spacing / norm_2(displacement);
487 }
488
489 // Get new node location
490 source_location += displacement;
491
492 // Account for periodic boundary
493 for (unsigned i = 0; i < DIM; ++i)
494 {
495 source_location[i] = fmod(source_location[i] + 1.0, 1.0);
496 }
497
498 // Move the node
499 combined_sources[source_idx]->rGetModifiableLocation() = source_location;
500 }
501 }
502
503 // Finally, call ReMesh if required
504 const auto num_time_steps = SimulationTime::Instance()->GetTimeStepsElapsed();
505 if (num_time_steps > 0 && num_time_steps % mReMeshFrequency == 0)
506 {
507 mpImmersedBoundaryMesh->ReMesh();
508 }
509 }
510 else
511 {
513 }
514}
515
516template <unsigned DIM>
517double ImmersedBoundaryCellPopulation<DIM>::Delta1D(double dist, double spacing)
518{
519 return (0.25 * (1.0 + cos(M_PI * dist / (2 * spacing))));
520}
521
522template <unsigned DIM>
524{
525 return GetElementCorrespondingToCell(pCell)->IsDeleted();
526}
527
528template <unsigned DIM>
529void ImmersedBoundaryCellPopulation<DIM>::Update(bool hasHadBirthsOrDeaths)
530{
531 // If the first cell has target atea property, assume there is a target area modifier in place
532 if (this->Begin()->GetCellData()->HasItem("target area"))
533 {
534 for (auto cell_iter = this->Begin();
535 cell_iter != this->End();
536 ++cell_iter)
537 {
538 double target_area = cell_iter->GetCellData()->GetItem("target area");
539 double actual_area = this->GetVolumeOfCell(*cell_iter);
540
541 double strength = 1e-2 * (target_area - actual_area) / target_area;
542
543 this->GetElementCorrespondingToCell(*cell_iter)->GetFluidSource()->SetStrength(strength);
544 }
545 }
546}
547
548template <unsigned DIM>
550{
551 // Check each element has only one cell attached
552 std::vector<unsigned> validated_element = std::vector<unsigned>(this->GetNumElements(), 0);
553 for (auto cell_iter = this->Begin();
554 cell_iter != this->End();
555 ++cell_iter)
556 {
557 unsigned elem_index = this->GetLocationIndexUsingCell(*cell_iter);
558 validated_element[elem_index]++;
559 }
560
561 for (unsigned i = 0; i < validated_element.size(); ++i)
562 {
563 if (validated_element[i] == 0)
564 {
565 EXCEPTION("At time " << SimulationTime::Instance()->GetTime() << ", Element " << i << " does not appear to have a cell associated with it");
566 }
567
568 if (validated_element[i] > 1)
569 {
570 // This should never be reached as you can only set one cell per element index
572 EXCEPTION("At time " << SimulationTime::Instance()->GetTime() << ", Element " << i << " appears to have " << validated_element[i] << " cells associated with it"); //LCOV_EXCL_LINE
573 }
574 }
575}
576
577template <unsigned DIM>
579 unsigned nodeIndex,
580 c_vector<double, DIM>& rDisplacement,
581 double dt)
582{
583 double length = boost::numeric::ublas::norm_2(rDisplacement);
584
585 if (length > 0.5*mpImmersedBoundaryMesh->GetCellRearrangementThreshold())
586 {
587 rDisplacement *= 0.5*mpImmersedBoundaryMesh->GetCellRearrangementThreshold()/length;
588
589 std::ostringstream message;
590 message << "Vertices are moving more than half the CellRearrangementThreshold. This could cause elements to become inverted ";
591 message << "so the motion has been restricted. Use a smaller timestep to avoid these warnings.";
592
593 double suggested_step = 0.95*dt*((0.5*mpImmersedBoundaryMesh->GetCellRearrangementThreshold())/length);
594
595 // The first time we see this behaviour, throw a StepSizeException, but not more than once
596 if(mThrowStepSizeException)
597 {
598 mThrowStepSizeException = false;
599 throw StepSizeException(suggested_step, message.str(), false);
600 }
601 }
602}
603
604template <unsigned DIM>
606 boost::shared_ptr<AbstractCellPopulationWriter<DIM, DIM> > pPopulationWriter)
607{
608 pPopulationWriter->Visit(this);
609}
610
611template <unsigned DIM>
613 boost::shared_ptr<AbstractCellPopulationEventWriter<DIM, DIM> > pPopulationEventWriter)
614{
615 pPopulationEventWriter->Visit(this);
616}
617
618template <unsigned DIM>
620 boost::shared_ptr<AbstractCellPopulationCountWriter<DIM, DIM> > pPopulationCountWriter)
621{
622 pPopulationCountWriter->Visit(this);
623}
624
625template <unsigned DIM>
627 boost::shared_ptr<AbstractCellWriter<DIM, DIM> > pCellWriter, CellPtr pCell)
628{
629 pCellWriter->VisitCell(pCell, this);
630}
631
632template <unsigned DIM>
634{
635 // Get the element index corresponding to this cell
636 unsigned elem_index = this->GetLocationIndexUsingCell(pCell);
637
638 // Get the cell's volume from the immersed boundary mesh
639 double cell_volume = mpImmersedBoundaryMesh->GetVolumeOfElement(elem_index);
640
641 return cell_volume;
642}
643
644template <unsigned DIM>
646 const std::string& rDirectory)
647{
648#ifdef CHASTE_VTK
649 // Create mesh writer for VTK output
650 ImmersedBoundaryMeshWriter<DIM, DIM> mesh_writer(rDirectory, "results", false);
651
652 // Find the cell overlap information, and get the number of cell parts needed for each element
653 mesh_writer.FindElementOverlaps(*mpImmersedBoundaryMesh);
654 const std::vector<std::vector<unsigned>>& r_elem_parts = mesh_writer.rGetElementParts();
655
656 // Iterate over any cell writers that are present
657 for (auto cell_writer_iter = this->mCellWriters.begin();
658 cell_writer_iter != this->mCellWriters.end();
659 ++cell_writer_iter)
660 {
661 // Create vector to store VTK cell data
662 std::vector<double> vtk_cell_data;
663
664 // Iterate over immersed boundary elements
665 for (auto elem_iter = mpImmersedBoundaryMesh->GetElementIteratorBegin();
666 elem_iter != mpImmersedBoundaryMesh->GetElementIteratorEnd();
667 ++elem_iter)
668 {
669 /*
670 * Get index of this element in the mesh, and the number of parts it
671 * is broken into for visualisation.
672 */
673 const unsigned elem_index = elem_iter->GetIndex();
674 const auto num_elem_parts = r_elem_parts[elem_index].empty() ? 1 : r_elem_parts[elem_index].size();
675
676 // Get the cell corresponding to this element
677 CellPtr p_cell = this->GetCellUsingLocationIndex(elem_index);
678 assert(p_cell);
679
680 /*
681 * Populate the vector of VTK cell data. We loop over the number of
682 * output cells as this takes into account that some elements will
683 * be broken into pieces for visualisation.
684 */
685 for (unsigned elem_part = 0; elem_part < num_elem_parts; ++elem_part)
686 {
687 vtk_cell_data.push_back((*cell_writer_iter)->GetCellDataForVtkOutput(p_cell, this));
688 }
689 }
690
691 /*
692 * Iterate over immersed boundary laminas (no associated cell) to ensure
693 * vtk_cell_data is the correct size.
694 */
695 for (auto lam_iter = mpImmersedBoundaryMesh->GetLaminaIteratorBegin();
696 lam_iter != mpImmersedBoundaryMesh->GetLaminaIteratorEnd();
697 ++lam_iter)
698 {
699 vtk_cell_data.push_back(-1.0);
700 }
701
702 mesh_writer.AddCellData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
703 }
704
705 /*
706 * When outputting any CellData, we assume that the first cell is
707 * representative of all cells.
708 */
709 const unsigned num_cell_data_items = this->Begin()->GetCellData()->GetNumItems();
710 std::vector<std::string> cell_data_names = this->Begin()->GetCellData()->GetKeys();
711
712 std::vector<std::vector<double>> cell_data;
713 for (unsigned var = 0; var < num_cell_data_items; ++var)
714 {
715 std::vector<double> cell_data_var;
716 cell_data.push_back(cell_data_var);
717 }
718
719 // Iterate over immersed boundary elements
720 for (auto elem_iter = mpImmersedBoundaryMesh->GetElementIteratorBegin();
721 elem_iter != mpImmersedBoundaryMesh->GetElementIteratorEnd();
722 ++elem_iter)
723 {
724 /*
725 * Get index of this element in the mesh, and the number of parts it is
726 * broken into for visualisation.
727 */
728 const unsigned elem_index = elem_iter->GetIndex();
729 const auto num_elem_parts = r_elem_parts[elem_index].empty() ? 1 : r_elem_parts[elem_index].size();
730
731 // Get the cell corresponding to this element
732 CellPtr p_cell = this->GetCellUsingLocationIndex(elem_index);
733 assert(p_cell);
734
735 for (unsigned var = 0; var < num_cell_data_items; var++)
736 {
737 /*
738 * Populate the vector of VTK cell data. We loop over the number of
739 * output cells as this takes into account that some elements will
740 * be broken into pieces for visualisation.
741 */
742 for (unsigned elem_part = 0; elem_part < num_elem_parts; ++elem_part)
743 {
744 cell_data[var].push_back(p_cell->GetCellData()->GetItem(cell_data_names[var]));
745 }
746 }
747 }
748
749 /*
750 * Iterate over immersed boundary laminas (no associated cell) to ensure
751 * cell_data is the correct size.
752 */
753 for (auto lam_iter = mpImmersedBoundaryMesh->GetLaminaIteratorBegin();
754 lam_iter != mpImmersedBoundaryMesh->GetLaminaIteratorEnd();
755 ++lam_iter)
756 {
757 for (unsigned var = 0; var < num_cell_data_items; ++var)
758 {
759 cell_data[var].push_back(DOUBLE_UNSET);
760 }
761 }
762
763 for (unsigned var = 0; var < num_cell_data_items; ++var)
764 {
765 mesh_writer.AddCellData(cell_data_names[var], cell_data[var]);
766 }
767
768 // Write node regions
769 if (mOutputNodeRegionToVtk)
770 {
771 std::vector<double> node_regions;
772 for (auto node_iter = mpImmersedBoundaryMesh->GetNodeIteratorBegin();
773 node_iter != mpImmersedBoundaryMesh->GetNodeIteratorEnd();
774 ++node_iter)
775 {
776 node_regions.push_back(static_cast<double>(node_iter->GetRegion()));
777 }
778 mesh_writer.AddPointData("Node Regions", node_regions);
779 }
780
781 unsigned num_timesteps = SimulationTime::Instance()->GetTimeStepsElapsed();
782 std::stringstream time;
783 time << num_timesteps;
784
785 mesh_writer.WriteVtkUsingMesh(*mpImmersedBoundaryMesh, time.str());
786
787 *(this->mpVtkMetaFile) << " <DataSet timestep=\"";
788 *(this->mpVtkMetaFile) << num_timesteps;
789 *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
790 *(this->mpVtkMetaFile) << num_timesteps;
791 *(this->mpVtkMetaFile) << ".vtu\"/>\n";
792#endif //CHASTE_VTK
793}
794
795template <unsigned DIM>
797 OutputFileHandler& rOutputFileHandler)
798{
799 if (this->mOutputResultsForChasteVisualizer)
800 {
801 if (!this->template HasWriter<CellPopulationElementWriter>())
802 {
803 this->template AddPopulationWriter<CellPopulationElementWriter>();
804 }
805 }
806
808}
809
810template <unsigned DIM>
812 out_stream& rParamsFile)
813{
814 // Add the division rule parameters
815 *rParamsFile << "\t\t<ImmersedBoundaryDivisionRule>\n";
816 mpImmersedBoundaryDivisionRule->OutputCellImmersedBoundaryDivisionRuleInfo(rParamsFile);
817 *rParamsFile << "\t\t</ImmersedBoundaryDivisionRule>\n";
818
819 // Call method on direct parent class
821}
822
823template <unsigned DIM>
824double ImmersedBoundaryCellPopulation<DIM>::GetWidth(const unsigned& rDimension)
825{
826 double width = this->mrMesh.GetWidth(rDimension);
827 return width;
828}
829
830template <unsigned DIM>
832 unsigned index)
833{
834 return mpImmersedBoundaryMesh->GetNeighbouringNodeIndices(index);
835}
836
837template <unsigned DIM>
839{
840 // This method only works in 2D sequential
841 assert(PetscTools::IsSequential());
842 if constexpr (DIM == 2)
843 {
844 unsigned num_vertex_nodes = mpImmersedBoundaryMesh->GetNumNodes();
845 unsigned num_vertex_elements = mpImmersedBoundaryMesh->GetNumElements();
846
847 std::string mesh_file_name = "mesh";
848
849 // Get a unique temporary foldername
850 std::stringstream pid;
851 pid << getpid();
852 OutputFileHandler output_file_handler("2D_temporary_tetrahedral_mesh_" + pid.str());
853 std::string output_dir = output_file_handler.GetOutputDirectoryFullPath();
854
855 // Compute the number of nodes in the TetrahedralMesh
856 unsigned num_tetrahedral_nodes = num_vertex_nodes + num_vertex_elements;
857
858 // Write node file
859 out_stream p_node_file = output_file_handler.OpenOutputFile(mesh_file_name+".node");
860 (*p_node_file) << std::scientific;
861 (*p_node_file) << std::setprecision(20);
862 (*p_node_file) << num_tetrahedral_nodes << "\t2\t0\t1" << std::endl;
863
864 // Begin by writing each node in the VertexMesh
865 auto nodes = mpImmersedBoundaryMesh->rGetNodes();
866 for (auto p_node : nodes)
867 {
868 unsigned index = p_node->GetIndex();
869 const c_vector<double, DIM>& r_location = p_node->rGetLocation();
870 unsigned is_boundary_node = p_node->IsBoundaryNode() ? 1 : 0;
871
872 (*p_node_file) << index << "\t" << r_location[0] << "\t" << r_location[1] << "\t" << is_boundary_node << std::endl;
873 }
874
875 // Now write an additional node at each ImmersedBoundaryElement's centroid
876 unsigned num_tetrahedral_elements = 0;
877 for (unsigned vertex_elem_index = 0;
878 vertex_elem_index < num_vertex_elements;
879 ++vertex_elem_index)
880 {
881 unsigned index = num_vertex_nodes + vertex_elem_index;
882
883 c_vector<double, DIM> location = mpImmersedBoundaryMesh->GetCentroidOfElement(vertex_elem_index);
884
885 // Any node located at a ImmersedBoundaryElement's centroid will not be a boundary node
886 unsigned is_boundary_node = 0;
887 (*p_node_file) << index << "\t" << location[0] << "\t" << location[1] << "\t" << is_boundary_node << std::endl;
888
889 // Also keep track of how many tetrahedral elements there will be
890 num_tetrahedral_elements += mpImmersedBoundaryMesh->GetElement(vertex_elem_index)->GetNumNodes();
891 }
892 p_node_file->close();
893
894 // Write element file
895 out_stream p_elem_file = output_file_handler.OpenOutputFile(mesh_file_name+".ele");
896 (*p_elem_file) << std::scientific;
897 (*p_elem_file) << num_tetrahedral_elements << "\t3\t0" << std::endl;
898
899 std::set<std::pair<unsigned, unsigned> > tetrahedral_edges;
900
901 unsigned tetrahedral_elem_index = 0;
902 for (unsigned vertex_elem_index = 0;
903 vertex_elem_index < num_vertex_elements;
904 ++vertex_elem_index)
905 {
906 ImmersedBoundaryElement<DIM, DIM>* p_vertex_element = mpImmersedBoundaryMesh->GetElement(vertex_elem_index);
907
908 // Iterate over nodes owned by this ImmersedBoundaryElement
909 unsigned num_nodes_in_vertex_element = p_vertex_element->GetNumNodes();
910 for (unsigned local_index = 0;
911 local_index < num_nodes_in_vertex_element;
912 ++local_index)
913 {
914 unsigned node_0_index = p_vertex_element->GetNodeGlobalIndex(local_index);
915 unsigned node_1_index = p_vertex_element->GetNodeGlobalIndex((local_index+1)%num_nodes_in_vertex_element);
916 unsigned node_2_index = num_vertex_nodes + vertex_elem_index;
917
918 (*p_elem_file) << tetrahedral_elem_index++ << "\t" << node_0_index << "\t" << node_1_index << "\t" << node_2_index << std::endl;
919
920 // Add edges to the set if they are not already present
921 std::pair<unsigned, unsigned> edge_0 = this->CreateOrderedPair(node_0_index, node_1_index);
922 std::pair<unsigned, unsigned> edge_1 = this->CreateOrderedPair(node_1_index, node_2_index);
923 std::pair<unsigned, unsigned> edge_2 = this->CreateOrderedPair(node_2_index, node_0_index);
924
925 tetrahedral_edges.insert(edge_0);
926 tetrahedral_edges.insert(edge_1);
927 tetrahedral_edges.insert(edge_2);
928 }
929 }
930 p_elem_file->close();
931
932 // Write edge file
933 out_stream p_edge_file = output_file_handler.OpenOutputFile(mesh_file_name+".edge");
934 (*p_edge_file) << std::scientific;
935 (*p_edge_file) << tetrahedral_edges.size() << "\t1" << std::endl;
936
937 unsigned edge_index = 0;
938 for (auto edge_iter = tetrahedral_edges.begin();
939 edge_iter != tetrahedral_edges.end();
940 ++edge_iter)
941 {
942 std::pair<unsigned, unsigned> this_edge = *edge_iter;
943
944 // To be a boundary edge both nodes need to be boundary nodes.
945 bool is_boundary_edge = false;
946 if (this_edge.first < mpImmersedBoundaryMesh->GetNumNodes() &&
947 this_edge.second < mpImmersedBoundaryMesh->GetNumNodes())
948 {
949 is_boundary_edge = (mpImmersedBoundaryMesh->GetNode(this_edge.first)->IsBoundaryNode() &&
950 mpImmersedBoundaryMesh->GetNode(this_edge.second)->IsBoundaryNode() );
951 }
952 unsigned is_boundary_edge_unsigned = is_boundary_edge ? 1 : 0;
953
954 (*p_edge_file) << edge_index++ << "\t" << this_edge.first << "\t" << this_edge.second << "\t" << is_boundary_edge_unsigned << std::endl;
955 }
956 p_edge_file->close();
957
958 // Having written the mesh to file, now construct it using TrianglesMeshReader
960
961 // Nested scope so reader is destroyed before we remove the temporary files
962 {
963 TrianglesMeshReader<DIM, DIM> mesh_reader(output_dir + mesh_file_name);
964 p_mesh->ConstructFromMeshReader(mesh_reader);
965 }
966
967 // Delete the temporary files
968 output_file_handler.FindFile("").Remove();
969
970 /*
971 * The original files have been deleted, it is better if the mesh object
972 * forgets about them.
973 */
975
976 return p_mesh;
977 }
978 else
979 {
980 EXCEPTION("ImmersedBoundaryCellPopulation::GetTetrahedralMeshForPDEModifier is only implemented in 2D");
981 }
982}
983
984template <unsigned DIM>
986{
987 bool non_apoptotic_cell_present = true;
988
989 if (pdeNodeIndex < this->GetNumNodes())
990 {
991 std::set<unsigned> containing_element_indices = this->GetNode(pdeNodeIndex)->rGetContainingElementIndices();
992
993 for (auto iter = containing_element_indices.begin();
994 iter != containing_element_indices.end();
995 iter++) //LCOV_EXCL_LINE
996 {
997 if (this->GetCellUsingLocationIndex(*iter)->template HasCellProperty<ApoptoticCellProperty>() )
998 {
999 non_apoptotic_cell_present = false;
1000 break;
1001 }
1002 }
1003 }
1004 else
1005 {
1006 /*
1007 * This node of the tetrahedral finite element mesh is in the centre of
1008 * the element of the immersed boundary-based cell population, so we can use an
1009 * offset to compute which cell to interrogate.
1010 */
1011 non_apoptotic_cell_present = !(this->GetCellUsingLocationIndex(pdeNodeIndex - this->GetNumNodes())->template HasCellProperty<ApoptoticCellProperty>());
1012 }
1013
1014 return non_apoptotic_cell_present;
1015}
1016
1017template <unsigned DIM>
1019 unsigned pdeNodeIndex,
1020 std::string& rVariableName,
1021 bool dirichletBoundaryConditionApplies,
1022 double dirichletBoundaryValue)
1023{
1024 unsigned num_nodes = this->GetNumNodes();
1025 double value = 0.0;
1026
1027 /*
1028 * Cells correspond to nodes in the centre of the vertex element; nodes on
1029 * vertices have averaged values from containing cells.
1030 */
1031 if (pdeNodeIndex >= num_nodes)
1032 {
1033 // Offset to relate elements in vertex mesh to nodes in tetrahedral mesh
1034 assert(pdeNodeIndex-num_nodes < num_nodes);
1035
1036 CellPtr p_cell = this->GetCellUsingLocationIndex(pdeNodeIndex - num_nodes);
1037 value = p_cell->GetCellData()->GetItem(rVariableName);
1038 }
1039 else
1040 {
1042 if (dirichletBoundaryConditionApplies)
1043 {
1044 // We need to impose the Dirichlet boundaries again here as not represented in cell data
1045 value = dirichletBoundaryValue;
1046 }
1047 else
1048 {
1049 assert(pdeNodeIndex < num_nodes);
1050 Node<DIM>* p_node = this->GetNode(pdeNodeIndex);
1051
1052 // Average over data from containing elements (cells)
1053 std::set<unsigned> containing_elements = p_node->rGetContainingElementIndices();
1054 for (auto index_iter = containing_elements.begin();
1055 index_iter != containing_elements.end();
1056 ++index_iter)
1057 {
1058 assert(*index_iter < num_nodes);
1059 CellPtr p_cell = this->GetCellUsingLocationIndex(*index_iter);
1060 value += p_cell->GetCellData()->GetItem(rVariableName);
1061 }
1062 value /= containing_elements.size();
1063 }
1064 }
1065
1066 return value;
1067}
1068
1069template <unsigned DIM>
1070boost::shared_ptr<AbstractImmersedBoundaryDivisionRule<DIM> > ImmersedBoundaryCellPopulation<DIM>::GetImmersedBoundaryDivisionRule()
1071{
1072 return mpImmersedBoundaryDivisionRule;
1073}
1074
1075template <unsigned DIM>
1077 boost::shared_ptr<AbstractImmersedBoundaryDivisionRule<DIM> > pImmersedBoundaryDivisionRule)
1078{
1079 mpImmersedBoundaryDivisionRule = pImmersedBoundaryDivisionRule;
1080}
1081
1082template <unsigned DIM>
1084{
1085 return mPopulationHasActiveSources;
1086}
1087
1088template <unsigned DIM>
1090{
1091 return this->GetElementCorrespondingToCell(pCell)->IsElementOnBoundary();
1092}
1093
1094template <unsigned DIM>
1096 bool hasActiveSources)
1097{
1098 mPopulationHasActiveSources = hasActiveSources;
1099}
1100
1101template <unsigned DIM>
1103 bool outputNodeRegionsToVtk)
1104{
1105 mOutputNodeRegionToVtk = outputNodeRegionsToVtk;
1106}
1107
1108template <unsigned DIM>
1110{
1111 return 0.002;
1112}
1113
1114// Explicit instantiation
1118
1119// Serialization for Boost >= 1.36
const double DOUBLE_UNSET
Definition Exception.hpp:57
#define EXCEPTION(message)
#define NEVER_REACHED
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
virtual void AddCellUsingLocationIndex(unsigned index, CellPtr pCell)
virtual void OpenWritersFiles(OutputFileHandler &rOutputFileHandler)
AbstractMesh< ELEMENT_DIM, SPACE_DIM > & mrMesh
unsigned GetNumNodes() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
void SetMeshHasChangedSinceLoading()
virtual void OutputCellPopulationParameters(out_stream &rParamsFile)
void Remove() const
virtual void OpenWritersFiles(OutputFileHandler &rOutputFileHandler)
ImmersedBoundaryCellPopulation(ImmersedBoundaryMesh< DIM, DIM > &rMesh, std::vector< CellPtr > &rCells, bool deleteMesh=false, bool validate=true, const std::vector< unsigned > locationIndices=std::vector< unsigned >())
double Delta1D(double dist, double spacing)
boost::shared_ptr< AbstractImmersedBoundaryDivisionRule< DIM > > mpImmersedBoundaryDivisionRule
virtual void AcceptCellWriter(boost::shared_ptr< AbstractCellWriter< DIM, DIM > > pCellWriter, CellPtr pCell)
virtual double GetCellDataItemAtPdeNode(unsigned pdeNodeIndex, std::string &rVariableName, bool dirichletBoundaryConditionApplies=false, double dirichletBoundaryValue=0.0)
void OutputCellPopulationParameters(out_stream &rParamsFile)
virtual void AcceptPopulationEventWriter(boost::shared_ptr< AbstractCellPopulationEventWriter< DIM, DIM > > pPopulationEventWriter)
ImmersedBoundaryMesh< DIM, DIM > & rGetMesh()
boost::shared_ptr< AbstractImmersedBoundaryDivisionRule< DIM > > GetImmersedBoundaryDivisionRule()
ImmersedBoundaryElement< DIM - 1, DIM > * GetLamina(unsigned laminaIndex)
ImmersedBoundaryElement< DIM, DIM > * GetElementCorrespondingToCell(CellPtr pCell)
ImmersedBoundaryElement< DIM, DIM > * GetElement(unsigned elementIndex)
std::set< unsigned > GetNeighbouringNodeIndices(unsigned index)
c_vector< double, DIM > GetLocationOfCellCentre(CellPtr pCell)
void Update(bool hasHadBirthsOrDeaths=true)
void SetCellRearrangementThreshold(double newThreshold)
virtual void AcceptPopulationWriter(boost::shared_ptr< AbstractCellPopulationWriter< DIM, DIM > > pPopulationWriter)
virtual void AcceptPopulationCountWriter(boost::shared_ptr< AbstractCellPopulationCountWriter< DIM, DIM > > pPopulationCountWriter)
virtual bool IsPdeNodeAssociatedWithNonApoptoticCell(unsigned pdeNodeIndex)
void SetImmersedBoundaryDivisionRule(boost::shared_ptr< AbstractImmersedBoundaryDivisionRule< DIM > > pImmersedBoundaryDivisionRule)
void SetOutputNodeRegionToVtk(bool outputNodeRegionsToVtk)
CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell=CellPtr())
virtual void WriteVtkResultsToFile(const std::string &rDirectory)
std::set< unsigned > GetNeighbouringLocationIndices(CellPtr pCell)
void SetNode(unsigned index, ChastePoint< DIM > &rNewLocation)
virtual TetrahedralMesh< DIM, DIM > * GetTetrahedralMeshForPdeModifier()
ImmersedBoundaryMesh< DIM, DIM > * mpImmersedBoundaryMesh
void SetIfPopulationHasActiveSources(bool hasActiveSources)
virtual void CheckForStepSizeException(unsigned nodeIndex, c_vector< double, DIM > &rDisplacement, double dt)
double GetWidth(const unsigned &rDimension)
void FindElementOverlaps(ImmersedBoundaryMesh< ELEMENT_DIM, SPACE_DIM > &rMesh)
const std::vector< std::vector< unsigned > > & rGetElementParts() const
void WriteVtkUsingMesh(ImmersedBoundaryMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, std::string stamp="")
void AddPointData(std::string dataName, std::vector< double > dataPayload)
void AddCellData(std::string dataName, std::vector< double > dataPayload)
void SetNeighbourDist(double neighbourDist)
void SetElementDivisionSpacing(double elementDivisionSpacing)
Definition Node.hpp:59
std::set< unsigned > & rGetContainingElementIndices()
Definition Node.cpp:300
std::string GetOutputDirectoryFullPath() const
FileFinder FindFile(std::string leafName) const
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
static bool IsSequential()
static SimulationTime * Instance()
unsigned GetTimeStepsElapsed() const
void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)