289 if (!(this->mUpdateRuleCollection.empty()))
293 for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
294 cell_iter != this->mCells.end();
298 unsigned node_index = this->GetLocationIndexUsingCell(*cell_iter);
302 std::vector<double> neighbouring_node_propensities;
303 std::vector<unsigned> neighbouring_node_indices_vector;
305 if (!neighbouring_node_indices.empty())
307 unsigned num_neighbours = neighbouring_node_indices.size();
308 double probability_of_not_moving = 1.0;
310 for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
311 iter != neighbouring_node_indices.end();
314 double probability_of_moving = 0.0;
316 neighbouring_node_indices_vector.push_back(*iter);
318 if (IsSiteAvailable(*iter, *cell_iter))
321 for (
typename std::vector<boost::shared_ptr<
AbstractUpdateRule<DIM> > >::iterator iter_rule = this->mUpdateRuleCollection.begin();
322 iter_rule != this->mUpdateRuleCollection.end();
326 double p = (boost::static_pointer_cast<AbstractCaUpdateRule<DIM> >(*iter_rule))->EvaluateProbability(node_index, *iter, *
this, dt, 1, *cell_iter);
327 probability_of_moving += p;
328 if (probability_of_moving < 0)
330 EXCEPTION(
"The probability of cellular movement is smaller than zero. In order to prevent it from happening you should change your time step and parameters");
333 if (probability_of_moving > 1)
335 EXCEPTION(
"The probability of the cellular movement is bigger than one. In order to prevent it from happening you should change your time step and parameters");
339 probability_of_not_moving -= probability_of_moving;
340 neighbouring_node_propensities.push_back(probability_of_moving);
344 neighbouring_node_propensities.push_back(0.0);
347 if (probability_of_not_moving < 0)
349 EXCEPTION(
"The probability of the cell not moving is smaller than zero. In order to prevent it from happening you should change your time step and parameters");
354 double random_number = p_gen->
ranf();
356 double total_probability = 0.0;
357 for (
unsigned counter=0; counter<num_neighbours; counter++)
359 total_probability += neighbouring_node_propensities[counter];
360 if (total_probability >= random_number)
363 unsigned chosen_neighbour_location_index = neighbouring_node_indices_vector[counter];
364 this->MoveCellInLocationMap((*cell_iter), node_index, chosen_neighbour_location_index);
384 if (!(mSwitchingUpdateRuleCollection.empty()))
386 assert(mLatticeCarryingCapacity == 1);
389 unsigned num_nodes = this->mrMesh.GetNumNodes();
392 if (this->mIterateRandomlyOverUpdateRuleCollection)
395 p_gen->
Shuffle(mSwitchingUpdateRuleCollection);
398 for (
unsigned i=0; i<num_nodes; i++)
402 if (this->mUpdateNodesInRandomOrder)
404 node_index = p_gen->
randMod(num_nodes);
409 node_index = i%num_nodes;
415 unsigned neighbour_location_index;
417 if (!neighbouring_node_indices.empty())
419 unsigned num_neighbours = neighbouring_node_indices.size();
420 unsigned chosen_neighbour = p_gen->
randMod(num_neighbours);
422 std::set<unsigned>::iterator neighbour_iter = neighbouring_node_indices.begin();
423 for (
unsigned j=0; j<chosen_neighbour; j++)
427 neighbour_location_index = *neighbour_iter;
429 bool is_cell_on_node_index = mAvailableSpaces[node_index] == 0 ? true :
false;
430 bool is_cell_on_neighbour_location_index = mAvailableSpaces[neighbour_location_index] == 0 ? true :
false;
432 if (is_cell_on_node_index || is_cell_on_neighbour_location_index)
434 double probability_of_switch = 0.0;
437 for (
typename std::vector<boost::shared_ptr<
AbstractUpdateRule<DIM> > >::iterator iter_rule = mSwitchingUpdateRuleCollection.begin();
438 iter_rule != mSwitchingUpdateRuleCollection.end();
442 double p = (boost::static_pointer_cast<AbstractCaSwitchingUpdateRule<DIM> >(*iter_rule))->EvaluateSwitchingProbability(node_index, neighbour_location_index, *
this, dt, 1);
443 probability_of_switch += p;
446 assert(probability_of_switch >= 0);
447 assert(probability_of_switch <= 1);
450 double random_number = p_gen->
ranf();
452 if (random_number < probability_of_switch)
454 if (is_cell_on_node_index && is_cell_on_neighbour_location_index)
457 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
458 CellPtr p_neighbour_cell = this->GetCellUsingLocationIndex(neighbour_location_index);
461 RemoveCellUsingLocationIndex(node_index, p_cell);
462 RemoveCellUsingLocationIndex(neighbour_location_index, p_neighbour_cell);
465 AddCellUsingLocationIndex(node_index, p_neighbour_cell);
466 AddCellUsingLocationIndex(neighbour_location_index, p_cell);
468 else if (is_cell_on_node_index && !is_cell_on_neighbour_location_index)
471 CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
472 RemoveCellUsingLocationIndex(node_index, p_cell);
473 AddCellUsingLocationIndex(neighbour_location_index, p_cell);
475 else if (!is_cell_on_node_index && is_cell_on_neighbour_location_index)
478 CellPtr p_neighbour_cell = this->GetCellUsingLocationIndex(neighbour_location_index);
479 RemoveCellUsingLocationIndex(neighbour_location_index, p_neighbour_cell);
480 AddCellUsingLocationIndex(node_index, p_neighbour_cell);
630 std::stringstream time;
631 time << num_timesteps;
634 unsigned num_cells = this->GetNumRealCells();
637 unsigned num_cell_data_items = 0u;
638 std::vector<std::string> cell_data_names;
641 num_cell_data_items = this->Begin()->GetCellData()->GetNumItems();
642 cell_data_names = this->Begin()->GetCellData()->GetKeys();
644 std::vector<std::vector<double> > cell_data;
645 for (
unsigned var=0; var<num_cell_data_items; var++)
647 std::vector<double> cell_data_var(num_cells);
648 cell_data.push_back(cell_data_var);
655 unsigned num_sites = this->mrMesh.GetNumNodes();
656 boost::scoped_array<unsigned> number_of_cells_at_site(
new unsigned[num_sites]);
657 for (
unsigned i=0; i<num_sites; i++)
659 number_of_cells_at_site[i] = 0;
663 std::vector<Node<DIM>*> nodes;
664 unsigned node_index = 0;
666 cell_iter != this->End();
670 unsigned location_index = this->GetLocationIndexUsingCell(*cell_iter);
671 number_of_cells_at_site[location_index]++;
672 assert(number_of_cells_at_site[location_index] <= mLatticeCarryingCapacity);
675 c_vector<double, DIM> coords;
676 coords = this->mrMesh.GetNode(location_index)->rGetLocation();
679 if (mLatticeCarryingCapacity > 1)
681 c_vector<double, DIM> offset;
685 double angle = (
double)number_of_cells_at_site[location_index]*2.0*M_PI/(
double)mLatticeCarryingCapacity;
686 offset[0] = 0.2*sin(angle);
687 offset[1] = 0.2*cos(angle);
692 for (
unsigned i=0; i<DIM; i++)
694 offset[i] = p_gen->
ranf();
698 for (
unsigned i=0; i<DIM; i++)
700 coords[i] += offset[i];
704 nodes.push_back(
new Node<DIM>(node_index, coords,
false));
710 cell_writer_iter != this->mCellWriters.end();
715 std::vector<double> vtk_cell_data(num_cells, -1.0);
718 unsigned cell_index = 0;
720 cell_iter != this->End();
724 vtk_cell_data[cell_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(*cell_iter,
this);
728 mesh_writer.
AddPointData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
732 unsigned cell_index = 0;
734 cell_iter != this->End();
737 for (
unsigned var=0; var<num_cell_data_items; var++)
739 cell_data[var][cell_index] = cell_iter->GetCellData()->GetItem(cell_data_names[var]);
743 for (
unsigned var=0; var<num_cell_data_items; var++)
745 mesh_writer.
AddPointData(cell_data_names[var], cell_data[var]);
757 double volume = this->mrMesh.
GetWidth(0);
758 for (
unsigned idx=1; idx<DIM; idx++)
760 volume *= this->mrMesh.GetWidth(idx);
764 if (this->mrMesh.GetNumNodes() >0 && volume > 0.0)
766 spacing = std::pow(volume /
double(this->mrMesh.GetNumNodes()), 1.0/
double(DIM));
776 *(this->mpVtkMetaFile) <<
" <DataSet timestep=\"";
777 *(this->mpVtkMetaFile) << num_timesteps;
778 *(this->mpVtkMetaFile) <<
"\" group=\"\" part=\"0\" file=\"results_";
779 *(this->mpVtkMetaFile) << num_timesteps;
780 *(this->mpVtkMetaFile) <<
".vtu\"/>\n";
783 for (
unsigned i=0; i<nodes.size(); i++)