MeshBasedCellPopulationWithGhostNodes.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 #include "MeshBasedCellPopulationWithGhostNodes.hpp"
00029 #include "CellwiseData.hpp"
00030 
00031 template<unsigned DIM>
00032 MeshBasedCellPopulationWithGhostNodes<DIM>::MeshBasedCellPopulationWithGhostNodes(
00033      MutableMesh<DIM, DIM>& rMesh,
00034      std::vector<CellPtr>& rCells,
00035      const std::vector<unsigned> locationIndices,
00036      bool deleteMesh,
00037      double ghostSpringStiffness)
00038              : MeshBasedCellPopulation<DIM>(rMesh, rCells, locationIndices, deleteMesh, false), // do not call the base class Validate()
00039                mGhostSpringStiffness(ghostSpringStiffness)
00040 {
00041     if (!locationIndices.empty())
00042     {
00043         // Create a set of node indices corresponding to ghost nodes
00044         std::set<unsigned> node_indices;
00045         std::set<unsigned> location_indices;
00046         std::set<unsigned> ghost_node_indices;
00047 
00048         for (unsigned i=0; i<this->GetNumNodes(); i++)
00049         {
00050             node_indices.insert(this->GetNode(i)->GetIndex());
00051         }
00052         for (unsigned i=0; i<locationIndices.size(); i++)
00053         {
00054             location_indices.insert(locationIndices[i]);
00055         }
00056 
00057         std::set_difference(node_indices.begin(), node_indices.end(),
00058                             location_indices.begin(), location_indices.end(),
00059                             std::inserter(ghost_node_indices, ghost_node_indices.begin()));
00060 
00061         // This method finishes and then calls Validate()
00062         SetGhostNodes(ghost_node_indices);
00063     }
00064     else
00065     {
00066         this->mIsGhostNode = std::vector<bool>(this->GetNumNodes(), false);
00067         Validate();
00068     }
00069 }
00070 
00071 template<unsigned DIM>
00072 MeshBasedCellPopulationWithGhostNodes<DIM>::MeshBasedCellPopulationWithGhostNodes(MutableMesh<DIM, DIM>& rMesh,
00073                                                                   double ghostSpringStiffness)
00074              : MeshBasedCellPopulation<DIM>(rMesh),
00075                mGhostSpringStiffness(ghostSpringStiffness)
00076 {
00077 }
00078 
00079 template<unsigned DIM>
00080 std::vector<bool>& MeshBasedCellPopulationWithGhostNodes<DIM>::rGetGhostNodes()
00081 {
00082     return this->mIsGhostNode;
00083 }
00084 
00085 template<unsigned DIM>
00086 bool MeshBasedCellPopulationWithGhostNodes<DIM>::IsGhostNode(unsigned index)
00087 {
00088     return this->mIsGhostNode[index];
00089 }
00090 
00091 template<unsigned DIM>
00092 std::set<unsigned> MeshBasedCellPopulationWithGhostNodes<DIM>::GetGhostNodeIndices()
00093 {
00094     std::set<unsigned> ghost_node_indices;
00095     for (unsigned i=0; i<this->mIsGhostNode.size(); i++)
00096     {
00097         if (this->mIsGhostNode[i])
00098         {
00099             ghost_node_indices.insert(i);
00100         }
00101     }
00102     return ghost_node_indices;
00103 }
00104 
00105 template<unsigned DIM>
00106 void MeshBasedCellPopulationWithGhostNodes<DIM>::SetGhostNodes(const std::set<unsigned>& rGhostNodeIndices)
00107 {
00108     // Reinitialise all entries of mIsGhostNode to false
00109     this->mIsGhostNode = std::vector<bool>(this->mrMesh.GetNumNodes(), false);
00110 
00111     // Update mIsGhostNode
00112     for (std::set<unsigned>::iterator iter=rGhostNodeIndices.begin(); iter!=rGhostNodeIndices.end(); ++iter)
00113     {
00114         this->mIsGhostNode[*iter] = true;
00115     }
00116 
00117     Validate();
00118 }
00119 
00120 template<unsigned DIM>
00121 void MeshBasedCellPopulationWithGhostNodes<DIM>::UpdateGhostPositions(double dt)
00122 {
00123     // Initialise vector of forces on ghost nodes
00124     std::vector<c_vector<double, DIM> > drdt(this->GetNumNodes());
00125     for (unsigned i=0; i<drdt.size(); i++)
00126     {
00127         drdt[i] = zero_vector<double>(DIM);
00128     }
00129 
00130     // Calculate forces on ghost nodes
00131     for (typename MutableMesh<DIM, DIM>::EdgeIterator edge_iterator = this->mrMesh.EdgesBegin();
00132         edge_iterator != this->mrMesh.EdgesEnd();
00133         ++edge_iterator)
00134     {
00135         unsigned nodeA_global_index = edge_iterator.GetNodeA()->GetIndex();
00136         unsigned nodeB_global_index = edge_iterator.GetNodeB()->GetIndex();
00137 
00138         c_vector<double, DIM> force = CalculateForceBetweenGhostNodes(nodeA_global_index, nodeB_global_index);
00139 
00140         double damping_constant = this->GetDampingConstantNormal();
00141 
00142         if (!this->mIsGhostNode[nodeA_global_index])
00143         {
00144             drdt[nodeB_global_index] -= force / damping_constant;
00145         }
00146         else
00147         {
00148             drdt[nodeA_global_index] += force / damping_constant;
00149 
00150             if (this->mIsGhostNode[nodeB_global_index])
00151             {
00152                 drdt[nodeB_global_index] -= force / damping_constant;
00153             }
00154         }
00155     }
00156 
00157     for (typename AbstractMesh<DIM,DIM>::NodeIterator node_iter = this->mrMesh.GetNodeIteratorBegin();
00158          node_iter != this->mrMesh.GetNodeIteratorEnd();
00159          ++node_iter)
00160     {
00161         unsigned node_index = node_iter->GetIndex();
00162         if (this->mIsGhostNode[node_index])
00163         {
00164             ChastePoint<DIM> new_point(node_iter->rGetLocation() + dt*drdt[node_index]);
00165             this->mrMesh.SetNode(node_index, new_point, false);
00166         }
00167     }
00168 }
00169 
00170 template<unsigned DIM>
00171 c_vector<double, DIM> MeshBasedCellPopulationWithGhostNodes<DIM>::CalculateForceBetweenGhostNodes(const unsigned& rNodeAGlobalIndex, const unsigned& rNodeBGlobalIndex)
00172 {
00173     assert(rNodeAGlobalIndex != rNodeBGlobalIndex);
00174     c_vector<double, DIM> unit_difference;
00175     c_vector<double, DIM> node_a_location = this->GetNode(rNodeAGlobalIndex)->rGetLocation();
00176     c_vector<double, DIM> node_b_location = this->GetNode(rNodeBGlobalIndex)->rGetLocation();
00177 
00178     // There is reason not to subtract one position from the other (cylindrical meshes)
00179     unit_difference = this->mrMesh.GetVectorFromAtoB(node_a_location, node_b_location);
00180 
00181     double distance_between_nodes = norm_2(unit_difference);
00182 
00183     unit_difference /= distance_between_nodes;
00184 
00185     double rest_length = 1.0;
00186 
00187     return mGhostSpringStiffness * unit_difference * (distance_between_nodes - rest_length);
00188 }
00189 
00190 template<unsigned DIM>
00191 CellPtr MeshBasedCellPopulationWithGhostNodes<DIM>::AddCell(CellPtr pNewCell, const c_vector<double,DIM>& rCellDivisionVector, CellPtr pParentCell)
00192 {
00193     // Add new cell to cell population
00194     CellPtr p_created_cell = MeshBasedCellPopulation<DIM>::AddCell(pNewCell, rCellDivisionVector, pParentCell);
00195     assert(p_created_cell == pNewCell);
00196 
00197     // Update size of mIsGhostNode if necessary
00198     unsigned new_node_index = this->mCellLocationMap[p_created_cell.get()];
00199 
00200     if (this->GetNumNodes() > this->mIsGhostNode.size())
00201     {
00202         this->mIsGhostNode.resize(this->GetNumNodes());
00203         this->mIsGhostNode[new_node_index] = false;
00204     }
00205 
00206     // Return pointer to new cell
00207     return p_created_cell;
00208 }
00209 
00210 template<unsigned DIM>
00211 void MeshBasedCellPopulationWithGhostNodes<DIM>::Validate()
00212 {
00213     // Get a list of all the nodes that are ghosts
00214     std::vector<bool> validated_node = mIsGhostNode;
00215     assert(mIsGhostNode.size()==this->GetNumNodes());
00216 
00217     // Look through all of the cells and record what node they are associated with.
00218     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter=this->Begin(); cell_iter!=this->End(); ++cell_iter)
00219     {
00220         unsigned node_index = this->mCellLocationMap[(*cell_iter).get()];
00221 
00222         // If the node attached to this cell is labelled as a ghost node, then throw an error
00223         if (mIsGhostNode[node_index])
00224         {
00225             std::stringstream ss;
00226             ss << "Node " << node_index << " is labelled as a ghost node and has a cell attached";
00227             EXCEPTION(ss.str());
00228         }
00229         validated_node[node_index] = true;
00230     }
00231 
00232     for (unsigned i=0; i<validated_node.size(); i++)
00233     {
00234         if (!validated_node[i])
00235         {
00236             std::stringstream ss;
00237             ss << "Node " << i << " does not appear to be a ghost node or have a cell associated with it";
00238             EXCEPTION(ss.str());
00239         }
00240     }
00241 }
00242 
00243 template<unsigned DIM>
00244 void MeshBasedCellPopulationWithGhostNodes<DIM>::UpdateGhostNodesAfterReMesh(NodeMap& rMap)
00245 {
00246     // Copy mIsGhostNode to a temporary vector
00247     std::vector<bool> ghost_nodes_before_remesh = mIsGhostNode;
00248 
00249     // Reinitialise mIsGhostNode
00250     mIsGhostNode.clear();
00251     mIsGhostNode.resize(this->GetNumNodes());
00252 
00253     // Update mIsGhostNode using the node map
00254     for (unsigned old_index=0; old_index<rMap.Size(); old_index++)
00255     {
00256         if (!rMap.IsDeleted(old_index))
00257         {
00258             unsigned new_index = rMap.GetNewIndex(old_index);
00259             mIsGhostNode[new_index] = ghost_nodes_before_remesh[old_index];
00260         }
00261     }
00262 }
00263 
00264 template<unsigned DIM>
00265 void MeshBasedCellPopulationWithGhostNodes<DIM>::UpdateNodeLocations(const std::vector< c_vector<double, DIM> >& rNodeForces, double dt)
00266 {
00267     // First update ghost positions first because they do not affect the real cells
00268     UpdateGhostPositions(dt);
00269 
00270     // Then call the base class method
00271     AbstractCentreBasedCellPopulation<DIM>::UpdateNodeLocations(rNodeForces, dt);
00272 }
00273 
00274 template<unsigned DIM>
00275 void MeshBasedCellPopulationWithGhostNodes<DIM>::WriteVtkResultsToFile()
00276 {
00277 #ifdef CHASTE_VTK
00278     if (this->mpVoronoiTessellation != NULL)
00279     {
00280         VertexMeshWriter<DIM, DIM> mesh_writer(this->mDirPath, "results", false);
00281         std::stringstream time;
00282         time << SimulationTime::Instance()->GetTimeStepsElapsed();
00283     
00284         unsigned num_elements = this->mpVoronoiTessellation->GetNumElements();
00285         std::vector<double> ghosts(num_elements);
00286         std::vector<double> cell_types(num_elements);
00287         std::vector<double> cell_ancestors(num_elements);
00288         std::vector<double> cell_mutation_states(num_elements);
00289         std::vector<double> cell_ages(num_elements);
00290         std::vector<double> cell_cycle_phases(num_elements);
00291         std::vector<double> cell_volumes(num_elements);
00292         std::vector<std::vector<double> > cellwise_data;
00293     
00294         if (CellwiseData<DIM>::Instance()->IsSetUp())
00295         {
00296             CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00297             unsigned num_variables = p_data->GetNumVariables();
00298             for (unsigned var=0; var<num_variables; var++)
00299             {
00300                 std::vector<double> cellwise_data_var(num_elements);
00301                 cellwise_data.push_back(cellwise_data_var);
00302             }
00303         }
00304     
00305         // Loop over Voronoi elements
00306         for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = this->mpVoronoiTessellation->GetElementIteratorBegin();
00307              elem_iter != this->mpVoronoiTessellation->GetElementIteratorEnd();
00308              ++elem_iter)
00309         {
00310             // Get index of this element in the Voronoi tessellation mesh
00311             unsigned elem_index = elem_iter->GetIndex();
00312     
00313             unsigned node_index = this->mpVoronoiTessellation->GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(elem_index);
00314     
00315             ghosts[elem_index] = (double)(this->IsGhostNode(node_index));
00316     
00317             if (!this->IsGhostNode(node_index))
00318             {
00319                 // Get the cell corresponding to this element
00320                 CellPtr p_cell = this->mLocationCellMap[node_index];
00321     
00322                 if (this->mOutputCellAncestors)
00323                 {
00324                     double ancestor_index = (p_cell->GetAncestor() == UNSIGNED_UNSET) ? (-1.0) : (double)p_cell->GetAncestor();
00325                     cell_ancestors[elem_index] = ancestor_index;
00326                 }
00327                 if (this->mOutputCellProliferativeTypes)
00328                 {
00329                     double cell_type = p_cell->GetCellCycleModel()->GetCellProliferativeType();
00330                     cell_types[elem_index] = cell_type;
00331                 }
00332                 if (this->mOutputCellMutationStates)
00333                 {
00334                     double mutation_state = p_cell->GetMutationState()->GetColour();
00335                     cell_mutation_states[elem_index] = mutation_state;
00336                 }
00337                 if (this->mOutputCellAges)
00338                 {
00339                     double age = p_cell->GetAge();
00340                     cell_ages[elem_index] = age;
00341                 }
00342                 if (this->mOutputCellCyclePhases)
00343                 {
00344                     double cycle_phase = p_cell->GetCellCycleModel()->GetCurrentCellCyclePhase();
00345                     cell_cycle_phases[elem_index] = cycle_phase;
00346                 }
00347                 if (this->mOutputCellVolumes)
00348                 {
00349                     double cell_volume = this->mpVoronoiTessellation->GetVolumeOfElement(elem_index);
00350                     cell_volumes[elem_index] = cell_volume;
00351                 }
00352                 if (CellwiseData<DIM>::Instance()->IsSetUp())
00353                 {
00354                     CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00355                     unsigned num_variables = p_data->GetNumVariables();
00356                     for (unsigned var=0; var<num_variables; var++)
00357                     {
00358                         cellwise_data[var][elem_index] = p_data->GetValue(p_cell, var);
00359                     }
00360                 }
00361             }
00362             else
00363             {
00364                 if (this->mOutputCellAncestors)
00365                 {
00366                     cell_ancestors[elem_index] = -1.0;
00367                 }
00368                 if (this->mOutputCellProliferativeTypes)
00369                 {
00370                     cell_types[elem_index] = -1.0;
00371                 }
00372                 if (this->mOutputCellMutationStates)
00373                 {
00374                     cell_mutation_states[elem_index] = -1.0;
00375                 }
00376                 if (this->mOutputCellAges)
00377                 {
00378                     cell_ages[elem_index] = -1.0;
00379                 }
00380                 if (this->mOutputCellCyclePhases)
00381                 {
00382                     cell_cycle_phases[elem_index] = -1.0;
00383                 }
00384                 if (this->mOutputCellVolumes)
00385                 {
00386                     cell_volumes[elem_index] = -1.0;
00387                 }
00388             }
00389         }
00390     
00391         mesh_writer.AddCellData("Non-ghosts", ghosts);
00392         if (this->mOutputCellProliferativeTypes)
00393         {
00394             mesh_writer.AddCellData("Cell types", cell_types);
00395         }
00396         if (this->mOutputCellAncestors)
00397         {
00398             mesh_writer.AddCellData("Ancestors", cell_ancestors);
00399         }
00400         if (this->mOutputCellMutationStates)
00401         {
00402             mesh_writer.AddCellData("Mutation states", cell_mutation_states);
00403         }
00404         if (this->mOutputCellAges)
00405         {
00406             mesh_writer.AddCellData("Ages", cell_ages);
00407         }
00408         if (this->mOutputCellCyclePhases)
00409         {
00410             mesh_writer.AddCellData("Cycle phases", cell_cycle_phases);
00411         }
00412         if (this->mOutputCellVolumes)
00413         {
00414             mesh_writer.AddCellData("Cell volumes", cell_volumes);
00415         }
00416         if (CellwiseData<DIM>::Instance()->IsSetUp())
00417         {
00418             for (unsigned var=0; var<cellwise_data.size(); var++)
00419             {
00420                 std::stringstream data_name;
00421                 data_name << "Cellwise data " << var;
00422                 std::vector<double> cellwise_data_var = cellwise_data[var];
00423                 mesh_writer.AddCellData(data_name.str(), cellwise_data_var);
00424             }
00425         }
00426     
00427         mesh_writer.WriteVtkUsingMesh(*(this->mpVoronoiTessellation), time.str());
00428         *(this->mpVtkMetaFile) << "        <DataSet timestep=\"";
00429         *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
00430         *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
00431         *(this->mpVtkMetaFile) << SimulationTime::Instance()->GetTimeStepsElapsed();
00432         *(this->mpVtkMetaFile) << ".vtu\"/>\n";
00433     }
00434 #endif //CHASTE_VTK
00435 }
00436 
00437 template<unsigned DIM>
00438 void MeshBasedCellPopulationWithGhostNodes<DIM>::OutputCellPopulationParameters(out_stream& rParamsFile)
00439 {
00440     *rParamsFile << "\t\t<GhostSpringStiffness>" << mGhostSpringStiffness << "</GhostSpringStiffness> \n";
00441 
00442     // Call method on direct parent class
00443     MeshBasedCellPopulation<DIM>::OutputCellPopulationParameters(rParamsFile);
00444 }
00445 
00447 // Explicit instantiation
00449 
00450 
00451 template class MeshBasedCellPopulationWithGhostNodes<1>;
00452 template class MeshBasedCellPopulationWithGhostNodes<2>;
00453 template class MeshBasedCellPopulationWithGhostNodes<3>;
00454 
00455 
00456 // Serialization for Boost >= 1.36
00457 #include "SerializationExportWrapperForCpp.hpp"
00458 EXPORT_TEMPLATE_CLASS_SAME_DIMS(MeshBasedCellPopulationWithGhostNodes)

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5