C++ typenameAbstractMesh类代码示例

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