void VolumeTrackingModifier<DIM>::UpdateCellData(AbstractCellPopulation<DIM,DIM>& rCellPopulation)
{
// Make sure the cell population is updated
rCellPopulation.Update();
/**
* This hack is needed because in the case of a MeshBasedCellPopulation in which
* multiple cell divisions have occurred over one time step, the Voronoi tessellation
* (while existing) is out-of-date. Thus, if we did not regenerate the Voronoi
* tessellation here, an assertion may trip as we try to access a Voronoi element
* whose index exceeds the number of elements in the out-of-date tessellation.
*
* \todo work out how to properly fix this (#1986)
*/
if (bool(dynamic_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation)))
{
static_cast<MeshBasedCellPopulation<DIM>*>(&(rCellPopulation))->CreateVoronoiTessellation();
}
// Iterate over cell population
for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = rCellPopulation.Begin();
cell_iter != rCellPopulation.End();
++cell_iter)
{
// Get the volume of this cell
double cell_volume = rCellPopulation.GetVolumeOfCell(*cell_iter);
// Store the cell's volume in CellData
cell_iter->GetCellData()->SetItem("volume", cell_volume);
}
}
void IsolatedLabelledCellKiller<DIM>::CheckAndLabelCellsForApoptosisOrDeath()
{
MutableVertexMesh<DIM, DIM>& vertex_mesh = static_cast<VertexBasedCellPopulation<DIM>*>(this->mpCellPopulation)->rGetMesh();
unsigned num_labelled_cells = this->mpCellPopulation->GetCellPropertyRegistry()->template Get<CellLabel>()->GetCellCount();
// If there is more than one labelled cell...
if (num_labelled_cells > 1)
{
// Iterate over cell population
for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mpCellPopulation->Begin();
cell_iter != this->mpCellPopulation->End();
++cell_iter)
{
// Only consider cells with the CellLabel property
if (cell_iter->template HasCellProperty<CellLabel>())
{
// Get the element index corresponding to this cell
unsigned elem_index = this->mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
// Get the set of neighbouring element indices
std::set<unsigned> neighbouring_elem_indices = vertex_mesh.GetNeighbouringElementIndices(elem_index);
// Check if any of the corresponding cells have the CellLabel property...
unsigned num_labelled_neighbours = 0;
for (std::set<unsigned>::iterator elem_iter = neighbouring_elem_indices.begin();
elem_iter != neighbouring_elem_indices.end();
++elem_iter)
{
if (this->mpCellPopulation->GetCellUsingLocationIndex(*elem_iter)->template HasCellProperty<CellLabel>())
{
num_labelled_neighbours++;
}
}
// ...and if none do, then kill this cell
if (num_labelled_neighbours == 0)
{
cell_iter->Kill();
}
}
}
}
}
void DeltaNotchTrackingModifier<DIM>::UpdateCellData(AbstractCellPopulation<DIM,DIM>& rCellPopulation)
{
// Make sure the cell population is updated
rCellPopulation.Update();
// First recover each cell's Notch and Delta concentrations from the ODEs and store in CellData
for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = rCellPopulation.Begin();
cell_iter != rCellPopulation.End();
++cell_iter)
{
DeltaNotchCellCycleModel* p_model = static_cast<DeltaNotchCellCycleModel*>(cell_iter->GetCellCycleModel());
double this_delta = p_model->GetDelta();
double this_notch = p_model->GetNotch();
// Note that the state variables must be in the same order as listed in DeltaNotchOdeSystem
cell_iter->GetCellData()->SetItem("notch", this_notch);
cell_iter->GetCellData()->SetItem("delta", this_delta);
}
// Next iterate over the population to compute and store each cell's neighbouring Delta concentration in CellData
for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = rCellPopulation.Begin();
cell_iter != rCellPopulation.End();
++cell_iter)
{
// Get the set of neighbouring location indices
std::set<unsigned> neighbour_indices = rCellPopulation.GetNeighbouringLocationIndices(*cell_iter);
// Compute this cell's average neighbouring Delta concentration and store in CellData
if (!neighbour_indices.empty())
{
double mean_delta = 0.0;
for (std::set<unsigned>::iterator iter = neighbour_indices.begin();
iter != neighbour_indices.end();
++iter)
{
CellPtr p_cell = rCellPopulation.GetCellUsingLocationIndex(*iter);
double this_delta = p_cell->GetCellData()->GetItem("delta");
mean_delta += this_delta/neighbour_indices.size();
}
cell_iter->GetCellData()->SetItem("mean delta", mean_delta);
}
else
{
// If this cell has no neighbours, such as an isolated cell in a CaBasedCellPopulation, store 0.0 for the cell data
cell_iter->GetCellData()->SetItem("mean delta", 0.0);
}
}
}
void MeshBasedCellPopulation<ELEMENT_DIM,SPACE_DIM>::WriteVtkResultsToFile(const std::string& rDirectory)
{
#ifdef CHASTE_VTK
// Store the present time as a string
unsigned num_timesteps = SimulationTime::Instance()->GetTimeStepsElapsed();
std::stringstream time;
time << num_timesteps;
// Store the number of cells for which to output data to VTK
unsigned num_cells_from_mesh = GetNumNodes();
if (!mWriteVtkAsPoints && (mpVoronoiTessellation != NULL))
{
num_cells_from_mesh = mpVoronoiTessellation->GetNumElements();
}
// When outputting any CellData, we assume that the first cell is representative of all cells
unsigned num_cell_data_items = this->Begin()->GetCellData()->GetNumItems();
std::vector<std::string> cell_data_names = this->Begin()->GetCellData()->GetKeys();
std::vector<std::vector<double> > cell_data;
for (unsigned var=0; var<num_cell_data_items; var++)
{
std::vector<double> cell_data_var(num_cells_from_mesh);
cell_data.push_back(cell_data_var);
}
if (mOutputMeshInVtk)
{
// Create mesh writer for VTK output
VtkMeshWriter<ELEMENT_DIM, SPACE_DIM> mesh_writer(rDirectory, "mesh_"+time.str(), false);
mesh_writer.WriteFilesUsingMesh(rGetMesh());
}
if (mWriteVtkAsPoints)
{
// Create mesh writer for VTK output
VtkMeshWriter<SPACE_DIM, SPACE_DIM> cells_writer(rDirectory, "results_"+time.str(), false);
// Iterate over any cell writers that are present
unsigned num_cells = this->GetNumAllCells();
for (typename std::vector<boost::shared_ptr<AbstractCellWriter<ELEMENT_DIM, SPACE_DIM> > >::iterator cell_writer_iter = this->mCellWriters.begin();
cell_writer_iter != this->mCellWriters.end();
++cell_writer_iter)
{
// Create vector to store VTK cell data
std::vector<double> vtk_cell_data(num_cells);
// Loop over cells
for (typename AbstractCellPopulation<ELEMENT_DIM,SPACE_DIM>::Iterator cell_iter = this->Begin();
cell_iter != this->End();
++cell_iter)
{
// Get the node index corresponding to this cell
unsigned node_index = this->GetLocationIndexUsingCell(*cell_iter);
// Populate the vector of VTK cell data
vtk_cell_data[node_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(*cell_iter, this);
}
cells_writer.AddPointData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
}
// Loop over cells
for (typename AbstractCellPopulation<ELEMENT_DIM,SPACE_DIM>::Iterator cell_iter = this->Begin();
cell_iter != this->End();
++cell_iter)
{
// Get the node index corresponding to this cell
unsigned node_index = this->GetLocationIndexUsingCell(*cell_iter);
for (unsigned var=0; var<num_cell_data_items; var++)
{
cell_data[var][node_index] = cell_iter->GetCellData()->GetItem(cell_data_names[var]);
}
}
for (unsigned var=0; var<num_cell_data_items; var++)
{
cells_writer.AddPointData(cell_data_names[var], cell_data[var]);
}
// Make a copy of the nodes in a disposable mesh for writing
{
std::vector<Node<SPACE_DIM>* > nodes;
for (unsigned index=0; index<this->mrMesh.GetNumNodes(); index++)
{
Node<SPACE_DIM>* p_node = this->mrMesh.GetNode(index);
nodes.push_back(p_node);
}
NodesOnlyMesh<SPACE_DIM> mesh;
mesh.ConstructNodesWithoutMesh(nodes, 1.5); // Arbitrary cut off as connectivity not used.
cells_writer.WriteFilesUsingMesh(mesh);
}
*(this->mpVtkMetaFile) << " <DataSet timestep=\"";
*(this->mpVtkMetaFile) << num_timesteps;
*(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
*(this->mpVtkMetaFile) << num_timesteps;
*(this->mpVtkMetaFile) << ".vtu\"/>\n";
}
//.........这里部分代码省略.........
void CombinedStaticGonadBoundaryCondition<DIM>::ImposeBoundaryCondition(const std::map<Node<DIM>*, c_vector<double, DIM> >& rOldLocations)
{
// Iterate over the cell population and get cell location
for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mpCellPopulation->Begin();
cell_iter != this->mpCellPopulation->End();
++cell_iter)
{
c_vector<double,DIM> cell_location = this->mpCellPopulation->GetLocationOfCellCentre(*cell_iter);
Node<DIM>* cell_centre_node = this->mpCellPopulation->GetNode(this->mpCellPopulation->GetLocationIndexUsingCell(*cell_iter));
double radius = cell_centre_node->GetRadius();
// Find C, the closest point on the growth path for this cell, and R the distance to it.
double R=0;
c_vector<double, DIM> C= zero_vector<double>(3);
//If all three parts of the path exist, measure distances to the lower+upper straights and the loop
//Find min distance to each of the three parts of the path
double R1=0; double R2=0; double R3=0;
c_vector<double, DIM> C1=zero_vector<double>(3);
c_vector<double, DIM> C2=zero_vector<double>(3);
c_vector<double, DIM> C3=zero_vector<double>(3);
GetClosestOnLowerStraight(C1,cell_location,R1);
GetClosestOnLoop(C2,cell_location,R2);
GetClosestOnUpperStraight(C3,cell_location,R3);
//Take min and record closest point
R=std::min(R1,std::min(R2,R3));
if(R==R1){
C=C1;
}else if(R==R2){
C=C2;
}else{
C=C3;
}
double distance;
HowFarAlongAreYou(C,distance);
double SyncytiumRadius = GetSyncytiumRadius(distance);
if(distance>mStraightLengthLower){
//Prevent a cell moving back down the tube to somewhere it doesn't fit
if(radius>(mCurrentTubeRadius-SyncytiumRadius)/2){
unsigned node_index = this->mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
Node<DIM>* p_node = this->mpCellPopulation->GetNode(node_index);
c_vector<double,DIM>& p_force = p_node->rGetAppliedForce();
double damping_const = dynamic_cast<NodeBasedCellPopulation<DIM>*>(this->mpCellPopulation)->GetDampingConstant(node_index);
p_node->rGetModifiableLocation() = cell_location-SimulationTime::Instance()->GetTimeStep()*(p_force)/damping_const;
}else{
// If the cell is too far inside the growth path, and therefore in the syncytium...
if (R-radius<SyncytiumRadius-mMaximumDistance)
{
// ...move the cell back onto the surface of the tube by translating away from C:
unsigned node_index = this->mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
Node<DIM>* p_node = this->mpCellPopulation->GetNode(node_index);
p_node->rGetModifiableLocation() = C+(SyncytiumRadius+radius)*(cell_location-C)/norm_2(cell_location-C);
}
// If the cell is too far from the growth path, and therefore outside the tube...
if (R+radius-mCurrentTubeRadius > mMaximumDistance)
{
// ...move the cell back onto the surface of the tube by translating toward C:
unsigned node_index = this->mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
Node<DIM>* p_node = this->mpCellPopulation->GetNode(node_index);
p_node->rGetModifiableLocation() = C+(mCurrentTubeRadius-radius)*(cell_location-C)/norm_2(cell_location-C);
}
}
}else{
// If the cell is too far from the growth path, and therefore outside the tube...
if (R+radius-mCurrentTubeRadius > mMaximumDistance)
{
// ...move the cell back onto the surface of the tube by translating toward C:
unsigned node_index = this->mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
Node<DIM>* p_node = this->mpCellPopulation->GetNode(node_index);
p_node->rGetModifiableLocation() = C+(mCurrentTubeRadius-radius)*(cell_location-C)/norm_2(cell_location-C);
}
}
/*Assuming all is now well, update the cell data to record how far along the gonad arm this cell is.
* Use the vector C, which stores the closest point on the growth path to this cell.*/
HowFarAlongAreYou(C,distance);
cell_iter->GetCellData()->SetItem("DistanceAwayFromDTC", mStraightLengthLower+mStraightLengthUpper+mTurnRadius*M_PI-distance);
if(distance<(mStraightLengthLower-1)){
cell_iter->GetCellData()->SetItem("MaxRadius",mCurrentTubeRadius);
}else{
cell_iter->GetCellData()->SetItem("MaxRadius",(mCurrentTubeRadius-SyncytiumRadius-0.1)/2);
}
}
}
void NodeBasedCellPopulationWithParticles<DIM>::WriteVtkResultsToFile(const std::string& rDirectory)
{
#ifdef CHASTE_VTK
// Store the present time as a string
std::stringstream time;
time << SimulationTime::Instance()->GetTimeStepsElapsed();
// Make sure the nodes are ordered contiguously in memory
NodeMap map(1 + this->mpNodesOnlyMesh->GetMaximumNodeIndex());
this->mpNodesOnlyMesh->ReMesh(map);
// Store the number of cells for which to output data to VTK
unsigned num_nodes = this->GetNumNodes();
std::vector<double> rank(num_nodes);
std::vector<double> particles(num_nodes);
unsigned num_cell_data_items = 0;
std::vector<std::string> cell_data_names;
// We assume that the first cell is representative of all cells
if (num_nodes > 0)
{
num_cell_data_items = this->Begin()->GetCellData()->GetNumItems();
cell_data_names = this->Begin()->GetCellData()->GetKeys();
}
std::vector<std::vector<double> > cell_data;
for (unsigned var=0; var<num_cell_data_items; var++)
{
std::vector<double> cell_data_var(num_nodes);
cell_data.push_back(cell_data_var);
}
// Create mesh writer for VTK output
VtkMeshWriter<DIM, DIM> mesh_writer(rDirectory, "results_"+time.str(), false);
mesh_writer.SetParallelFiles(*(this->mpNodesOnlyMesh));
// Iterate over any cell writers that are present
for (typename std::vector<boost::shared_ptr<AbstractCellWriter<DIM, DIM> > >::iterator cell_writer_iter = this->mCellWriters.begin();
cell_writer_iter != this->mCellWriters.end();
++cell_writer_iter)
{
// Create vector to store VTK cell data
std::vector<double> vtk_cell_data(num_nodes);
// Loop over nodes
for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
node_iter != this->mrMesh.GetNodeIteratorEnd();
++node_iter)
{
unsigned node_index = node_iter->GetIndex();
// If this node is a particle (not a cell), then we set the 'dummy' VTK cell data for this to be -2.0...
if (this->IsParticle(node_index))
{
vtk_cell_data[node_index] = -2.0;
}
else
{
// ...otherwise we populate the vector of VTK cell data as usual
CellPtr p_cell = this->GetCellUsingLocationIndex(node_index);
vtk_cell_data[node_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(p_cell, this);
}
}
mesh_writer.AddPointData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
}
// Loop over cells
for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
cell_iter != this->End();
++cell_iter)
{
// Get the node index corresponding to this cell
unsigned global_index = this->GetLocationIndexUsingCell(*cell_iter);
unsigned node_index = this->rGetMesh().SolveNodeMapping(global_index);
for (unsigned var=0; var<num_cell_data_items; var++)
{
cell_data[var][node_index] = cell_iter->GetCellData()->GetItem(cell_data_names[var]);
}
rank[node_index] = (PetscTools::GetMyRank());
}
mesh_writer.AddPointData("Process rank", rank);
// Loop over nodes
for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
node_iter != this->mrMesh.GetNodeIteratorEnd();
++node_iter)
{
unsigned node_index = node_iter->GetIndex();
particles[node_index] = (double) (this->IsParticle(node_index));
}
mesh_writer.AddPointData("Non-particles", particles);
if (num_cell_data_items > 0)
{
//.........这里部分代码省略.........
unsigned AbstractCellBasedSimulation<ELEMENT_DIM,SPACE_DIM>::DoCellBirth()
{
if (mNoBirth)
{
return 0;
}
unsigned num_births_this_step = 0;
// Iterate over all cells, seeing if each one can be divided
for (typename AbstractCellPopulation<ELEMENT_DIM,SPACE_DIM>::Iterator cell_iter = mrCellPopulation.Begin();
cell_iter != mrCellPopulation.End();
++cell_iter)
{
// Check if this cell is ready to divide
double cell_age = cell_iter->GetAge();
if (cell_age > 0.0)
{
if (cell_iter->ReadyToDivide())
{
// Check if there is room into which the cell may divide
if (mrCellPopulation.IsRoomToDivide(*cell_iter))
{
// Create a new cell
CellPtr p_new_cell = cell_iter->Divide();
// Call method that determines how cell division occurs and returns a vector
c_vector<double, SPACE_DIM> new_location = CalculateCellDivisionVector(*cell_iter);
// If required, output this location to file
/**
* \todo (#2441)
*
* For consistency with the rest of the output code, consider removing the
* AbstractCellBasedSimulation member mOutputDivisionLocations, adding a new
* member mAgesAndLocationsOfDividingCells to AbstractCellPopulation, adding
* a new class CellDivisionLocationsWriter to the CellPopulationWriter hierarchy
* to output the content of mAgesAndLocationsOfDividingCells to file (remembering
* to clear mAgesAndLocationsOfDividingCells at each timestep), and replacing the
* following conditional statement with something like
*
* if (mrCellPopulation.HasWriter<CellDivisionLocationsWriter>())
* {
* mCellDivisionLocations.push_back(new_location);
* }
*/
if (mOutputDivisionLocations)
{
*mpDivisionLocationFile << SimulationTime::Instance()->GetTime() << "\t";
for (unsigned i=0; i<SPACE_DIM; i++)
{
*mpDivisionLocationFile << new_location[i] << "\t";
}
*mpDivisionLocationFile << "\t" << cell_age << "\n";
}
// Add new cell to the cell population
mrCellPopulation.AddCell(p_new_cell, new_location, *cell_iter);
// Update counter
num_births_this_step++;
}
}
}
}
return num_births_this_step;
}
void AbstractCellBasedSimulation<ELEMENT_DIM,SPACE_DIM>::Solve()
{
CellBasedEventHandler::BeginEvent(CellBasedEventHandler::EVERYTHING);
CellBasedEventHandler::BeginEvent(CellBasedEventHandler::SETUP);
// Set up the simulation time
SimulationTime* p_simulation_time = SimulationTime::Instance();
double current_time = p_simulation_time->GetTime();
assert(mDt != DOUBLE_UNSET); //Subclass constructors take care of this
if (mEndTime == DOUBLE_UNSET)
{
EXCEPTION("SetEndTime has not yet been called.");
}
/*
* Note that mDt is used here for "ideal time step". If this step doesn't divide the time remaining
* then a *different* time step will be taken by the time-stepper. The real time-step (used in the
* SimulationTime singleton) is currently not available to this class.
*
* \todo Should we over-write the value of mDt, or change this behaviour? (see #2159)
*/
unsigned num_time_steps = (unsigned) ((mEndTime-current_time)/mDt+0.5);
if (current_time > 0) // use the reset function if necessary
{
p_simulation_time->ResetEndTimeAndNumberOfTimeSteps(mEndTime, num_time_steps);
}
else
{
if (p_simulation_time->IsEndTimeAndNumberOfTimeStepsSetUp())
{
EXCEPTION("End time and number of timesteps already setup. You should not use SimulationTime::SetEndTimeAndNumberOfTimeSteps in cell-based tests.");
}
else
{
p_simulation_time->SetEndTimeAndNumberOfTimeSteps(mEndTime, num_time_steps);
}
}
if (mOutputDirectory == "")
{
EXCEPTION("OutputDirectory not set");
}
double time_now = p_simulation_time->GetTime();
std::ostringstream time_string;
time_string << time_now;
std::string results_directory = mOutputDirectory +"/results_from_time_" + time_string.str();
mSimulationOutputDirectory = results_directory;
// Set up simulation
// Create output files for the visualizer
OutputFileHandler output_file_handler(results_directory+"/", true);
mrCellPopulation.OpenWritersFiles(results_directory+"/");
if (mOutputDivisionLocations)
{
mpDivisionLocationFile = output_file_handler.OpenOutputFile("divisions.dat");
}
if (PetscTools::AmMaster())
{
mpVizSetupFile = output_file_handler.OpenOutputFile("results.vizsetup");
}
// If any PDEs have been defined, set up results files to store their solution
if (mpCellBasedPdeHandler != NULL)
{
mpCellBasedPdeHandler->OpenResultsFiles(this->mSimulationOutputDirectory);
if (PetscTools::AmMaster())
{
*this->mpVizSetupFile << "PDE \n";
}
/*
* If any PDEs have been defined, solve them here before updating cells and store
* their solution in results files. This also initializes the relevant CellData.
* NOTE that this works as the PDEs are elliptic.
*/
CellBasedEventHandler::BeginEvent(CellBasedEventHandler::PDE);
mpCellBasedPdeHandler->SolvePdeAndWriteResultsToFile(this->mSamplingTimestepMultiple);
CellBasedEventHandler::EndEvent(CellBasedEventHandler::PDE);
}
SetupSolve();
// Call SetupSolve() on each modifier
for (typename std::vector<boost::shared_ptr<AbstractCellBasedSimulationModifier<ELEMENT_DIM, SPACE_DIM> > >::iterator iter = mSimulationModifiers.begin();
iter != mSimulationModifiers.end();
++iter)
{
(*iter)->SetupSolve(this->mrCellPopulation,this->mSimulationOutputDirectory);
}
/*
* Age the cells to the correct time. Note that cells are created with
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论