C++ typenameAbstractCellPopulation类代码示例

void GonadArmPositionTrackerModifier<DIM>::UpdateAtEndOfTimeStep(AbstractCellPopulation<DIM,DIM>& rCellPopulation)
{

	if((SimulationTime::Instance()->GetTime())-(int)SimulationTime::Instance()->GetTime()==(1/250))

    for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = rCellPopulation.Begin();
         cell_iter != rCellPopulation.End();
         ++cell_iter)
    {

    	double id = cell_iter->GetCellId();
    	double position = cell_iter->GetCellData()->GetItem("DistanceAwayFromDTC");
    	double prolif = cell_iter->GetCellData()->GetItem("Proliferating");

    	fprintf(OutputPositionFile,"%e\t",SimulationTime::Instance()->GetTime());
    	fprintf(OutputPositionFile,"%f\t",id);
    	fprintf(OutputPositionFile,"%e\t",position);
    	fprintf(OutputPositionFile,"%e\t",prolif);
    	fprintf(OutputPositionFile,"\n");
    }

	if(SimulationTime::Instance()->IsFinished()){
		//Close output file
		fclose(OutputPositionFile);
	}
}

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 AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::SetCellAncestorsToLocationIndices()
{
    for (typename AbstractCellPopulation<ELEMENT_DIM, SPACE_DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
    {
        MAKE_PTR_ARGS(CellAncestor, p_cell_ancestor, (mCellLocationMap[(*cell_iter).get()]));
        cell_iter->SetAncestor(p_cell_ancestor);
    }
}

void PlaneBasedCellKiller<DIM>::CheckAndLabelCellsForApoptosisOrDeath()
{
    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);

        if (inner_prod(cell_location - mPointOnPlane, mNormalToPlane) > 0.0)
        {
            cell_iter->Kill();
        }
    }
}

void TimedPlaneBasedCellKiller<DIM>::CheckAndLabelCellsForApoptosisOrDeath()
{
    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);

        if (inner_prod(cell_location - mPointOnPlane, mNormalToPlane) > 0.0 && SimulationTime::Instance()->GetTime()>=mStartTime)
        {
        	 OUTPUT<<SimulationTime::Instance()->GetTime()<<"	"<<(*cell_iter)->GetCellId()<<std::endl;
            cell_iter->Kill();
        }
    }
}

void SloughingCellKiller<DIM>::CheckAndLabelCellsForApoptosisOrDeath()
{
    switch (DIM)
    {
        case 1:
        {
            for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mpCellPopulation->Begin();
                 cell_iter != this->mpCellPopulation->End();
                 ++cell_iter)
            {
                double x = this->mpCellPopulation->GetLocationOfCellCentre(*cell_iter)[0];

                if (x > mSloughHeight)
                {
                    cell_iter->Kill();
                }
            }
            break;
        }
        case 2:
        {
            for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->mpCellPopulation->Begin();
                 cell_iter != this->mpCellPopulation->End();
                 ++cell_iter)
            {
                c_vector<double, 2> location;
                location = this->mpCellPopulation->GetLocationOfCellCentre(*cell_iter);
                double x = location[0];
                double y = location[1];

                if ((y>mSloughHeight) || (mSloughSides && ((x<0.0) || (x>mSloughWidth))))
                {
                    cell_iter->Kill();
                }
            }
            break;
        }
        case 3:
        {
            EXCEPTION("SloughingCellKiller is not yet implemented in 3D");
            break;
        }
        default:
            // This can't happen
            NEVER_REACHED;
    }
}

