C++ typenameAbstractMesh::GetRadius方法代码示例

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