VertexBasedCellPopulation.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
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 
00029 #include "VertexBasedCellPopulation.hpp"
00030 #include "CellwiseData.hpp"
00031 #include "VertexMeshWriter.hpp"
00032 #include "Warnings.hpp"
00033 
00034 template<unsigned DIM>
00035 VertexBasedCellPopulation<DIM>::VertexBasedCellPopulation(MutableVertexMesh<DIM, DIM>& rMesh,
00036                                           std::vector<CellPtr>& rCells,
00037                                           bool deleteMesh,
00038                                           bool validate,
00039                                           const std::vector<unsigned> locationIndices)
00040     : AbstractCellPopulation<DIM>(rCells, locationIndices),
00041       mrMesh(rMesh),
00042       mDeleteMesh(deleteMesh)
00043 {
00044     // Check the mesh contains boundary nodes
00045     bool contains_boundary_nodes = false;
00046     for (typename MutableVertexMesh<DIM,DIM>::NodeIterator node_iter = mrMesh.GetNodeIteratorBegin();
00047              node_iter != mrMesh.GetNodeIteratorEnd();
00048              ++node_iter)
00049     {
00050         if (node_iter->IsBoundaryNode())
00051         {
00052             contains_boundary_nodes = true;
00053         }
00054     }
00055     //if (mrMesh.GetNumBoundaryNodes()==0)
00057     if (!contains_boundary_nodes)
00058     {
00059         EXCEPTION("No boundary nodes are defined in the supplied vertex mesh which are needed for vertex based simulations.");
00060     }
00061 
00062     // Check each element has only one cell attached
00063      if (validate)
00064     {
00065         Validate();
00066     }
00067 }
00068 
00069 
00070 template<unsigned DIM>
00071 VertexBasedCellPopulation<DIM>::VertexBasedCellPopulation(MutableVertexMesh<DIM, DIM>& rMesh)
00072              : mrMesh(rMesh)
00073 {
00074     mDeleteMesh = true;
00075 }
00076 
00077 
00078 template<unsigned DIM>
00079 VertexBasedCellPopulation<DIM>::~VertexBasedCellPopulation()
00080 {
00081     if (mDeleteMesh)
00082     {
00083         delete &mrMesh;
00084     }
00085 }
00086 
00087 
00088 template<unsigned DIM>
00089 double VertexBasedCellPopulation<DIM>::GetDampingConstant(unsigned nodeIndex)
00090 {
00091     // Take the average of the cells containing this vertex
00092     double average_damping_constant = 0.0;
00093 
00094     std::set<unsigned> containing_elements = GetNode(nodeIndex)->rGetContainingElementIndices();
00095     double temp = 1.0/((double) containing_elements.size());
00096 
00097     for (std::set<unsigned>::iterator iter = containing_elements.begin();
00098          iter != containing_elements.end();
00099          ++iter)
00100     {
00101         CellPtr p_cell = this->GetCellUsingLocationIndex(*iter);
00102         bool cell_is_wild_type = p_cell->GetMutationState()->IsType<WildTypeCellMutationState>();
00103         bool cell_is_labelled = p_cell->HasCellProperty<CellLabel>();
00104 
00105         if (cell_is_wild_type && !cell_is_labelled)
00106         {
00107             average_damping_constant += this->GetDampingConstantNormal()*temp;
00108         }
00109         else
00110         {
00111             average_damping_constant += this->GetDampingConstantMutant()*temp;
00112         }
00113     }
00114 
00115     return average_damping_constant;
00116 }
00117 
00118 
00119 template<unsigned DIM>
00120 MutableVertexMesh<DIM, DIM>& VertexBasedCellPopulation<DIM>::rGetMesh()
00121 {
00122     return mrMesh;
00123 }
00124 
00125 
00126 template<unsigned DIM>
00127 const MutableVertexMesh<DIM, DIM>& VertexBasedCellPopulation<DIM>::rGetMesh() const
00128 {
00129     return mrMesh;
00130 }
00131 
00132 
00133 template<unsigned DIM>
00134 VertexElement<DIM, DIM>* VertexBasedCellPopulation<DIM>::GetElement(unsigned elementIndex)
00135 {
00136     return mrMesh.GetElement(elementIndex);
00137 }
00138 
00139 
00140 template<unsigned DIM>
00141 unsigned VertexBasedCellPopulation<DIM>::GetNumNodes()
00142 {
00143     return mrMesh.GetNumNodes();
00144 }
00145 
00146 
00147 template<unsigned DIM>
00148 c_vector<double, DIM> VertexBasedCellPopulation<DIM>::GetLocationOfCellCentre(CellPtr pCell)
00149 {
00150     return mrMesh.GetCentroidOfElement(this->mCellLocationMap[pCell.get()]);
00151 }
00152 
00153 
00154 template<unsigned DIM>
00155 Node<DIM>* VertexBasedCellPopulation<DIM>::GetNode(unsigned index)
00156 {
00157     return mrMesh.GetNode(index);
00158 }
00159 
00160 
00161 template<unsigned DIM>
00162 unsigned VertexBasedCellPopulation<DIM>::AddNode(Node<DIM>* pNewNode)
00163 {
00164     return mrMesh.AddNode(pNewNode);
00165 }
00166 
00167 
00168 template<unsigned DIM>
00169 void VertexBasedCellPopulation<DIM>::SetNode(unsigned nodeIndex, ChastePoint<DIM>& rNewLocation)
00170 {
00171     mrMesh.SetNode(nodeIndex, rNewLocation);
00172 }
00173 
00174 
00175 template<unsigned DIM>
00176 VertexElement<DIM, DIM>* VertexBasedCellPopulation<DIM>::GetElementCorrespondingToCell(CellPtr pCell)
00177 {
00178     return mrMesh.GetElement(this->mCellLocationMap[pCell.get()]);
00179 }
00180 
00181 
00182 template<unsigned DIM>
00183 unsigned VertexBasedCellPopulation<DIM>::GetNumElements()
00184 {
00185     return mrMesh.GetNumElements();
00186 }
00187 
00188 
00189 template<unsigned DIM>
00190 CellPtr VertexBasedCellPopulation<DIM>::AddCell(CellPtr pNewCell, const c_vector<double,DIM>& rCellDivisionVector, CellPtr pParentCell)
00191 {
00192     // Get the element associated with this cell
00193     VertexElement<DIM, DIM>* p_element = GetElementCorrespondingToCell(pParentCell);
00194 
00195     // Divide the element
00196     unsigned new_element_index;
00197     if (norm_2(rCellDivisionVector) < DBL_EPSILON)
00198     {
00199         // If the cell division vector is the default zero vector, divide the element along the short axis
00200         new_element_index = mrMesh.DivideElementAlongShortAxis(p_element, true);
00201     }
00202     else
00203     {
00204         // If the cell division vector has any non-zero component, divide the element along this axis
00205         new_element_index = mrMesh.DivideElementAlongGivenAxis(p_element, rCellDivisionVector, true);
00206     }
00207 
00208     // Associate the new cell with the element
00209     this->mCells.push_back(pNewCell);
00210 
00211     // Update location cell map
00212     CellPtr p_created_cell = this->mCells.back();
00213     this->mLocationCellMap[new_element_index] = p_created_cell;
00214     this->mCellLocationMap[p_created_cell.get()] = new_element_index;
00215     return p_created_cell;
00216 }
00217 
00218 
00219 template<unsigned DIM>
00220 unsigned VertexBasedCellPopulation<DIM>::RemoveDeadCells()
00221 {
00222     unsigned num_removed = 0;
00223 
00224     for (std::list<CellPtr>::iterator it = this->mCells.begin();
00225          it != this->mCells.end();
00226          ++it)
00227     {
00228         if ((*it)->IsDead())
00229         {
00230             // Remove the element from the mesh
00231             num_removed++;
00232             mrMesh.DeleteElementPriorToReMesh(this->mCellLocationMap[(*it).get()]);
00233             it = this->mCells.erase(it);
00234             --it;
00235         }
00236     }
00237     return num_removed;
00238 }
00239 
00240 
00241 template<unsigned DIM>
00242 void VertexBasedCellPopulation<DIM>::UpdateNodeLocations(const std::vector< c_vector<double, DIM> >& rNodeForces, double dt)
00243 {
00244     // Iterate over all nodes associated with real cells to update their positions
00245     for (unsigned node_index=0; node_index<GetNumNodes(); node_index++)
00246     {
00247         // Get the damping constant for this node
00248         double damping_const = this->GetDampingConstant(node_index);
00249 
00250         // Compute the displacement of this node
00251         c_vector<double, DIM> displacement = dt*rNodeForces[node_index]/damping_const;
00252 
00253         /*
00254          * If the displacement of this node is greater than half the cell rearrangement threshold,
00255          * this could result in nodes moving into the interior of other elements, which should not
00256          * be possible. Therefore in this case we restrict the displacement of the node to the cell
00257          * rearrangement threshold and warn the user that a smaller timestep should be used. This
00258          * restriction ensures that vertex elements remain well defined (see #1376).
00259          */
00260         if (norm_2(displacement)>0.5*mrMesh.GetCellRearrangementThreshold())
00261         {
00262             WARN_ONCE_ONLY("Vertices are moving more than half the CellRearrangementThreshold this could cause elements to become inverted the motion has been restricted: - To avoid these warnings use a smaller timestep");
00263             displacement *= 0.5*mrMesh.GetCellRearrangementThreshold()/norm_2(displacement);
00264         }
00265 
00266         // Get new node location
00267         c_vector<double, DIM> new_node_location = this->GetNode(node_index)->rGetLocation() + displacement;
00268 
00269         // Create ChastePoint for new node location
00270         ChastePoint<DIM> new_point(new_node_location);
00271 
00272         // Move the node
00273         this->SetNode(node_index, new_point);
00274     }
00275 }
00276 
00277 
00278 template<unsigned DIM>
00279 bool VertexBasedCellPopulation<DIM>::IsCellAssociatedWithADeletedLocation(CellPtr pCell)
00280 {
00281     return GetElementCorrespondingToCell(pCell)->IsDeleted();;
00282 }
00283 
00284 
00285 template<unsigned DIM>
00286 void VertexBasedCellPopulation<DIM>::Update(bool hasHadBirthsOrDeaths)
00287 {
00288     VertexElementMap element_map(mrMesh.GetNumAllElements());
00289 
00290     mrMesh.ReMesh(element_map);
00291 
00292     if (!element_map.IsIdentityMap())
00293     {
00294         // Fix up the mappings between CellPtrs and VertexElements
00295         std::map<Cell*, unsigned> old_map = this->mCellLocationMap;
00296 
00297         this->mCellLocationMap.clear();
00298         this->mLocationCellMap.clear();
00299 
00300         for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
00301              cell_iter != this->mCells.end();
00302              ++cell_iter)
00303         {
00304             // This shouldn't ever happen, as the cell vector only contains living cells
00305             unsigned old_elem_index = old_map[(*cell_iter).get()];
00306 
00307             if (element_map.IsDeleted(old_elem_index))
00308             {
00309                 /*\todo this is a kludge to remove the cell once a T2Swap occurs this is not included in the dead cells counter.
00310                  * This should be included in the RemoveDeadCells method so the death is counted
00311                  */
00312                 WARNING("Cell removed due to T2Swap this is not counted in the dead cells counter");
00313                 cell_iter = this->mCells.erase(cell_iter);
00314                 --cell_iter;
00315             }
00316             else
00317             {
00318                 unsigned new_elem_index = element_map.GetNewIndex(old_elem_index);
00319 
00320                 this->mLocationCellMap[new_elem_index] = *cell_iter;
00321                 this->mCellLocationMap[(*cell_iter).get()] = new_elem_index;
00322             }
00323         }
00324 
00325         // Check that each VertexElement has only one CellPtr associated with it in the updated cell population
00326         Validate();
00327     }
00328 
00329     element_map.ResetToIdentity();
00330 }
00331 
00332 
00333 template<unsigned DIM>
00334 void VertexBasedCellPopulation<DIM>::Validate()
00335 {
00336     // Check each element has only one cell attached
00337     std::vector<unsigned> validated_element = std::vector<unsigned>(this->GetNumElements(), 0);
00338 
00339     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
00340          cell_iter != this->End();
00341          ++cell_iter)
00342     {
00343         unsigned elem_index = GetLocationIndexUsingCell(*cell_iter);
00344         validated_element[elem_index]++;
00345     }
00346 
00347     for (unsigned i=0; i<validated_element.size(); i++)
00348     {
00349         if (validated_element[i] == 0)
00350         {
00351             std::stringstream ss;
00352             ss << "Element " << i << " does not appear to have a cell associated with it";
00353             EXCEPTION(ss.str());
00354         }
00355 
00356         if (validated_element[i] > 1)
00357         {
00358             std::stringstream ss;
00359             ss << "Element " << i << " appears to have " << validated_element[i] << " cells associated with it";
00360             EXCEPTION(ss.str());
00361         }
00362     }
00363 }
00364 
00365 
00366 template<unsigned DIM>
00367 void VertexBasedCellPopulation<DIM>::WriteResultsToFiles()
00368 {
00369     AbstractCellPopulation<DIM>::WriteResultsToFiles();
00370 
00371     SimulationTime* p_time = SimulationTime::Instance();
00372 
00373 
00374     // Write Locations of T1Swaps to file
00375     *mpT1SwapLocationsFile << p_time->GetTime() << "\t";
00376     std::vector< c_vector<double, DIM> > t1_swap_locations = mrMesh.GetLocationsOfT1Swaps();
00377     *mpT1SwapLocationsFile << t1_swap_locations.size() << "\t";
00378     for (unsigned index = 0;  index < t1_swap_locations.size(); index++)
00379     {
00380         for (unsigned i=0; i<DIM; i++)
00381         {
00382             *mpT1SwapLocationsFile << t1_swap_locations[index][i] << "\t";
00383         }
00384     }
00385     *mpT1SwapLocationsFile << "\n";
00386     mrMesh.ClearLocationsOfT1Swaps();
00387 
00388 
00389     // Write Locations of T3Swaps to file
00390     *mpT3SwapLocationsFile << p_time->GetTime() << "\t";
00391     std::vector< c_vector<double, DIM> > t3_swap_locations = mrMesh.GetLocationsOfT3Swaps();
00392     *mpT3SwapLocationsFile << t3_swap_locations.size() << "\t";
00393     for (unsigned index = 0;  index < t3_swap_locations.size(); index++)
00394     {
00395         for (unsigned i=0; i<DIM; i++)
00396         {
00397             *mpT3SwapLocationsFile << t3_swap_locations[index][i] << "\t";
00398         }
00399     }
00400     *mpT3SwapLocationsFile << "\n";
00401     mrMesh.ClearLocationsOfT3Swaps();
00402 
00403 
00404     // Write element data to file
00405     *mpVizElementsFile << p_time->GetTime() << "\t";
00406 
00407     // Loop over cells and find associated elements so in the same order as the cells in output files
00408     for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
00409          cell_iter != this->mCells.end();
00410          ++cell_iter)
00411     {
00412         unsigned elem_index = this->GetLocationIndexUsingCell(*cell_iter);
00413 
00414         // Hack that covers the case where the element is associated with a cell that has just been killed (#1129)
00415         bool elem_corresponds_to_dead_cell = false;
00416 
00417         if (this->mLocationCellMap[elem_index])
00418         {
00419             elem_corresponds_to_dead_cell = this->mLocationCellMap[elem_index]->IsDead();
00420         }
00421 
00422         // Write node data to file
00423         if ( !(GetElement(elem_index)->IsDeleted()) && !elem_corresponds_to_dead_cell)
00424         {
00425             VertexElement<DIM, DIM>* p_element = mrMesh.GetElement(elem_index);
00426 
00427             unsigned num_nodes_in_element = p_element->GetNumNodes();
00428 
00429             // First write the number of Nodes belonging to this VertexElement
00430             *mpVizElementsFile << num_nodes_in_element << " ";
00431 
00432             // Then write the global index of each Node in this element
00433             for (unsigned i=0; i<num_nodes_in_element; i++)
00434             {
00435                 *mpVizElementsFile << p_element->GetNodeGlobalIndex(i) << " ";
00436             }
00437         }
00438     }
00439     *mpVizElementsFile << "\n";
00440 
00441 #ifdef CHASTE_VTK
00442     VertexMeshWriter<DIM, DIM> mesh_writer(mDirPath, "results", false);
00443     std::stringstream time;
00444     time << p_time->GetTimeStepsElapsed();
00445 
00446     unsigned num_elements = mrMesh.GetNumElements();
00447     std::vector<double> cell_types(num_elements);
00448     std::vector<double> cell_ancestors(num_elements);
00449     std::vector<double> cell_mutation_states(num_elements);
00450     std::vector<double> cell_ages(num_elements);
00451     std::vector<double> cell_cycle_phases(num_elements);
00452     std::vector<double> cell_volumes(num_elements);
00453     std::vector<std::vector<double> > cellwise_data;
00454 
00455     if (CellwiseData<DIM>::Instance()->IsSetUp())
00456     {
00457         CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00458         unsigned num_variables = p_data->GetNumVariables();
00459         for (unsigned var=0; var<num_variables; var++)
00460         {
00461             std::vector<double> cellwise_data_var(num_elements);
00462             cellwise_data.push_back(cellwise_data_var);
00463         }
00464     }
00465 
00466     // Loop over vertex elements
00467     for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = mrMesh.GetElementIteratorBegin();
00468          elem_iter != mrMesh.GetElementIteratorEnd();
00469          ++elem_iter)
00470     {
00471         // Get index of this element in the vertex mesh
00472         unsigned elem_index = elem_iter->GetIndex();
00473 
00474         // Get the cell corresponding to this element
00475         CellPtr p_cell = this->mLocationCellMap[elem_index];
00476         assert(p_cell);
00477 
00478         if (this->mOutputCellAncestors)
00479         {
00480             double ancestor_index = (p_cell->GetAncestor() == UNSIGNED_UNSET) ? (-1.0) : (double)p_cell->GetAncestor();
00481             cell_ancestors[elem_index] = ancestor_index;
00482         }
00483         if (this->mOutputCellProliferativeTypes)
00484         {
00485             double cell_type = p_cell->GetCellCycleModel()->GetCellProliferativeType();
00486             cell_types[elem_index] = cell_type;
00487         }
00488         if (this->mOutputCellMutationStates)
00489         {
00490             double mutation_state = p_cell->GetMutationState()->GetColour();
00491             cell_mutation_states[elem_index] = mutation_state;
00492         }
00493         if (this->mOutputCellAges)
00494         {
00495             double age = p_cell->GetAge();
00496             cell_ages[elem_index] = age;
00497         }
00498         if (this->mOutputCellCyclePhases)
00499         {
00500             double cycle_phase = p_cell->GetCellCycleModel()->GetCurrentCellCyclePhase();
00501             cell_cycle_phases[elem_index] = cycle_phase;
00502         }
00503         if (this->mOutputCellVolumes)
00504         {
00505             double cell_volume = mrMesh.GetVolumeOfElement(elem_index);
00506             cell_volumes[elem_index] = cell_volume;
00507         }
00508         if (CellwiseData<DIM>::Instance()->IsSetUp())
00509         {
00510             CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00511             unsigned num_variables = p_data->GetNumVariables();
00512 
00513             for (unsigned var=0; var<num_variables; var++)
00514             {
00515                 cellwise_data[var][elem_index] = p_data->GetValue(p_cell, var);
00516             }
00517         }
00518     }
00519 
00520     if (this->mOutputCellProliferativeTypes)
00521     {
00522         mesh_writer.AddCellData("Cell types", cell_types);
00523     }
00524     if (this->mOutputCellAncestors)
00525     {
00526         mesh_writer.AddCellData("Ancestors", cell_ancestors);
00527     }
00528     if (this->mOutputCellMutationStates)
00529     {
00530         mesh_writer.AddCellData("Mutation states", cell_mutation_states);
00531     }
00532     if (this->mOutputCellAges)
00533     {
00534         mesh_writer.AddCellData("Ages", cell_ages);
00535     }
00536     if (this->mOutputCellCyclePhases)
00537     {
00538         mesh_writer.AddCellData("Cycle phases", cell_cycle_phases);
00539     }
00540     if (this->mOutputCellVolumes)
00541     {
00542         mesh_writer.AddCellData("Cell volumes", cell_volumes);
00543     }
00544     if (CellwiseData<DIM>::Instance()->IsSetUp())
00545     {
00546         for (unsigned var=0; var<cellwise_data.size(); var++)
00547         {
00548             std::stringstream data_name;
00549             data_name << "Cellwise data " << var;
00550             std::vector<double> cellwise_data_var = cellwise_data[var];
00551             mesh_writer.AddCellData(data_name.str(), cellwise_data_var);
00552         }
00553     }
00554 
00555     mesh_writer.WriteVtkUsingMesh(mrMesh, time.str());
00556     *mpVtkMetaFile << "        <DataSet timestep=\"";
00557     *mpVtkMetaFile << p_time->GetTimeStepsElapsed();
00558     *mpVtkMetaFile << "\" group=\"\" part=\"0\" file=\"results_";
00559     *mpVtkMetaFile << p_time->GetTimeStepsElapsed();
00560     *mpVtkMetaFile << ".vtu\"/>\n";
00561 #endif //CHASTE_VTK
00562 }
00563 
00564 
00565 template<unsigned DIM>
00566 void VertexBasedCellPopulation<DIM>::CreateOutputFiles(const std::string& rDirectory, bool cleanOutputDirectory)
00567 {
00568     AbstractCellPopulation<DIM>::CreateOutputFiles(rDirectory, cleanOutputDirectory);
00569 
00570     OutputFileHandler output_file_handler(rDirectory, cleanOutputDirectory);
00571     mpVizElementsFile = output_file_handler.OpenOutputFile("results.vizelements");
00572     mpT1SwapLocationsFile = output_file_handler.OpenOutputFile("T1SwapLocations.dat");
00573     mpT3SwapLocationsFile = output_file_handler.OpenOutputFile("T3SwapLocations.dat");
00574     mDirPath = rDirectory;
00575 #ifdef CHASTE_VTK
00576     mpVtkMetaFile = output_file_handler.OpenOutputFile("results.pvd");
00577     *mpVtkMetaFile << "<?xml version=\"1.0\"?>\n";
00578     *mpVtkMetaFile << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n";
00579     *mpVtkMetaFile << "    <Collection>\n";
00580 #endif //CHASTE_VTK
00581 }
00582 
00583 
00584 template<unsigned DIM>
00585 void VertexBasedCellPopulation<DIM>::CloseOutputFiles()
00586 {
00587     AbstractCellPopulation<DIM>::CloseOutputFiles();
00588     mpVizElementsFile->close();
00589     mpT1SwapLocationsFile->close();
00590     mpT3SwapLocationsFile->close();
00591 #ifdef CHASTE_VTK
00592     *mpVtkMetaFile << "    </Collection>\n";
00593     *mpVtkMetaFile << "</VTKFile>\n";
00594     mpVtkMetaFile->close();
00595 #endif //CHASTE_VTK
00596 }
00597 
00598 
00599 template<unsigned DIM>
00600 void VertexBasedCellPopulation<DIM>::GenerateCellResultsAndWriteToFiles()
00601 {
00602     // Set up cell type counter
00603     unsigned num_cell_types = this->mCellProliferativeTypeCount.size();
00604     std::vector<unsigned> cell_type_counter(num_cell_types);
00605     for (unsigned i=0; i<num_cell_types; i++)
00606     {
00607         cell_type_counter[i] = 0;
00608     }
00609 
00610     // Set up cell cycle phase counter
00611     unsigned num_cell_cycle_phases = this->mCellCyclePhaseCount.size();
00612     std::vector<unsigned> cell_cycle_phase_counter(num_cell_cycle_phases);
00613     for (unsigned i=0; i<num_cell_cycle_phases; i++)
00614     {
00615         cell_cycle_phase_counter[i] = 0;
00616     }
00617 
00618     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
00619          cell_iter != this->End();
00620          ++cell_iter)
00621     {
00622         this->GenerateCellResults(this->GetLocationIndexUsingCell(*cell_iter), cell_type_counter, cell_cycle_phase_counter);
00623     }
00624 
00625     this->WriteCellResultsToFiles(cell_type_counter, cell_cycle_phase_counter);
00626 }
00627 
00628 template<unsigned DIM>
00629 void VertexBasedCellPopulation<DIM>::OutputCellPopulationParameters(out_stream& rParamsFile)
00630 {
00631     *rParamsFile <<  "\t\t<CellRearrangementThreshold>"<<  mrMesh.GetCellRearrangementThreshold() << "</CellRearrangementThreshold> \n" ;
00632     *rParamsFile <<  "\t\t<T2Threshold>"<<  mrMesh.GetT2Threshold() << "</T2Threshold> \n" ;
00633     *rParamsFile <<  "\t\t<CellRearrangementRatio>"<<  mrMesh.GetCellRearrangementRatio() << "</CellRearrangementRatio> \n" ;
00634 
00635     // Call direct parent class method
00636     AbstractCellPopulation<DIM>::OutputCellPopulationParameters(rParamsFile);
00637 
00638 }
00639 
00640 template<unsigned DIM>
00641 double VertexBasedCellPopulation<DIM>::GetWidth(const unsigned& rDimension)
00642 {
00643     // Call GetWidth() on the mesh
00644     double width = mrMesh.GetWidth(rDimension);
00645 
00646     return width;
00647 }
00648 
00650 // Explicit instantiation
00652 
00653 template class VertexBasedCellPopulation<1>;
00654 template class VertexBasedCellPopulation<2>;
00655 template class VertexBasedCellPopulation<3>;
00656 
00657 // Serialization for Boost >= 1.36
00658 #include "SerializationExportWrapperForCpp.hpp"
00659 EXPORT_TEMPLATE_CLASS_SAME_DIMS(VertexBasedCellPopulation)

Generated on Mon Nov 1 12:35:16 2010 for Chaste by  doxygen 1.5.5