void NodeBasedCellPopulation<DIM>::Update(bool hasHadBirthsOrDeaths)
{
    UpdateCellProcessLocation();

    mpNodesOnlyMesh->UpdateBoxCollection();

    if (mLoadBalanceMesh)
    {
        if ((SimulationTime::Instance()->GetTimeStepsElapsed() % mLoadBalanceFrequency) == 0)
        {
            mpNodesOnlyMesh->LoadBalanceMesh();

            UpdateCellProcessLocation();

            mpNodesOnlyMesh->UpdateBoxCollection();
        }
    }

    RefreshHaloCells();

    mpNodesOnlyMesh->CalculateInteriorNodePairs(mNodePairs, mNodeNeighbours);

    AddReceivedHaloCells();

    mpNodesOnlyMesh->CalculateBoundaryNodePairs(mNodePairs, mNodeNeighbours);

    /*
     * Update cell radii based on CellData
     */
    if (mUseVariableRadii)
    {
        for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
             cell_iter != this->End();
             ++cell_iter)
        {
            double cell_radius = cell_iter->GetCellData()->GetItem("Radius");
            unsigned node_index = this->GetLocationIndexUsingCell(*cell_iter);
            this->GetNode(node_index)->SetRadius(cell_radius);
        }
    }

    // Make sure that everyone exits update together so that all asynchronous communications are complete.
    PetscTools::Barrier("Update");
}

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 CellBasedPdeHandler<DIM>::WritePdeSolution(double time)
{
    if (PetscTools::AmMaster())
    {
        (*mpVizPdeSolutionResultsFile) << time << "\t";

#ifdef CHASTE_VTK
        // Note that this mesh writer is only constructed and used if mpCoarsePdeMesh exists
        VtkMeshWriter<DIM,DIM>* p_vtk_mesh_writer = NULL;
        if (DIM>1 && mpCoarsePdeMesh != NULL )
        {
            std::ostringstream time_string;
            time_string << SimulationTime::Instance()->GetTimeStepsElapsed()+1;
            std::string results_file = "pde_results_"+time_string.str();
            // Note that this mesh writer is always constructed, but is only used if mpCoarsePdeMesh exists
            p_vtk_mesh_writer = new VtkMeshWriter<DIM,DIM>(mDirPath, results_file, false);
        }
#endif //CHASTE_VTK
        for (unsigned pde_index=0; pde_index<mPdeAndBcCollection.size(); pde_index++)
        {
            if (mpCoarsePdeMesh != NULL)
            {
                PdeAndBoundaryConditions<DIM>* p_pde_and_bc = mPdeAndBcCollection[pde_index];
                assert( p_pde_and_bc->rGetDependentVariableName()!="");

#ifdef CHASTE_VTK
                if (p_pde_and_bc->GetSolution())
                {
                    if (DIM>1)
                    {
                        ReplicatableVector solution_repl(p_pde_and_bc->GetSolution());
                        std::vector<double> pde_solution;
                        for (unsigned i=0; i<mpCoarsePdeMesh->GetNumNodes(); i++)
                        {
                            pde_solution.push_back(solution_repl[i]);
                        }

                        p_vtk_mesh_writer->AddPointData(p_pde_and_bc->rGetDependentVariableName(),pde_solution);
                    }
                }

#endif //CHASTE_VTK

                for (unsigned i=0; i<mpCoarsePdeMesh->GetNumNodes(); i++)
                {
                    (*mpVizPdeSolutionResultsFile) << i << " ";
                    c_vector<double,DIM> location = mpCoarsePdeMesh->GetNode(i)->rGetLocation();
                    for (unsigned k=0; k<DIM; k++)
                    {
                        (*mpVizPdeSolutionResultsFile) << location[k] << " ";
                    }

                    if (p_pde_and_bc->GetSolution())
                    {
                        ReplicatableVector solution_repl(p_pde_and_bc->GetSolution());
                        (*mpVizPdeSolutionResultsFile) << solution_repl[i] << " ";
                    }
                    else
                    {
                        // should only come into this method AFTER solving the PDE
                        NEVER_REACHED;
                    }
                }
            }
            else // Not coarse mesh
            {
                for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = mpCellPopulation->Begin();
                     cell_iter != mpCellPopulation->End();
                     ++cell_iter)
                {
                    unsigned node_index = mpCellPopulation->GetLocationIndexUsingCell(*cell_iter);
                    (*mpVizPdeSolutionResultsFile) << node_index << " ";
                    const c_vector<double,DIM>& position = mpCellPopulation->GetLocationOfCellCentre(*cell_iter);
                    for (unsigned i=0; i<DIM; i++)
                    {
                        (*mpVizPdeSolutionResultsFile) << position[i] << " ";
                    }
                    double solution = cell_iter->GetCellData()->GetItem(mPdeAndBcCollection[pde_index]->rGetDependentVariableName());
                    (*mpVizPdeSolutionResultsFile) << solution << " ";
                }
            }
        }
        (*mpVizPdeSolutionResultsFile) << "\n";
#ifdef CHASTE_VTK
        if (p_vtk_mesh_writer != NULL)
        {
            p_vtk_mesh_writer->WriteFilesUsingMesh(*mpCoarsePdeMesh);
            delete p_vtk_mesh_writer;
        }
#endif //CHASTE_VTK
    }
}

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)
    {
//.........这里部分代码省略.........

void NodeBasedCellPopulation<DIM>::WriteVtkResultsToFile(const std::string& rDirectory)
{
#ifdef CHASTE_VTK
    std::stringstream time;
    time << SimulationTime::Instance()->GetTimeStepsElapsed();
    VtkMeshWriter<DIM, DIM> mesh_writer(rDirectory, "results_"+time.str(), false);

    // Make sure the nodes are ordered contiguously in memory.
    NodeMap map(1 + this->mpNodesOnlyMesh->GetMaximumNodeIndex());
    this->mpNodesOnlyMesh->ReMesh(map);

    mesh_writer.SetParallelFiles(*mpNodesOnlyMesh);

    unsigned num_nodes = GetNumNodes();
    std::vector<double> cell_types(num_nodes);
    std::vector<double> cell_ancestors(num_nodes);
    std::vector<double> cell_mutation_states(num_nodes);
    std::vector<double> cell_ages(num_nodes);
    std::vector<double> cell_cycle_phases(num_nodes);
    std::vector<double> cell_radii(num_nodes);
    std::vector<std::vector<double> > cellwise_data;
    std::vector<double> rank(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();
    }

    for (unsigned var=0; var<num_cell_data_items; var++)
    {
        std::vector<double> cellwise_data_var(num_nodes);
        cellwise_data.push_back(cellwise_data_var);
    }

    // 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);

        Node<DIM>* p_node = this->GetNode(global_index);

        unsigned node_index = this->rGetMesh().SolveNodeMapping(global_index);

        if (this-> template HasWriter<CellAncestorWriter>())
        {
            double ancestor_index = (cell_iter->GetAncestor() == UNSIGNED_UNSET) ? (-1.0) : (double)cell_iter->GetAncestor();
            cell_ancestors[node_index] = ancestor_index;
        }
        if (this-> template HasWriter<CellProliferativeTypesWriter>())
        {
            double cell_type = cell_iter->GetCellProliferativeType()->GetColour();
            cell_types[node_index] = cell_type;
        }
        if (this-> template HasWriter<CellMutationStatesCountWriter>())
        {
            double mutation_state = cell_iter->GetMutationState()->GetColour();

            CellPropertyCollection collection = cell_iter->rGetCellPropertyCollection();
            CellPropertyCollection label_collection = collection.GetProperties<CellLabel>();

            if (label_collection.GetSize() == 1)
            {
                boost::shared_ptr<CellLabel> p_label = boost::static_pointer_cast<CellLabel>(label_collection.GetProperty());
                mutation_state = p_label->GetColour();
            }

            cell_mutation_states[node_index] = mutation_state;
        }
        if (this-> template HasWriter<CellAgesWriter>())
        {
            double age = cell_iter->GetAge();
            cell_ages[node_index] = age;
        }
        if (this-> template HasWriter<CellProliferativePhasesWriter>())
        {
            double cycle_phase = cell_iter->GetCellCycleModel()->GetCurrentCellCyclePhase();
            cell_cycle_phases[node_index] = cycle_phase;
        }
        if (this-> template HasWriter<CellVolumesWriter>())
        {
            double cell_radius = p_node->GetRadius();
            cell_radii[node_index] = cell_radius;
        }

        for (unsigned var=0; var<num_cell_data_items; var++)
        {
            cellwise_data[var][node_index] = cell_iter->GetCellData()->GetItem(cell_data_names[var]);
        }

        rank[node_index] = (PetscTools::GetMyRank());
    }

//.........这里部分代码省略.........

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;
}