void NodeBasedCellPopulationWithParticles<DIM>::UpdateParticlePositions(double dt)
{
    // Initialise vector of forces on particles
    std::vector<c_vector<double, DIM> > drdt(this->GetNumNodes());
    for (unsigned i=0; i<drdt.size(); i++)
    {
        drdt[i] = zero_vector<double>(DIM);
    }

    // Calculate forces on particles
    double damping_constant = this->GetDampingConstantNormal();
    for (unsigned i=0; i<drdt.size(); i++)
    {
        drdt[i] = this->GetNode(i)->rGetAppliedForce()/damping_constant;
    }

    for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
         node_iter != this->mrMesh.GetNodeIteratorEnd();
         ++node_iter)
    {
        if (node_iter->IsParticle())
        {
            ChastePoint<DIM> new_point(node_iter->rGetLocation() + dt*drdt[node_iter->GetIndex()]);
            node_iter->SetPoint(new_point);
        }
    }
}

void VertexCryptBoundaryForce<DIM>::AddForceContribution(AbstractCellPopulation<DIM>& rCellPopulation)
{
    // Helper variable that is a static cast of the cell population
    VertexBasedCellPopulation<DIM>* p_cell_population = static_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation);

    // Throw an exception message if not using a VertexBasedCellPopulation
    if (dynamic_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation) == nullptr)
    {
        EXCEPTION("VertexCryptBoundaryForce is to be used with VertexBasedCellPopulations only");
    }

    // Iterate over nodes
    for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = p_cell_population->rGetMesh().GetNodeIteratorBegin();
         node_iter != p_cell_population->rGetMesh().GetNodeIteratorEnd();
         ++node_iter)
    {
        double y = node_iter->rGetLocation()[1]; // y-coordinate of node

        // If the node lies below the line y=0, then add the boundary force contribution to the node forces
        if (y < 0.0)
        {
            c_vector<double, DIM> boundary_force = zero_vector<double>(DIM);
            boundary_force[1] = mForceStrength*SmallPow(y, 2);

            node_iter->AddAppliedForceContribution(boundary_force);
        }
    }
}

