127 if (!this->mSetBcsOnBoxBoundary)
130 if (this->mSetBcsOnBoundingSphere)
133 c_vector<double, DIM> tissue_maxima = zero_vector<double>(DIM);
134 c_vector<double, DIM> tissue_minima = zero_vector<double>(DIM);
135 for (
unsigned i = 0; i < DIM; i++)
137 tissue_maxima[i] = -DBL_MAX;
138 tissue_minima[i] = DBL_MAX;
143 cell_iter != rCellPopulation.
End();
148 for (
unsigned i = 0; i < DIM; i++)
150 if (r_position_of_cell[i] > tissue_maxima[i])
152 tissue_maxima[i] = r_position_of_cell[i];
154 if (r_position_of_cell[i] < tissue_minima[i])
156 tissue_minima[i] = r_position_of_cell[i];
161 c_vector<double, DIM> tissue_centre = 0.5*(tissue_maxima + tissue_minima);
164 double tissue_radius = 0.0;
166 cell_iter != rCellPopulation.
End();
171 double radius = norm_2(tissue_centre - r_position_of_cell.
rGetLocation());
173 if (tissue_radius < radius)
175 tissue_radius = radius;
180 if (this->IsNeumannBoundaryCondition())
187 for (
unsigned i=0; i<this->mpFeMesh->GetNumNodes(); i++)
189 double radius = norm_2(tissue_centre - this->mpFeMesh->GetNode(i)->rGetLocation());
191 if (radius > tissue_radius)
193 pBcc->AddDirichletBoundaryCondition(this->mpFeMesh->GetNode(i), this->mpBoundaryCondition.get(), 0,
false);
194 this->mIsDirichletBoundaryNode[i] = 1.0;
202 std::set<unsigned> coarse_element_indices_in_map;
204 cell_iter != rCellPopulation.
End();
207 coarse_element_indices_in_map.insert(this->mCellPdeElementMap[*cell_iter]);
211 std::set<unsigned> coarse_mesh_boundary_node_indices;
212 for (
unsigned i=0; i<this->mpFeMesh->GetNumElements(); i++)
214 if (coarse_element_indices_in_map.find(i) == coarse_element_indices_in_map.end())
217 for (
unsigned j=0; j<DIM+1; j++)
220 coarse_mesh_boundary_node_indices.insert(node_index);
226 std::set<unsigned> nearby_node_indices;
227 for (
unsigned int coarse_mesh_boundary_node_idx : coarse_mesh_boundary_node_indices)
229 bool remove_node =
false;
231 c_vector<double,DIM> node_location = this->mpFeMesh->GetNode(coarse_mesh_boundary_node_idx)->rGetLocation();
234 cell_iter != rCellPopulation.
End();
239 double separation = norm_2(node_location - r_position_of_cell.
rGetLocation());
241 if (separation <= mTypicalCellRadius)
251 nearby_node_indices.insert(coarse_mesh_boundary_node_idx);
256 for (
unsigned int nearby_node_idx : nearby_node_indices)
258 coarse_mesh_boundary_node_indices.erase(nearby_node_idx);
263 if (this->IsNeumannBoundaryCondition())
270 for (
unsigned int coarse_mesh_boundary_node_idx : coarse_mesh_boundary_node_indices)
272 pBcc->AddDirichletBoundaryCondition(this->mpFeMesh->GetNode(coarse_mesh_boundary_node_idx), this->mpBoundaryCondition.get(), 0,
false);
273 this->mIsDirichletBoundaryNode[coarse_mesh_boundary_node_idx] = 1.0;
280 if (this->IsNeumannBoundaryCondition())
283 for (
auto elem_iter = this->mpFeMesh->GetBoundaryElementIteratorBegin();
284 elem_iter != this->mpFeMesh->GetBoundaryElementIteratorEnd();
287 pBcc->AddNeumannBoundaryCondition(*elem_iter, this->mpBoundaryCondition.get());
293 for (
auto node_iter = this->mpFeMesh->GetBoundaryNodeIteratorBegin();
294 node_iter != this->mpFeMesh->GetBoundaryNodeIteratorEnd();
297 pBcc->AddDirichletBoundaryCondition(*node_iter, this->mpBoundaryCondition.get());
298 this->mIsDirichletBoundaryNode[(*node_iter)->GetIndex()] = 1.0;
336 pMesh->
ConstructRegularSlabMesh(stepSize, pMeshCuboid->GetWidth(0), pMeshCuboid->GetWidth(1), pMeshCuboid->GetWidth(2));
348 c_vector<double,DIM> centre_of_coarse_mesh = zero_vector<double>(DIM);
351 centre_of_coarse_mesh += pMesh->
GetNode(i)->rGetLocation();
356 pMesh->
Translate(centre_of_cuboid - centre_of_coarse_mesh);
365 if (mUseVoronoiCellsForInterpolation)
367 unsigned num_nodes = rCellPopulation.
GetNumNodes();
369 std::vector<double> cell_data(num_nodes, -1);
370 std::vector<unsigned> num_cells(num_nodes, -1);
373 cell_iter != rCellPopulation.
End();
377 cell_data[cell_location_index]=0.0;
378 num_cells[cell_location_index]=0;
383 node_iter != this->mpFeMesh->GetNodeIteratorEnd();
386 unsigned node_index = node_iter->GetIndex();
388 c_vector<double,DIM> node_location = node_iter->rGetLocation();
390 double closest_separation = DBL_MAX;
394 cell_iter != rCellPopulation.
End();
399 double separation = norm_2(node_location - cell_location);
401 if (separation < closest_separation)
403 closest_separation = separation;
407 assert(closest_separation<DBL_MAX);
409 cell_data[nearest_cell] = cell_data[nearest_cell] + solution_repl[node_index];
410 num_cells[nearest_cell] = num_cells[nearest_cell] + 1;
415 cell_iter != rCellPopulation.
End();
421 if (num_cells[cell_location_index]==0)
423 EXCEPTION(
"One or more of the cells doesnt contain any pde nodes so cant use Voronoi CellData calculation in the ");
426 double solution_at_cell = cell_data[cell_location_index]/num_cells[cell_location_index];
428 cell_iter->GetCellData()->SetItem(this->mDependentVariableName, solution_at_cell);
431 if (this->mOutputGradient)
440 cell_iter != rCellPopulation.
End();
444 double solution_at_cell = 0.0;
447 unsigned elem_index = mCellPdeElementMap[*cell_iter];
454 for (
unsigned i=0; i<DIM+1; i++)
457 solution_at_cell += nodal_value * weights(i);
460 cell_iter->GetCellData()->SetItem(this->mDependentVariableName, solution_at_cell);
462 if (this->mOutputGradient)
465 c_vector<double, DIM> solution_gradient = zero_vector<double>(DIM);
468 c_matrix<double, DIM, DIM> jacobian, inverse_jacobian;
470 this->mpFeMesh->GetInverseJacobianForElement(elem_index, jacobian, jacobian_det, inverse_jacobian);
472 c_matrix<double, DIM, DIM+1> grad_phi;
475 for (
unsigned node_index=0; node_index<DIM+1; node_index++)
479 for (
unsigned j=0; j<DIM; j++)
481 solution_gradient(j) += nodal_value* grad_phi(j, node_index);
488 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+
"_grad_x", solution_gradient(0));
491 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+
"_grad_x", solution_gradient(0));
492 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+
"_grad_y", solution_gradient(1));
495 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+
"_grad_x", solution_gradient(0));
496 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+
"_grad_y", solution_gradient(1));
497 cell_iter->GetCellData()->SetItem(this->mDependentVariableName+
"_grad_z", solution_gradient(2));
AbstractBoxDomainPdeModifier(boost::shared_ptr< AbstractLinearPde< DIM, DIM > > pPde=boost::shared_ptr< AbstractLinearPde< DIM, DIM > >(), boost::shared_ptr< AbstractBoundaryCondition< DIM > > pBoundaryCondition=boost::shared_ptr< AbstractBoundaryCondition< DIM > >(), bool isNeumannBoundaryCondition=true, boost::shared_ptr< ChasteCuboid< DIM > > pMeshCuboid=boost::shared_ptr< ChasteCuboid< DIM > >(), double stepSize=1.0, Vec solution=nullptr)
void OutputSimulationModifierParameters(out_stream &rParamsFile) override
void GenerateFeMesh(boost::shared_ptr< ChasteCuboid< DIM > > pMeshCuboid, double stepSize)
void GenerateAndReturnFeMesh(boost::shared_ptr< ChasteCuboid< DIM > > pMeshCuboid, double stepSize, TetrahedralMesh< DIM, DIM > *pMesh)
void SetupSolve(AbstractCellPopulation< DIM, DIM > &rCellPopulation, std::string outputDirectory) override
boost::shared_ptr< ChasteCuboid< DIM > > mpMeshCuboid
void ConstructBoundaryConditionsContainerHelper(AbstractCellPopulation< DIM, DIM > &rCellPopulation, std::shared_ptr< BoundaryConditionsContainer< DIM, DIM, 1 > > pBcc)
void SetUseVoronoiCellsForInterpolation(bool useVoronoiCellsForInterpolation)
void UpdateCellData(AbstractCellPopulation< DIM, DIM > &rCellPopulation)
void UpdateCellPdeElementMap(AbstractCellPopulation< DIM, DIM > &rCellPopulation)
void InitialiseCellPdeElementMap(AbstractCellPopulation< DIM, DIM > &rCellPopulation)