36#include "ImmersedBoundaryCellPopulation.hpp"
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"
48template <
unsigned DIM>
51 std::vector<CellPtr>& rCells,
54 const std::vector<unsigned> 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)
71 std::list<CellPtr>::iterator it = this->
mCells.begin();
72 for (
unsigned i = 0; it != this->
mCells.end(); ++it, ++i)
75 unsigned index = locationIndices.empty() ? i : locationIndices[i];
92template <
unsigned DIM>
97 mIntrinsicSpacing(0.01),
98 mPopulationHasActiveSources(false),
99 mOutputNodeRegionToVtk(false),
100 mReMeshFrequency(UINT_MAX),
101 mThrowStepSizeException(true),
102 mCellRearrangementThreshold(0.5)
107template <
unsigned DIM>
112 delete &this->mrMesh;
116template <
unsigned DIM>
119 double average_intrinsic_size = 0.0;
121 for (
auto elem_iter = mpImmersedBoundaryMesh->GetElementIteratorBegin();
122 elem_iter != mpImmersedBoundaryMesh->GetElementIteratorEnd();
125 average_intrinsic_size += mpImmersedBoundaryMesh->GetVolumeOfElement(elem_iter->GetIndex());
128 return sqrt(average_intrinsic_size / mpImmersedBoundaryMesh->GetNumElements());
131template <
unsigned DIM>
137template <
unsigned DIM>
140 return *mpImmersedBoundaryMesh;
143template <
unsigned DIM>
146 return *mpImmersedBoundaryMesh;
149template <
unsigned DIM>
152 return mpImmersedBoundaryMesh->GetElement(elementIndex);
155template <
unsigned DIM>
158 return mpImmersedBoundaryMesh->GetLamina(laminaIndex);
161template <
unsigned DIM>
164 return this->mrMesh.GetNumNodes();
167template <
unsigned DIM>
170 return mpImmersedBoundaryMesh->GetCentroidOfElement(this->mCellLocationMap[pCell.get()]);
173template <
unsigned DIM>
176 return this->mrMesh.GetNode(index);
179template <
unsigned DIM>
182 assert(newDistance >= 0.0);
183 mInteractionDistance = newDistance;
186template <
unsigned DIM>
189 return mInteractionDistance;
192template <
unsigned DIM>
195 mReMeshFrequency = newFrequency;
198template <
unsigned DIM>
201 return mReMeshFrequency;
204template <
unsigned DIM>
207 mCellRearrangementThreshold = newThreshold;
210template <
unsigned DIM>
213 return mCellRearrangementThreshold;
216template <
unsigned DIM>
219 mThrowStepSizeException =
throws;
222template <
unsigned DIM>
225 return mThrowStepSizeException;
228template <
unsigned DIM>
231 return mIntrinsicSpacing;
234template <
unsigned DIM>
238 return mpImmersedBoundaryMesh->GetNeighbouringElementIndices(this->GetLocationIndexUsingCell(pCell));
241template <
unsigned DIM>
244 if (pNewNode ==
nullptr)
248 return mpImmersedBoundaryMesh->AddNode(pNewNode);
251template <
unsigned DIM>
254 mpImmersedBoundaryMesh->SetNode(nodeIndex, rNewLocation);
257template <
unsigned DIM>
261 return mpImmersedBoundaryMesh->GetElement(this->GetLocationIndexUsingCell(pCell));
264template <
unsigned DIM>
267 return mpImmersedBoundaryMesh->GetNumElements();
270template <
unsigned DIM>
273 return mpImmersedBoundaryMesh->GetNumLaminas();
276template <
unsigned DIM>
285 c_vector<double, DIM> division_vector = mpImmersedBoundaryDivisionRule->CalculateCellDivisionVector(pParentCell, *
this);
288 unsigned new_elem_idx = mpImmersedBoundaryMesh->DivideElementAlongGivenAxis(p_element, division_vector,
true);
291 this->mCells.push_back(pNewCell);
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;
298 return p_created_cell;
301template <
unsigned DIM>
304 unsigned num_removed = 0;
306 for (std::list<CellPtr>::iterator it = this->mCells.begin();
307 it != this->mCells.end();
315 if (!(this->GetElement(this->GetLocationIndexUsingCell((*it)))->IsDeleted()))
317 this->GetElement(this->GetLocationIndexUsingCell(*it))->MarkAsDeleted();
321 it = this->mCells.erase(it);
331template <
unsigned DIM>
333 [[maybe_unused]]
double dt)
335 if constexpr (DIM == 2)
338 unsigned num_grid_pts_x = this->rGetMesh().GetNumGridPtsX();
339 unsigned num_grid_pts_y = this->rGetMesh().GetNumGridPtsY();
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;
345 unsigned first_idx_x;
346 unsigned first_idx_y;
348 std::vector<unsigned> x_indices(4);
349 std::vector<unsigned> y_indices(4);
351 std::vector<double> x_deltas(4);
352 std::vector<double> y_deltas(4);
356 c_vector<double, DIM> displacement = zero_vector<double>(DIM);
359 const multi_array<double, 3>& vel_grids = this->rGetMesh().rGet2dVelocityGrids();
362 c_vector<double, DIM> node_location;
363 for (
auto node_iter = this->rGetMesh().GetNodeIteratorBegin(
false);
364 node_iter != this->rGetMesh().GetNodeIteratorEnd();
368 node_location = node_iter->rGetLocation();
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;
375 for (
unsigned i = 0; i < 4; i++)
377 x_indices[i] = (first_idx_x + i) % num_grid_pts_x;
378 y_indices[i] = (first_idx_y + i) % num_grid_pts_y;
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);
385 for (
unsigned x_idx = 0; x_idx < 4; ++x_idx)
387 for (
unsigned y_idx = 0; y_idx < 4; ++y_idx)
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;
400 if (norm_2(displacement) > characteristic_spacing)
402 if (norm_2(displacement) > 10.0 * characteristic_spacing)
404 EXCEPTION(
"Nodes are moving more than 10x CharacteristicNodeSpacing. Aborting.");
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);
411 CheckForStepSizeException(node_iter->GetIndex(), displacement, dt);
414 node_location += displacement;
417 for (
unsigned i = 0; i < DIM; ++i)
419 node_location[i] = fmod(node_location[i] + 1.0, 1.0);
426 this->SetNode(node_iter->GetIndex(), new_point);
430 if (this->DoesPopulationHaveActiveSources())
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();
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());
440 c_vector<double, DIM> source_location;
443 for (
unsigned source_idx = 0; source_idx < combined_sources.size(); source_idx++)
446 source_location = combined_sources[source_idx]->rGetLocation();
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;
453 for (
unsigned i = 0; i < 4; ++i)
455 x_indices[i] = (first_idx_x + i) % num_grid_pts_x;
456 y_indices[i] = (first_idx_y + i) % num_grid_pts_y;
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);
463 for (
unsigned x_idx = 0; x_idx < 4; ++x_idx)
465 for (
unsigned y_idx = 0; y_idx < 4; ++y_idx)
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;
478 if (norm_2(displacement) > characteristic_spacing)
480 if (norm_2(displacement) > 10.0 * characteristic_spacing)
482 EXCEPTION(
"Sources are moving more than 10x CharacteristicNodeSpacing. Aborting.");
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);
490 source_location += displacement;
493 for (
unsigned i = 0; i < DIM; ++i)
495 source_location[i] = fmod(source_location[i] + 1.0, 1.0);
499 combined_sources[source_idx]->rGetModifiableLocation() = source_location;
505 if (num_time_steps > 0 && num_time_steps % mReMeshFrequency == 0)
507 mpImmersedBoundaryMesh->ReMesh();
516template <
unsigned DIM>
519 return (0.25 * (1.0 + cos(M_PI * dist / (2 * spacing))));
522template <
unsigned DIM>
525 return GetElementCorrespondingToCell(pCell)->IsDeleted();
528template <
unsigned DIM>
532 if (this->Begin()->GetCellData()->HasItem(
"target area"))
534 for (
auto cell_iter = this->Begin();
535 cell_iter != this->End();
538 double target_area = cell_iter->GetCellData()->GetItem(
"target area");
539 double actual_area = this->GetVolumeOfCell(*cell_iter);
541 double strength = 1e-2 * (target_area - actual_area) / target_area;
543 this->GetElementCorrespondingToCell(*cell_iter)->GetFluidSource()->SetStrength(strength);
548template <
unsigned DIM>
552 std::vector<unsigned> validated_element = std::vector<unsigned>(this->GetNumElements(), 0);
553 for (
auto cell_iter = this->Begin();
554 cell_iter != this->End();
557 unsigned elem_index = this->GetLocationIndexUsingCell(*cell_iter);
558 validated_element[elem_index]++;
561 for (
unsigned i = 0; i < validated_element.size(); ++i)
563 if (validated_element[i] == 0)
568 if (validated_element[i] > 1)
572 EXCEPTION(
"At time " <<
SimulationTime::Instance()->GetTime() <<
", Element " << i <<
" appears to have " << validated_element[i] <<
" cells associated with it");
577template <
unsigned DIM>
580 c_vector<double, DIM>& rDisplacement,
583 double length = boost::numeric::ublas::norm_2(rDisplacement);
585 if (length > 0.5*mpImmersedBoundaryMesh->GetCellRearrangementThreshold())
587 rDisplacement *= 0.5*mpImmersedBoundaryMesh->GetCellRearrangementThreshold()/length;
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.";
593 double suggested_step = 0.95*dt*((0.5*mpImmersedBoundaryMesh->GetCellRearrangementThreshold())/length);
596 if(mThrowStepSizeException)
598 mThrowStepSizeException =
false;
604template <
unsigned DIM>
608 pPopulationWriter->Visit(
this);
611template <
unsigned DIM>
615 pPopulationEventWriter->Visit(
this);
618template <
unsigned DIM>
622 pPopulationCountWriter->Visit(
this);
625template <
unsigned DIM>
629 pCellWriter->VisitCell(pCell,
this);
632template <
unsigned DIM>
636 unsigned elem_index = this->GetLocationIndexUsingCell(pCell);
639 double cell_volume = mpImmersedBoundaryMesh->GetVolumeOfElement(elem_index);
644template <
unsigned DIM>
646 const std::string& rDirectory)
654 const std::vector<std::vector<unsigned>>& r_elem_parts = mesh_writer.
rGetElementParts();
657 for (
auto cell_writer_iter = this->mCellWriters.begin();
658 cell_writer_iter != this->mCellWriters.end();
662 std::vector<double> vtk_cell_data;
665 for (
auto elem_iter = mpImmersedBoundaryMesh->GetElementIteratorBegin();
666 elem_iter != mpImmersedBoundaryMesh->GetElementIteratorEnd();
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();
677 CellPtr p_cell = this->GetCellUsingLocationIndex(elem_index);
685 for (
unsigned elem_part = 0; elem_part < num_elem_parts; ++elem_part)
687 vtk_cell_data.push_back((*cell_writer_iter)->GetCellDataForVtkOutput(p_cell,
this));
695 for (
auto lam_iter = mpImmersedBoundaryMesh->GetLaminaIteratorBegin();
696 lam_iter != mpImmersedBoundaryMesh->GetLaminaIteratorEnd();
699 vtk_cell_data.push_back(-1.0);
702 mesh_writer.
AddCellData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
709 const unsigned num_cell_data_items = this->Begin()->GetCellData()->GetNumItems();
710 std::vector<std::string> cell_data_names = this->Begin()->GetCellData()->GetKeys();
712 std::vector<std::vector<double>> cell_data;
713 for (
unsigned var = 0; var < num_cell_data_items; ++var)
715 std::vector<double> cell_data_var;
716 cell_data.push_back(cell_data_var);
720 for (
auto elem_iter = mpImmersedBoundaryMesh->GetElementIteratorBegin();
721 elem_iter != mpImmersedBoundaryMesh->GetElementIteratorEnd();
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();
732 CellPtr p_cell = this->GetCellUsingLocationIndex(elem_index);
735 for (
unsigned var = 0; var < num_cell_data_items; var++)
742 for (
unsigned elem_part = 0; elem_part < num_elem_parts; ++elem_part)
744 cell_data[var].push_back(p_cell->GetCellData()->GetItem(cell_data_names[var]));
753 for (
auto lam_iter = mpImmersedBoundaryMesh->GetLaminaIteratorBegin();
754 lam_iter != mpImmersedBoundaryMesh->GetLaminaIteratorEnd();
757 for (
unsigned var = 0; var < num_cell_data_items; ++var)
763 for (
unsigned var = 0; var < num_cell_data_items; ++var)
765 mesh_writer.
AddCellData(cell_data_names[var], cell_data[var]);
769 if (mOutputNodeRegionToVtk)
771 std::vector<double> node_regions;
772 for (
auto node_iter = mpImmersedBoundaryMesh->GetNodeIteratorBegin();
773 node_iter != mpImmersedBoundaryMesh->GetNodeIteratorEnd();
776 node_regions.push_back(
static_cast<double>(node_iter->GetRegion()));
782 std::stringstream time;
783 time << num_timesteps;
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";
795template <
unsigned DIM>
799 if (this->mOutputResultsForChasteVisualizer)
801 if (!this->
template HasWriter<CellPopulationElementWriter>())
803 this->
template AddPopulationWriter<CellPopulationElementWriter>();
810template <
unsigned DIM>
812 out_stream& rParamsFile)
815 *rParamsFile <<
"\t\t<ImmersedBoundaryDivisionRule>\n";
816 mpImmersedBoundaryDivisionRule->OutputCellImmersedBoundaryDivisionRuleInfo(rParamsFile);
817 *rParamsFile <<
"\t\t</ImmersedBoundaryDivisionRule>\n";
823template <
unsigned DIM>
826 double width = this->mrMesh.GetWidth(rDimension);
830template <
unsigned DIM>
834 return mpImmersedBoundaryMesh->GetNeighbouringNodeIndices(index);
837template <
unsigned DIM>
842 if constexpr (DIM == 2)
844 unsigned num_vertex_nodes = mpImmersedBoundaryMesh->GetNumNodes();
845 unsigned num_vertex_elements = mpImmersedBoundaryMesh->GetNumElements();
847 std::string mesh_file_name =
"mesh";
850 std::stringstream pid;
852 OutputFileHandler output_file_handler(
"2D_temporary_tetrahedral_mesh_" + pid.str());
856 unsigned num_tetrahedral_nodes = num_vertex_nodes + num_vertex_elements;
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;
865 auto nodes = mpImmersedBoundaryMesh->rGetNodes();
866 for (
auto p_node : nodes)
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;
872 (*p_node_file) << index <<
"\t" << r_location[0] <<
"\t" << r_location[1] <<
"\t" << is_boundary_node << std::endl;
876 unsigned num_tetrahedral_elements = 0;
877 for (
unsigned vertex_elem_index = 0;
878 vertex_elem_index < num_vertex_elements;
881 unsigned index = num_vertex_nodes + vertex_elem_index;
883 c_vector<double, DIM> location = mpImmersedBoundaryMesh->GetCentroidOfElement(vertex_elem_index);
886 unsigned is_boundary_node = 0;
887 (*p_node_file) << index <<
"\t" << location[0] <<
"\t" << location[1] <<
"\t" << is_boundary_node << std::endl;
890 num_tetrahedral_elements += mpImmersedBoundaryMesh->GetElement(vertex_elem_index)->GetNumNodes();
892 p_node_file->close();
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;
899 std::set<std::pair<unsigned, unsigned> > tetrahedral_edges;
901 unsigned tetrahedral_elem_index = 0;
902 for (
unsigned vertex_elem_index = 0;
903 vertex_elem_index < num_vertex_elements;
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;
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;
918 (*p_elem_file) << tetrahedral_elem_index++ <<
"\t" << node_0_index <<
"\t" << node_1_index <<
"\t" << node_2_index << std::endl;
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);
925 tetrahedral_edges.insert(edge_0);
926 tetrahedral_edges.insert(edge_1);
927 tetrahedral_edges.insert(edge_2);
930 p_elem_file->close();
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;
937 unsigned edge_index = 0;
938 for (
auto edge_iter = tetrahedral_edges.begin();
939 edge_iter != tetrahedral_edges.end();
942 std::pair<unsigned, unsigned> this_edge = *edge_iter;
945 bool is_boundary_edge =
false;
946 if (this_edge.first < mpImmersedBoundaryMesh->GetNumNodes() &&
947 this_edge.second < mpImmersedBoundaryMesh->GetNumNodes())
949 is_boundary_edge = (mpImmersedBoundaryMesh->GetNode(this_edge.first)->IsBoundaryNode() &&
950 mpImmersedBoundaryMesh->GetNode(this_edge.second)->IsBoundaryNode() );
952 unsigned is_boundary_edge_unsigned = is_boundary_edge ? 1 : 0;
954 (*p_edge_file) << edge_index++ <<
"\t" << this_edge.first <<
"\t" << this_edge.second <<
"\t" << is_boundary_edge_unsigned << std::endl;
956 p_edge_file->close();
980 EXCEPTION(
"ImmersedBoundaryCellPopulation::GetTetrahedralMeshForPDEModifier is only implemented in 2D");
984template <
unsigned DIM>
987 bool non_apoptotic_cell_present =
true;
989 if (pdeNodeIndex < this->GetNumNodes())
991 std::set<unsigned> containing_element_indices = this->GetNode(pdeNodeIndex)->rGetContainingElementIndices();
993 for (
auto iter = containing_element_indices.begin();
994 iter != containing_element_indices.end();
997 if (this->GetCellUsingLocationIndex(*iter)->template HasCellProperty<ApoptoticCellProperty>() )
999 non_apoptotic_cell_present =
false;
1011 non_apoptotic_cell_present = !(this->GetCellUsingLocationIndex(pdeNodeIndex - this->GetNumNodes())->template HasCellProperty<ApoptoticCellProperty>());
1014 return non_apoptotic_cell_present;
1017template <
unsigned DIM>
1019 unsigned pdeNodeIndex,
1020 std::string& rVariableName,
1021 bool dirichletBoundaryConditionApplies,
1022 double dirichletBoundaryValue)
1024 unsigned num_nodes = this->GetNumNodes();
1031 if (pdeNodeIndex >= num_nodes)
1034 assert(pdeNodeIndex-num_nodes < num_nodes);
1036 CellPtr p_cell = this->GetCellUsingLocationIndex(pdeNodeIndex - num_nodes);
1037 value = p_cell->GetCellData()->GetItem(rVariableName);
1042 if (dirichletBoundaryConditionApplies)
1045 value = dirichletBoundaryValue;
1049 assert(pdeNodeIndex < num_nodes);
1050 Node<DIM>* p_node = this->GetNode(pdeNodeIndex);
1054 for (
auto index_iter = containing_elements.begin();
1055 index_iter != containing_elements.end();
1058 assert(*index_iter < num_nodes);
1059 CellPtr p_cell = this->GetCellUsingLocationIndex(*index_iter);
1060 value += p_cell->GetCellData()->GetItem(rVariableName);
1062 value /= containing_elements.size();
1069template <
unsigned DIM>
1072 return mpImmersedBoundaryDivisionRule;
1075template <
unsigned DIM>
1079 mpImmersedBoundaryDivisionRule = pImmersedBoundaryDivisionRule;
1082template <
unsigned DIM>
1085 return mPopulationHasActiveSources;
1088template <
unsigned DIM>
1091 return this->GetElementCorrespondingToCell(pCell)->IsElementOnBoundary();
1094template <
unsigned DIM>
1096 bool hasActiveSources)
1098 mPopulationHasActiveSources = hasActiveSources;
1101template <
unsigned DIM>
1103 bool outputNodeRegionsToVtk)
1105 mOutputNodeRegionToVtk = outputNodeRegionsToVtk;
1108template <
unsigned DIM>
const double DOUBLE_UNSET
#define EXCEPTION(message)
#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
std::list< CellPtr > mCells
unsigned GetNumNodes() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
void SetMeshHasChangedSinceLoading()
virtual void OutputCellPopulationParameters(out_stream &rParamsFile)
virtual void OpenWritersFiles(OutputFileHandler &rOutputFileHandler)
double GetCellRearrangementThreshold() const
bool IsCellOnBoundary(CellPtr pCell)
void SetReMeshFrequency(unsigned newFrequency)
ImmersedBoundaryCellPopulation(ImmersedBoundaryMesh< DIM, DIM > &rMesh, std::vector< CellPtr > &rCells, bool deleteMesh=false, bool validate=true, const std::vector< unsigned > locationIndices=std::vector< unsigned >())
unsigned RemoveDeadCells()
double Delta1D(double dist, double spacing)
double mInteractionDistance
boost::shared_ptr< AbstractImmersedBoundaryDivisionRule< DIM > > mpImmersedBoundaryDivisionRule
double GetIntrinsicSpacing() const
unsigned GetReMeshFrequency() const
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)
double GetInteractionDistance() const
virtual void AcceptPopulationEventWriter(boost::shared_ptr< AbstractCellPopulationEventWriter< DIM, DIM > > pPopulationEventWriter)
ImmersedBoundaryMesh< DIM, DIM > & rGetMesh()
double CalculateIntrinsicCellSize()
boost::shared_ptr< AbstractImmersedBoundaryDivisionRule< DIM > > GetImmersedBoundaryDivisionRule()
unsigned GetNumElements()
ImmersedBoundaryElement< DIM - 1, DIM > * GetLamina(unsigned laminaIndex)
virtual ~ImmersedBoundaryCellPopulation()
ImmersedBoundaryElement< DIM, DIM > * GetElementCorrespondingToCell(CellPtr pCell)
unsigned AddNode(Node< DIM > *pNewNode)
bool DoesPopulationHaveActiveSources() const
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)
void UpdateNodeLocations(double dt)
bool ThrowsStepSizeException() const
bool IsCellAssociatedWithADeletedLocation(CellPtr pCell)
double GetVolumeOfCell(CellPtr pCell)
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 SetThrowsStepSizeException(bool throws)
void SetImmersedBoundaryDivisionRule(boost::shared_ptr< AbstractImmersedBoundaryDivisionRule< DIM > > pImmersedBoundaryDivisionRule)
void SetOutputNodeRegionToVtk(bool outputNodeRegionsToVtk)
virtual double GetDefaultTimeStep()
CellPtr AddCell(CellPtr pNewCell, CellPtr pParentCell=CellPtr())
virtual void WriteVtkResultsToFile(const std::string &rDirectory)
double GetDampingConstant(unsigned nodeIndex)
Node< DIM > * GetNode(unsigned index)
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 SetInteractionDistance(double newDistance)
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)
std::set< unsigned > & rGetContainingElementIndices()
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 SimulationTime * Instance()
unsigned GetTimeStepsElapsed() const
void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)