void NodeBasedCellPopulationWithParticles<DIM>::Validate()
{
    std::map<unsigned, bool> validated_nodes;
    for (typename AbstractMesh<DIM, DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
         node_iter != this->mrMesh.GetNodeIteratorEnd();
         ++node_iter)
    {
        validated_nodes[node_iter->GetIndex()] = node_iter->IsParticle();
    }

    // Look through all of the cells and record what node they are associated with.
    for (typename AbstractCellPopulation<DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
    {
        unsigned node_index = this->GetLocationIndexUsingCell((*cell_iter));

        // If the node attached to this cell is labelled as a particle, then throw an error
        if (this->GetNode(node_index)->IsParticle())
        {
            EXCEPTION("Node " << node_index << " is labelled as a particle and has a cell attached");
        }
        validated_nodes[node_index] = true;
    }

    for (std::map<unsigned, bool>::iterator map_iter = validated_nodes.begin();
         map_iter != validated_nodes.end();
         map_iter++)
    {
        if (!map_iter->second)
        {
            EXCEPTION("Node " << map_iter->first << " does not appear to be a particle or has a cell associated with it");
        }
    }
}

void ForwardEulerNumericalMethod<ELEMENT_DIM,SPACE_DIM>::UpdateAllNodePositions(double dt)
{
    if (!this->mUseUpdateNodeLocation)
    {
        // Apply forces to each cell, and save a vector of net forces F
        std::vector<c_vector<double, SPACE_DIM> > forces = this->ComputeForcesIncludingDamping();

        unsigned index = 0;
        for (typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::NodeIterator node_iter = this->mpCellPopulation->rGetMesh().GetNodeIteratorBegin();
             node_iter != this->mpCellPopulation->rGetMesh().GetNodeIteratorEnd();
             ++node_iter, ++index)
        {
            // Get the current node location and calculate the new location according to the forward Euler method
            const c_vector<double, SPACE_DIM>& r_old_location = node_iter->rGetLocation();
            c_vector<double, SPACE_DIM> displacement = dt * forces[index];

            // In the vertex-based case, the displacement may be scaled if the cell rearrangement threshold is exceeded
            this->DetectStepSizeExceptions(node_iter->GetIndex(), displacement, dt);

            c_vector<double, SPACE_DIM> new_location = r_old_location + displacement;
            this->SafeNodePositionUpdate(node_iter->GetIndex(), new_location);
        }
    }
    else
    {
        /*
         * If this type of cell population does not support the new numerical methods, delegate
         * updating node positions to the population itself.
         *
         * This only applies to NodeBasedCellPopulationWithBuskeUpdates.
         */
        this->mpCellPopulation->UpdateNodeLocations(dt);
    }
}

void OffLatticeSimulation<ELEMENT_DIM,SPACE_DIM>::UpdateCellLocationsAndTopology()
{
    // Calculate forces
    CellBasedEventHandler::BeginEvent(CellBasedEventHandler::FORCE);

    // Clear all forces
    for (typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::NodeIterator node_iter = this->mrCellPopulation.rGetMesh().GetNodeIteratorBegin();
         node_iter != this->mrCellPopulation.rGetMesh().GetNodeIteratorEnd();
         ++node_iter)
    {
        node_iter->ClearAppliedForce();
    }

    // Now add force contributions from each AbstractForce
    for (typename std::vector<boost::shared_ptr<AbstractForce<ELEMENT_DIM, SPACE_DIM> > >::iterator iter = mForceCollection.begin();
         iter != mForceCollection.end();
         ++iter)
    {
        (*iter)->AddForceContribution(this->mrCellPopulation);
    }
    CellBasedEventHandler::EndEvent(CellBasedEventHandler::FORCE);

    // Update node positions
    CellBasedEventHandler::BeginEvent(CellBasedEventHandler::POSITION);
    UpdateNodePositions();
    CellBasedEventHandler::EndEvent(CellBasedEventHandler::POSITION);
}

MeshBasedCellPopulation<ELEMENT_DIM,SPACE_DIM>::MeshBasedCellPopulation(MutableMesh<ELEMENT_DIM,SPACE_DIM>& rMesh,
                                      std::vector<CellPtr>& rCells,
                                      const std::vector<unsigned> locationIndices,
                                      bool deleteMesh,
                                      bool validate)
    : AbstractCentreBasedCellPopulation<ELEMENT_DIM,SPACE_DIM>(rMesh, rCells, locationIndices),
      mpVoronoiTessellation(NULL),
      mDeleteMesh(deleteMesh),
      mUseAreaBasedDampingConstant(false),
      mAreaBasedDampingConstantParameter(0.1),
      mWriteVtkAsPoints(false),
      mOutputMeshInVtk(false),
      mHasVariableRestLength(false)
{
    mpMutableMesh = static_cast<MutableMesh<ELEMENT_DIM,SPACE_DIM>* >(&(this->mrMesh));

    assert(this->mCells.size() <= this->mrMesh.GetNumNodes());

    if (validate)
    {
        Validate();
    }

    // Initialise the applied force at each node to zero
    for (typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::NodeIterator node_iter = this->rGetMesh().GetNodeIteratorBegin();
         node_iter != this->rGetMesh().GetNodeIteratorEnd();
         ++node_iter)
    {
        node_iter->ClearAppliedForce();
    }
}

void TCellDiffusionForce<DIM>::AddForceContribution(AbstractCellPopulation<DIM>& rCellPopulation)
{
    double dt = SimulationTime::Instance()->GetTimeStep();

    // Iterate over the nodes
    for (typename AbstractMesh<DIM, DIM>::NodeIterator node_iter = rCellPopulation.rGetMesh().GetNodeIteratorBegin();
         node_iter != rCellPopulation.rGetMesh().GetNodeIteratorEnd();
         ++node_iter)
    {
        // Get the radius of this node
        unsigned node_index = node_iter->GetIndex();
        double node_radius = node_iter->GetRadius();
        
        // Get cell associated with this index
        CellPtr p_cell = rCellPopulation.GetCellUsingLocationIndex(node_index);

        // Reject if no radius has been set
        if (node_radius == 0.0)
        {
            EXCEPTION("SetRadius() must be called on each Node before calling TCellDiffusionForce::AddForceContribution() to avoid a division by zero error");
        }
        
        //   If the selected cell is a Unlabelled Differentiated T Cell, apply diffusion force contribution.
        if (  (p_cell->GetMutationState()->IsType<TCellMutationState>()) && (p_cell->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>()) 
            && !(p_cell->HasCellProperty<CellLabel>())  )
        {
            double nu = dynamic_cast<AbstractOffLatticeCellPopulation<DIM>*>(&rCellPopulation)->GetDampingConstant(node_index);

            /* Compute the diffusion coefficient D as D = k*T/(6*pi*eta*r), where
             *
             * k = Boltzmann's constant,
             * T = absolute temperature,
             * eta = dynamic viscosity,
             * r = cell radius. */
            double diffusion_const_scaling = GetDiffusionScalingConstant();
            double diffusion_constant = diffusion_const_scaling/node_radius;

            c_vector<double, DIM> force_contribution;
            for (unsigned i=0; i<DIM; i++)
            {
                /* The force on this cell is scaled with the timestep such that when it is
                 * used in the discretised equation of motion for the cell, we obtain the
                 * correct formula
                 *
                 * x_new = x_old + sqrt(2*D*dt)*W
                 *
                 * where W is a standard normal random variable. */
                 
                double xi = RandomNumberGenerator::Instance()->StandardNormalRandomDeviate();
                force_contribution[i] = mStrengthParameter * ((nu*sqrt(2.0*diffusion_constant*dt)/dt)*xi);
            }
            node_iter->AddAppliedForceContribution(force_contribution);
        }
    }
    
    
}

void OffLatticeSimulation<ELEMENT_DIM,SPACE_DIM>::SetupSolve()
{
    // Clear all forces
    for (typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::NodeIterator node_iter = this->mrCellPopulation.rGetMesh().GetNodeIteratorBegin();
         node_iter != this->mrCellPopulation.rGetMesh().GetNodeIteratorEnd();
         ++node_iter)
    {
        node_iter->ClearAppliedForce();
    }
}

void NodeBasedCellPopulationWithParticles<DIM>::SetParticles(const std::set<unsigned>& rParticleIndices)
{
    for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
         node_iter != this->mrMesh.GetNodeIteratorEnd();
         ++node_iter)
    {
        if (rParticleIndices.find(node_iter->GetIndex()) != rParticleIndices.end())
        {
            node_iter->SetIsParticle(true);
        }
    }
    NodeBasedCellPopulationWithParticles::Validate();
}

