VertexBasedCellPopulation.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 
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 so the motion has been restricted. Use a smaller timestep to avoid these warnings.");
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             {
00313                 WARNING("Cell removed due to T2Swap this is not counted in the dead cells counter");
00314                 cell_iter = this->mCells.erase(cell_iter);
00315                 --cell_iter;
00316             }
00317             else
00318             {
00319                 unsigned new_elem_index = element_map.GetNewIndex(old_elem_index);
00320 
00321                 this->mLocationCellMap[new_elem_index] = *cell_iter;
00322                 this->mCellLocationMap[(*cell_iter).get()] = new_elem_index;
00323             }
00324         }
00325 
00326         // Check that each VertexElement has only one CellPtr associated with it in the updated cell population
00327         Validate();
00328     }
00329 
00330     element_map.ResetToIdentity();
00331 }
00332 
00333 
00334 template<unsigned DIM>
00335 void VertexBasedCellPopulation<DIM>::Validate()
00336 {
00337     // Check each element has only one cell attached
00338     std::vector<unsigned> validated_element = std::vector<unsigned>(this->GetNumElements(), 0);
00339 
00340     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
00341          cell_iter != this->End();
00342          ++cell_iter)
00343     {
00344         unsigned elem_index = GetLocationIndexUsingCell(*cell_iter);
00345         validated_element[elem_index]++;
00346     }
00347 
00348     for (unsigned i=0; i<validated_element.size(); i++)
00349     {
00350         if (validated_element[i] == 0)
00351         {
00352             std::stringstream ss;
00353             ss << "Element " << i << " does not appear to have a cell associated with it";
00354             EXCEPTION(ss.str());
00355         }
00356 
00357         if (validated_element[i] > 1)
00358         {
00359             std::stringstream ss;
00360             ss << "Element " << i << " appears to have " << validated_element[i] << " cells associated with it";
00361             EXCEPTION(ss.str());
00362         }
00363     }
00364 }
00365 
00366 
00367 template<unsigned DIM>
00368 void VertexBasedCellPopulation<DIM>::WriteResultsToFiles()
00369 {
00370     AbstractCellPopulation<DIM>::WriteResultsToFiles();
00371 
00372     SimulationTime* p_time = SimulationTime::Instance();
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     // Write Locations of T3Swaps to file
00389     *mpT3SwapLocationsFile << p_time->GetTime() << "\t";
00390     std::vector< c_vector<double, DIM> > t3_swap_locations = mrMesh.GetLocationsOfT3Swaps();
00391     *mpT3SwapLocationsFile << t3_swap_locations.size() << "\t";
00392     for (unsigned index = 0;  index < t3_swap_locations.size(); index++)
00393     {
00394         for (unsigned i=0; i<DIM; i++)
00395         {
00396             *mpT3SwapLocationsFile << t3_swap_locations[index][i] << "\t";
00397         }
00398     }
00399     *mpT3SwapLocationsFile << "\n";
00400     mrMesh.ClearLocationsOfT3Swaps();
00401 
00402     // Write element data to file
00403     *mpVizElementsFile << p_time->GetTime() << "\t";
00404 
00405     // Loop over cells and find associated elements so in the same order as the cells in output files
00406     for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
00407          cell_iter != this->mCells.end();
00408          ++cell_iter)
00409     {
00410         unsigned elem_index = this->GetLocationIndexUsingCell(*cell_iter);
00411 
00412         // Hack that covers the case where the element is associated with a cell that has just been killed (#1129)
00413         bool elem_corresponds_to_dead_cell = false;
00414 
00415         if (this->mLocationCellMap[elem_index])
00416         {
00417             elem_corresponds_to_dead_cell = this->mLocationCellMap[elem_index]->IsDead();
00418         }
00419 
00420         // Write node data to file
00421         if ( !(GetElement(elem_index)->IsDeleted()) && !elem_corresponds_to_dead_cell)
00422         {
00423             VertexElement<DIM, DIM>* p_element = mrMesh.GetElement(elem_index);
00424 
00425             unsigned num_nodes_in_element = p_element->GetNumNodes();
00426 
00427             // First write the number of Nodes belonging to this VertexElement
00428             *mpVizElementsFile << num_nodes_in_element << " ";
00429 
00430             // Then write the global index of each Node in this element
00431             for (unsigned i=0; i<num_nodes_in_element; i++)
00432             {
00433                 *mpVizElementsFile << p_element->GetNodeGlobalIndex(i) << " ";
00434             }
00435         }
00436     }
00437     *mpVizElementsFile << "\n";
00438 
00439 
00440     if (this->mOutputCellVolumes)
00441     {
00442         WriteCellVolumeResultsToFile();
00443     }
00444 
00445 }
00446 
00447 template<unsigned DIM>
00448 void VertexBasedCellPopulation<DIM>::WriteCellVolumeResultsToFile()
00449 {
00450     assert(DIM==2);
00451 
00452      // Write time to file
00453     *mpCellVolumesFile << SimulationTime::Instance()->GetTime() << " ";
00454 
00455     // Loop over cells and find associated elements so in the same order as the cells in output files
00456     for (std::list<CellPtr>::iterator cell_iter = this->mCells.begin();
00457          cell_iter != this->mCells.end();
00458          ++cell_iter)
00459     {
00460         unsigned elem_index = this->GetLocationIndexUsingCell(*cell_iter);
00461 
00462         // Hack that covers the case where the element is associated with a cell that has just been killed (#1129)
00463         bool elem_corresponds_to_dead_cell = false;
00464 
00465         if (this->mLocationCellMap[elem_index])
00466         {
00467             elem_corresponds_to_dead_cell = this->mLocationCellMap[elem_index]->IsDead();
00468         }
00469 
00470         // Write node data to file
00471         if ( !(GetElement(elem_index)->IsDeleted()) && !elem_corresponds_to_dead_cell)
00472         {
00473             // Write element index to file
00474             *mpCellVolumesFile << elem_index << " ";
00475 
00476             // Write cell ID to file
00477             unsigned cell_index = (*cell_iter)->GetCellId();
00478             *mpCellVolumesFile << cell_index << " ";
00479 
00480             // Write location of element centroid to file
00481             c_vector<double, DIM> centre_location = GetLocationOfCellCentre(*cell_iter);
00482             for (unsigned i=0; i<DIM; i++)
00483             {
00484                 *mpCellVolumesFile << centre_location[i] << " ";
00485             }
00486 
00487             // Write cell volume (in 3D) or area (in 2D) to file
00488             double cell_volume = mrMesh.GetVolumeOfElement(elem_index);
00489             *mpCellVolumesFile << cell_volume << " ";
00490         }
00491     }
00492     *mpCellVolumesFile << "\n";
00493 }
00494 
00495 
00496 
00497 template<unsigned DIM>
00498 void VertexBasedCellPopulation<DIM>::WriteVtkResultsToFile()
00499 {
00500 #ifdef CHASTE_VTK
00501     SimulationTime* p_time = SimulationTime::Instance();
00502 
00503     VertexMeshWriter<DIM, DIM> mesh_writer(this->mDirPath, "results", false);
00504     std::stringstream time;
00505     time << p_time->GetTimeStepsElapsed();
00506 
00507     unsigned num_elements = mrMesh.GetNumElements();
00508     std::vector<double> cell_types(num_elements);
00509     std::vector<double> cell_ancestors(num_elements);
00510     std::vector<double> cell_mutation_states(num_elements);
00511     std::vector<double> cell_ages(num_elements);
00512     std::vector<double> cell_cycle_phases(num_elements);
00513     std::vector<double> cell_volumes(num_elements);
00514     std::vector<std::vector<double> > cellwise_data;
00515 
00516     if (CellwiseData<DIM>::Instance()->IsSetUp())
00517     {
00518         CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00519         unsigned num_variables = p_data->GetNumVariables();
00520         for (unsigned var=0; var<num_variables; var++)
00521         {
00522             std::vector<double> cellwise_data_var(num_elements);
00523             cellwise_data.push_back(cellwise_data_var);
00524         }
00525     }
00526 
00527     // Loop over vertex elements
00528     for (typename VertexMesh<DIM,DIM>::VertexElementIterator elem_iter = mrMesh.GetElementIteratorBegin();
00529          elem_iter != mrMesh.GetElementIteratorEnd();
00530          ++elem_iter)
00531     {
00532         // Get index of this element in the vertex mesh
00533         unsigned elem_index = elem_iter->GetIndex();
00534 
00535         // Get the cell corresponding to this element
00536         CellPtr p_cell = this->mLocationCellMap[elem_index];
00537         assert(p_cell);
00538 
00539         if (this->mOutputCellAncestors)
00540         {
00541             double ancestor_index = (p_cell->GetAncestor() == UNSIGNED_UNSET) ? (-1.0) : (double)p_cell->GetAncestor();
00542             cell_ancestors[elem_index] = ancestor_index;
00543         }
00544         if (this->mOutputCellProliferativeTypes)
00545         {
00546             double cell_type = p_cell->GetCellCycleModel()->GetCellProliferativeType();
00547             cell_types[elem_index] = cell_type;
00548         }
00549         if (this->mOutputCellMutationStates)
00550         {
00551             double mutation_state = p_cell->GetMutationState()->GetColour();
00552             cell_mutation_states[elem_index] = mutation_state;
00553         }
00554         if (this->mOutputCellAges)
00555         {
00556             double age = p_cell->GetAge();
00557             cell_ages[elem_index] = age;
00558         }
00559         if (this->mOutputCellCyclePhases)
00560         {
00561             double cycle_phase = p_cell->GetCellCycleModel()->GetCurrentCellCyclePhase();
00562             cell_cycle_phases[elem_index] = cycle_phase;
00563         }
00564         if (this->mOutputCellVolumes)
00565         {
00566             double cell_volume = mrMesh.GetVolumeOfElement(elem_index);
00567             cell_volumes[elem_index] = cell_volume;
00568         }
00569         if (CellwiseData<DIM>::Instance()->IsSetUp())
00570         {
00571             CellwiseData<DIM>* p_data = CellwiseData<DIM>::Instance();
00572             unsigned num_variables = p_data->GetNumVariables();
00573 
00574             for (unsigned var=0; var<num_variables; var++)
00575             {
00576                 cellwise_data[var][elem_index] = p_data->GetValue(p_cell, var);
00577             }
00578         }
00579     }
00580 
00581     if (this->mOutputCellProliferativeTypes)
00582     {
00583         mesh_writer.AddCellData("Cell types", cell_types);
00584     }
00585     if (this->mOutputCellAncestors)
00586     {
00587         mesh_writer.AddCellData("Ancestors", cell_ancestors);
00588     }
00589     if (this->mOutputCellMutationStates)
00590     {
00591         mesh_writer.AddCellData("Mutation states", cell_mutation_states);
00592     }
00593     if (this->mOutputCellAges)
00594     {
00595         mesh_writer.AddCellData("Ages", cell_ages);
00596     }
00597     if (this->mOutputCellCyclePhases)
00598     {
00599         mesh_writer.AddCellData("Cycle phases", cell_cycle_phases);
00600     }
00601     if (this->mOutputCellVolumes)
00602     {
00603         mesh_writer.AddCellData("Cell volumes", cell_volumes);
00604     }
00605     if (CellwiseData<DIM>::Instance()->IsSetUp())
00606     {
00607         for (unsigned var=0; var<cellwise_data.size(); var++)
00608         {
00609             std::stringstream data_name;
00610             data_name << "Cellwise data " << var;
00611             std::vector<double> cellwise_data_var = cellwise_data[var];
00612             mesh_writer.AddCellData(data_name.str(), cellwise_data_var);
00613         }
00614     }
00615 
00616     mesh_writer.WriteVtkUsingMesh(mrMesh, time.str());
00617     *(this->mpVtkMetaFile) << "        <DataSet timestep=\"";
00618     *(this->mpVtkMetaFile) << p_time->GetTimeStepsElapsed();
00619     *(this->mpVtkMetaFile) << "\" group=\"\" part=\"0\" file=\"results_";
00620     *(this->mpVtkMetaFile) << p_time->GetTimeStepsElapsed();
00621     *(this->mpVtkMetaFile) << ".vtu\"/>\n";
00622 #endif //CHASTE_VTK
00623 }
00624 
00625 template<unsigned DIM>
00626 void VertexBasedCellPopulation<DIM>::CreateOutputFiles(const std::string& rDirectory, bool cleanOutputDirectory)
00627 {
00628     AbstractCellPopulation<DIM>::CreateOutputFiles(rDirectory, cleanOutputDirectory);
00629 
00630     OutputFileHandler output_file_handler(rDirectory, cleanOutputDirectory);
00631     mpVizElementsFile = output_file_handler.OpenOutputFile("results.vizelements");
00632     mpT1SwapLocationsFile = output_file_handler.OpenOutputFile("T1SwapLocations.dat");
00633     mpT3SwapLocationsFile = output_file_handler.OpenOutputFile("T3SwapLocations.dat");
00634 
00635     if (this->mOutputCellVolumes)
00636     {
00637         mpCellVolumesFile = output_file_handler.OpenOutputFile("cellareas.dat");
00638     }
00639 }
00640 
00641 template<unsigned DIM>
00642 void VertexBasedCellPopulation<DIM>::CloseOutputFiles()
00643 {
00644     AbstractCellPopulation<DIM>::CloseOutputFiles();
00645     mpVizElementsFile->close();
00646     mpT1SwapLocationsFile->close();
00647     mpT3SwapLocationsFile->close();
00648 
00649     if (this->mOutputCellVolumes)
00650     {
00651         mpCellVolumesFile->close();
00652     }
00653 }
00654 
00655 
00656 template<unsigned DIM>
00657 void VertexBasedCellPopulation<DIM>::GenerateCellResultsAndWriteToFiles()
00658 {
00659     // Set up cell type counter
00660     unsigned num_cell_types = this->mCellProliferativeTypeCount.size();
00661     std::vector<unsigned> cell_type_counter(num_cell_types);
00662     for (unsigned i=0; i<num_cell_types; i++)
00663     {
00664         cell_type_counter[i] = 0;
00665     }
00666 
00667     // Set up cell cycle phase counter
00668     unsigned num_cell_cycle_phases = this->mCellCyclePhaseCount.size();
00669     std::vector<unsigned> cell_cycle_phase_counter(num_cell_cycle_phases);
00670     for (unsigned i=0; i<num_cell_cycle_phases; i++)
00671     {
00672         cell_cycle_phase_counter[i] = 0;
00673     }
00674 
00675     for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
00676          cell_iter != this->End();
00677          ++cell_iter)
00678     {
00679         this->GenerateCellResults(this->GetLocationIndexUsingCell(*cell_iter), cell_type_counter, cell_cycle_phase_counter);
00680     }
00681 
00682     this->WriteCellResultsToFiles(cell_type_counter, cell_cycle_phase_counter);
00683 }
00684 
00685 template<unsigned DIM>
00686 void VertexBasedCellPopulation<DIM>::OutputCellPopulationParameters(out_stream& rParamsFile)
00687 {
00688     *rParamsFile << "\t\t<CellRearrangementThreshold>" << mrMesh.GetCellRearrangementThreshold() << "</CellRearrangementThreshold> \n";
00689     *rParamsFile << "\t\t<T2Threshold>" <<  mrMesh.GetT2Threshold() << "</T2Threshold> \n";
00690     *rParamsFile << "\t\t<CellRearrangementRatio>" << mrMesh.GetCellRearrangementRatio() << "</CellRearrangementRatio> \n";
00691 
00692     // Call method on direct parent class
00693     AbstractCellPopulation<DIM>::OutputCellPopulationParameters(rParamsFile);
00694 }
00695 
00696 template<unsigned DIM>
00697 double VertexBasedCellPopulation<DIM>::GetWidth(const unsigned& rDimension)
00698 {
00699     // Call GetWidth() on the mesh
00700     double width = mrMesh.GetWidth(rDimension);
00701 
00702     return width;
00703 }
00704 
00706 // Explicit instantiation
00708 
00709 template class VertexBasedCellPopulation<1>;
00710 template class VertexBasedCellPopulation<2>;
00711 template class VertexBasedCellPopulation<3>;
00712 
00713 // Serialization for Boost >= 1.36
00714 #include "SerializationExportWrapperForCpp.hpp"
00715 EXPORT_TEMPLATE_CLASS_SAME_DIMS(VertexBasedCellPopulation)

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