void NodeBasedCellPopulation<DIM>::Validate()
{
    for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
         node_iter != this->mrMesh.GetNodeIteratorEnd();
         ++node_iter)
    {
        try
        {
            this->GetCellUsingLocationIndex(node_iter->GetIndex());
        }
        catch (Exception&)
        {
            EXCEPTION("Node " << node_iter->GetIndex() << " does not appear to have a cell associated with it");
        }
    }
}

void MeshBasedCellPopulationWithGhostNodes<DIM>::ApplyGhostForces(){

    // Initialise vector of forces on ghost nodes
    std::vector<c_vector<double, DIM> > drdt(this->GetNumNodes());
    for (unsigned i=0; i<drdt.size(); i++)
    {
        drdt[i] = zero_vector<double>(DIM);
    }

    // Calculate forces on ghost nodes
    for (typename MutableMesh<DIM, DIM>::EdgeIterator edge_iterator = static_cast<MutableMesh<DIM, DIM>&>((this->mrMesh)).EdgesBegin();
        edge_iterator != static_cast<MutableMesh<DIM, DIM>&>((this->mrMesh)).EdgesEnd();
        ++edge_iterator)
    {
        unsigned nodeA_global_index = edge_iterator.GetNodeA()->GetIndex();
        unsigned nodeB_global_index = edge_iterator.GetNodeB()->GetIndex();

        c_vector<double, DIM> force = CalculateForceBetweenGhostNodes(nodeA_global_index, nodeB_global_index);

        if (!this->mIsGhostNode[nodeA_global_index])
        {
            drdt[nodeB_global_index] -= force;
        }
        else
        {
            drdt[nodeA_global_index] += force;

            if (this->mIsGhostNode[nodeB_global_index])
            {
                drdt[nodeB_global_index] -= force;
            }
        }
    }

    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->mIsGhostNode[node_index])
        {
            node_iter->ClearAppliedForce();
            node_iter->AddAppliedForceContribution(drdt[node_index]);
        }
    }

};

void OffLatticeSimulation<ELEMENT_DIM,SPACE_DIM>::SetupSolve()
{
    // Clear all forces
    for (typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::NodeIterator node_iter = this->mrCellPopulation.rGetMesh().GetNodeIteratorBegin();
         node_iter != this->mrCellPopulation.rGetMesh().GetNodeIteratorEnd();
         ++node_iter)
    {
        node_iter->ClearAppliedForce();
    }

    // Use a forward Euler method by default, unless a numerical method has been specified already
    if (mpNumericalMethod == nullptr)
    {
        mpNumericalMethod = boost::make_shared<ForwardEulerNumericalMethod<ELEMENT_DIM, SPACE_DIM> >();
    }
    mpNumericalMethod->SetCellPopulation(dynamic_cast<AbstractOffLatticeCellPopulation<ELEMENT_DIM,SPACE_DIM>*>(&(this->mrCellPopulation)));
    mpNumericalMethod->SetForceCollection(&mForceCollection);
}

void MeshBasedCellPopulationWithGhostNodes<DIM>::AcceptCellWritersAcrossPopulation()
{
    for (typename AbstractMesh<DIM, DIM>::NodeIterator node_iter = this->rGetMesh().GetNodeIteratorBegin();
         node_iter != this->rGetMesh().GetNodeIteratorEnd();
         ++node_iter)
    {
        // If it isn't a ghost node then there might be cell writers attached
        if (! this->IsGhostNode(node_iter->GetIndex()))
        {
            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)
            {
                CellPtr cell_from_node = this->GetCellUsingLocationIndex(node_iter->GetIndex());
                this->AcceptCellWriter(*cell_writer_iter, cell_from_node);
            }
        }
    }
}

NodeBasedCellPopulationWithParticles<DIM>::NodeBasedCellPopulationWithParticles(NodesOnlyMesh<DIM>& rMesh,
                                      std::vector<CellPtr>& rCells,
                                      const std::vector<unsigned> locationIndices,
                                      bool deleteMesh)
    : NodeBasedCellPopulation<DIM>(rMesh, rCells, locationIndices, deleteMesh, false)
{
    EXCEPT_IF_NOT(PetscTools::IsSequential());
    if (!locationIndices.empty())
    {
        // Create a set of node indices corresponding to particles
        std::set<unsigned> node_indices;
        std::set<unsigned> location_indices;
        std::set<unsigned> particle_indices;

        for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = rMesh.GetNodeIteratorBegin();
             node_iter != rMesh.GetNodeIteratorEnd();
             ++node_iter)
        {
            node_indices.insert(node_iter->GetIndex());
        }
        for (unsigned i=0; i<locationIndices.size(); i++)
        {
            location_indices.insert(locationIndices[i]);
        }

        std::set_difference(node_indices.begin(), node_indices.end(),
                            location_indices.begin(), location_indices.end(),
                            std::inserter(particle_indices, particle_indices.begin()));

        // This method finishes and then calls Validate()
        SetParticles(particle_indices);
    }
    else
    {
        for (typename NodesOnlyMesh<DIM>::NodeIterator node_iter = rMesh.GetNodeIteratorBegin();
             node_iter != rMesh.GetNodeIteratorEnd();
             ++node_iter)
        {
            (*node_iter).SetIsParticle(false);
        }
        NodeBasedCellPopulationWithParticles::Validate();
    }
}

void MeshBasedCellPopulation<ELEMENT_DIM,SPACE_DIM>::Update(bool hasHadBirthsOrDeaths)
{
    ///\todo check if there is a more efficient way of keeping track of node velocity information (#2404)
    bool output_node_velocities = (this-> template HasWriter<NodeVelocityWriter>());

    /**
     * If node radii are set, then we must keep a record of these, since they will be cleared during
     * the remeshing process. We then restore these attributes to the nodes after calling ReMesh().
     *
     * At present, we check whether node radii are set by interrogating the radius of the first node
     * in the mesh and asking if it is strictly greater than zero (the default value, as set in the
     * NodeAttributes constructor). Hence, we assume that either ALL node radii are set, or NONE are.
     *
     * \todo There may be a better way of checking if node radii are set (#2694)
     */
    std::map<unsigned, double> old_node_radius_map;
    old_node_radius_map.clear();
    if (this->mrMesh.GetNodeIteratorBegin()->HasNodeAttributes())
    {
        if (this->mrMesh.GetNodeIteratorBegin()->GetRadius() > 0.0)
        {
            for (typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
                 node_iter != this->mrMesh.GetNodeIteratorEnd();
                 ++node_iter)
            {
                unsigned node_index = node_iter->GetIndex();
                old_node_radius_map[node_index] = node_iter->GetRadius();
            }
        }
    }

    std::map<unsigned, c_vector<double, SPACE_DIM> > old_node_applied_force_map;
    old_node_applied_force_map.clear();
    if (output_node_velocities)
    {
        /*
         * If outputting node velocities, we must keep a record of the applied force at each
         * node, since this will be cleared during the remeshing process. We then restore
         * these attributes to the nodes after calling ReMesh().
         */
        for (typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
             node_iter != this->mrMesh.GetNodeIteratorEnd();
             ++node_iter)
        {
            unsigned node_index = node_iter->GetIndex();
            old_node_applied_force_map[node_index] = node_iter->rGetAppliedForce();
        }
    }

    NodeMap node_map(this->mrMesh.GetNumAllNodes());

    // We must use a static_cast to call ReMesh() as this method is not defined in parent mesh classes
    static_cast<MutableMesh<ELEMENT_DIM,SPACE_DIM>&>((this->mrMesh)).ReMesh(node_map);

    if (!node_map.IsIdentityMap())
    {
        UpdateGhostNodesAfterReMesh(node_map);

        // Update the mappings between cells and location indices
        std::map<Cell*, unsigned> old_cell_location_map = this->mCellLocationMap;

        // Remove any dead pointers from the maps (needed to avoid archiving errors)
        this->mLocationCellMap.clear();
        this->mCellLocationMap.clear();

        for (std::list<CellPtr>::iterator it = this->mCells.begin(); it != this->mCells.end(); ++it)
        {
            unsigned old_node_index = old_cell_location_map[(*it).get()];

            // This shouldn't ever happen, as the cell vector only contains living cells
            assert(!node_map.IsDeleted(old_node_index));

            unsigned new_node_index = node_map.GetNewIndex(old_node_index);
            this->SetCellUsingLocationIndex(new_node_index,*it);

            if (old_node_radius_map[old_node_index] > 0.0)
            {
                this->GetNode(new_node_index)->SetRadius(old_node_radius_map[old_node_index]);
            }
            if (output_node_velocities)
            {
                this->GetNode(new_node_index)->AddAppliedForceContribution(old_node_applied_force_map[old_node_index]);
            }
        }

        this->Validate();
    }
    else
    {
        if (old_node_radius_map[this->mCellLocationMap[(*(this->mCells.begin())).get()]] > 0.0)
        {
            for (std::list<CellPtr>::iterator it = this->mCells.begin(); it != this->mCells.end(); ++it)
            {
                unsigned node_index = this->mCellLocationMap[(*it).get()];
                this->GetNode(node_index)->SetRadius(old_node_radius_map[node_index]);
            }
        }
        if (output_node_velocities)
        {
            for (std::list<CellPtr>::iterator it = this->mCells.begin(); it != this->mCells.end(); ++it)
//.........这里部分代码省略.........

void PottsBasedCellPopulation<DIM>::WriteVtkResultsToFile(const std::string& rDirectory)
{
#ifdef CHASTE_VTK
    unsigned num_timesteps = SimulationTime::Instance()->GetTimeStepsElapsed();
    std::stringstream time;
    time << num_timesteps;

    // Create mesh writer for VTK output
    VtkMeshWriter<DIM, DIM> mesh_writer(rDirectory, "results_"+time.str(), false);

    // Iterate over any cell writers that are present
    unsigned num_nodes = GetNumNodes();
    for (typename std::set<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);

        // Iterate over nodes in the mesh
        for (typename AbstractMesh<DIM,DIM>::NodeIterator iter = mpPottsMesh->GetNodeIteratorBegin();
             iter != mpPottsMesh->GetNodeIteratorEnd();
             ++iter)
        {
            // Get the index of this node in the mesh and those elements (i.e. cells) that contain this node
            unsigned node_index = iter->GetIndex();
            std::set<unsigned> element_indices = iter->rGetContainingElementIndices();

            // If there are no elements associated with this node, then we set the value of any VTK cell data to be -1 at this node...
            if (element_indices.empty())
            {
                // Populate the vector of VTK cell data
                vtk_cell_data[node_index] = -1.0;
            }
            else
            {
                // ... otherwise there should be exactly one element (i.e. cell) containing this node
                assert(element_indices.size() == 1);
                unsigned elem_index = *(element_indices.begin());
                CellPtr p_cell = this->GetCellUsingLocationIndex(elem_index);

                // Populate the vector of VTK cell data
                vtk_cell_data[node_index] = (*cell_writer_iter)->GetCellDataForVtkOutput(p_cell, this);
            }
        }

        mesh_writer.AddPointData((*cell_writer_iter)->GetVtkCellDataName(), vtk_cell_data);
    }


    // 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_nodes);
        cell_data.push_back(cell_data_var);
    }

    for (typename AbstractMesh<DIM,DIM>::NodeIterator iter = mpPottsMesh->GetNodeIteratorBegin();
         iter != mpPottsMesh->GetNodeIteratorEnd();
         ++iter)
    {
        // Get the index of this node in the mesh and those elements (i.e. cells) that contain this node
        unsigned node_index = iter->GetIndex();
        std::set<unsigned> element_indices = iter->rGetContainingElementIndices();

        // If there are no elements associated with this node, then we set the value of any VTK cell data to be -1 at this node...
        if (element_indices.empty())
        {
            for (unsigned var=0; var<num_cell_data_items; var++)
            {
                cell_data[var][node_index] = -1.0;
            }
        }
        else
        {
            // ... otherwise there should be exactly one element (i.e. cell) containing this node
            assert(element_indices.size() == 1);
            unsigned elem_index = *(element_indices.begin());
            CellPtr p_cell = this->GetCellUsingLocationIndex(elem_index);

            for (unsigned var=0; var<num_cell_data_items; var++)
            {
                cell_data[var][node_index] = p_cell->GetCellData()->GetItem(cell_data_names[var]);
            }
        }
    }
    for (unsigned var=0; var<cell_data.size(); var++)
    {
        mesh_writer.AddPointData(cell_data_names[var], cell_data[var]);
    }

    /*
     * The current VTK writer can only write things which inherit from AbstractTetrahedralMeshWriter.
     * For now, we do an explicit conversion to NodesOnlyMesh. This can be written to VTK then visualized as glyphs.
     */
    NodesOnlyMesh<DIM> temp_mesh;
//.........这里部分代码省略.........

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 AbstractContinuumMechanicsSolver<DIM>::AllocateMatrixMemory()
{
    Vec template_vec = mrQuadMesh.GetDistributedVectorFactory()->CreateVec(mProblemDimension);

    ///////////////////////////
    // three vectors
    ///////////////////////////
    VecDuplicate(template_vec, &mResidualVector);
    VecDuplicate(mResidualVector, &mLinearSystemRhsVector);
    // the one is only allocated if it will be needed (in ApplyDirichletBoundaryConditions),
    // depending on whether the matrix is kept symmetric.
    mDirichletBoundaryConditionsVector = NULL;
    PetscTools::Destroy(template_vec);

    ///////////////////////////
    // two matrices
    ///////////////////////////

    int lo, hi;
    VecGetOwnershipRange(mResidualVector, &lo, &hi);
    PetscInt local_size = hi - lo;


    if (DIM==2)
    {
        // 2D: N elements around a point => 7N+3 non-zeros in that row? Assume N<=10 (structured mesh would have N_max=6) => 73.
        unsigned num_non_zeros = std::min(75u, mNumDofs);

        PetscTools::SetupMat(mSystemLhsMatrix, mNumDofs, mNumDofs, num_non_zeros, local_size, local_size);
        PetscTools::SetupMat(mPreconditionMatrix, mNumDofs, mNumDofs, num_non_zeros, local_size, local_size);
    }
    else
    {
        assert(DIM==3);

        // in 3d we get the number of containing elements for each node and use that to obtain an upper bound
        // for the number of non-zeros for each DOF associated with that node.

        int* num_non_zeros_each_row = new int[mNumDofs];
        for (unsigned i=0; i<mNumDofs; i++)
        {
            num_non_zeros_each_row[i] = 0;
        }

        for (typename AbstractMesh<DIM,DIM>::NodeIterator iter = mrQuadMesh.GetNodeIteratorBegin();
                iter != mrQuadMesh.GetNodeIteratorEnd();
                ++iter)
        {
            // this upper bound neglects the fact that two containing elements will share the same nodes..
            // 4 = max num dofs associated with this node
            // 30 = 3*9+3 = 3 dimensions x 9 other nodes on this element   +  3 vertices with a pressure unknown
            unsigned num_non_zeros_upper_bound = 4 + 30*iter->GetNumContainingElements();

            num_non_zeros_upper_bound = std::min(num_non_zeros_upper_bound, mNumDofs);

            unsigned i = iter->GetIndex();

            num_non_zeros_each_row[mProblemDimension*i + 0] = num_non_zeros_upper_bound;
            num_non_zeros_each_row[mProblemDimension*i + 1] = num_non_zeros_upper_bound;
            num_non_zeros_each_row[mProblemDimension*i + 2] = num_non_zeros_upper_bound;

            if (mCompressibilityType==INCOMPRESSIBLE)
            {
                if(!iter->IsInternal())
                {
                    num_non_zeros_each_row[mProblemDimension*i + 3] = num_non_zeros_upper_bound;
                }
                else
                {
                    num_non_zeros_each_row[mProblemDimension*i + 3] = 1;
                }
            }
        }

        // NOTE: PetscTools::SetupMat() or the below creates a MATAIJ matrix, which means the matrix will
        // be of type MATSEQAIJ if num_procs=1 and MATMPIAIJ otherwise. In the former case
        // MatSeqAIJSetPreallocation MUST be called [MatMPIAIJSetPreallocation will have
        // no effect (silently)], and vice versa in the latter case

        /// We want to allocate different numbers of non-zeros per row, which means
        /// PetscTools::SetupMat isn't that useful. We could call
        //PetscTools::SetupMat(mSystemLhsMatrix, mNumDofs, mNumDofs, 0, PETSC_DECIDE, PETSC_DECIDE);
        //PetscTools::SetupMat(mPreconditionMatrix, mNumDofs, mNumDofs, 0, PETSC_DECIDE, PETSC_DECIDE);
        /// but we would get warnings due to the lack allocation

        // possible todo: create a PetscTools::SetupMatNoAllocation()


#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
        MatCreate(PETSC_COMM_WORLD,local_size,local_size,mNumDofs,mNumDofs,&mSystemLhsMatrix);
        MatCreate(PETSC_COMM_WORLD,local_size,local_size,mNumDofs,mNumDofs,&mPreconditionMatrix);
#else //New API
        MatCreate(PETSC_COMM_WORLD,&mSystemLhsMatrix);
        MatCreate(PETSC_COMM_WORLD,&mPreconditionMatrix);
        MatSetSizes(mSystemLhsMatrix,local_size,local_size,mNumDofs,mNumDofs);
        MatSetSizes(mPreconditionMatrix,local_size,local_size,mNumDofs,mNumDofs);
#endif

        if (PetscTools::IsSequential())
        {
//.........这里部分代码省略.........