Chaste  Release::3.4
NodeBasedCellPopulationWithBuskeUpdate.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 #include "NodeBasedCellPopulationWithBuskeUpdate.hpp"
36 
37 #include "ReplicatableVector.hpp"
38 #include "OdeLinearSystemSolver.hpp"
39 
40 template<unsigned DIM>
42  std::vector<CellPtr>& rCells,
43  const std::vector<unsigned> locationIndices,
44  bool deleteMesh)
45  : NodeBasedCellPopulation<DIM>(rMesh, rCells, locationIndices, deleteMesh)
46 {
47 }
48 
49 template<unsigned DIM>
51  : NodeBasedCellPopulation<DIM>(rMesh)
52 {
53  // No Validate() because the cells are not associated with the cell population yet in archiving
54 }
55 
56 template<unsigned DIM>
58 {
59  // Declare solver and give the size of the system and timestep
60  unsigned system_size = this->GetNumNodes()*DIM;
61 
62  OdeLinearSystemSolver solver(system_size, dt);
63 
64  // Set up the matrix
65  Mat& r_matrix = solver.rGetLhsMatrix();
66 
67  // Initial condition
68  Vec initial_condition = PetscTools::CreateAndSetVec(system_size, 0.0);
69 
70  // Then an rGetForceVector for RHS
71  Vec& r_vector = solver.rGetForceVector();
72 
73  // Iterate over all nodes associated with real cells to construct the matrix A.
74  for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
75  cell_iter != this->End();
76  ++cell_iter)
77  {
78  // Get index of node associated with cell
79  unsigned global_node_index = this->GetLocationIndexUsingCell((*cell_iter));
80 
81  // Get the local index using the mesh
82  unsigned node_index = this->rGetMesh().SolveNodeMapping(global_node_index);
83 
84  Node<DIM>* p_node_i = this->GetNode(global_node_index);
85 
86  // Get the location of this node
87  c_vector<double, DIM> node_i_location = p_node_i->rGetLocation();
88 
89  // Get the radius of this cell
90  double radius_of_cell_i = p_node_i->GetRadius();
91 
92  // Get damping constant for node
93  double damping_const = this->GetDampingConstant(global_node_index);
94 
95  // loop over neighbours to add contribution
96 
97  // Get the set of node indices corresponding to this cell's neighbours
98  std::set<unsigned> neighbouring_node_indices = this->GetNeighbouringNodeIndices(global_node_index);
99 
100  for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
101  iter != neighbouring_node_indices.end();
102  ++iter)
103  {
104  unsigned neighbour_node_global_index = *iter;
105 
106  unsigned neighbour_node_index = this->rGetMesh().SolveNodeMapping(neighbour_node_global_index);
107 
108  // Calculate Aij
109  double Aij = 0.0;
110 
111  Node<DIM>* p_node_j = this->GetNode(neighbour_node_global_index);
112 
113  // Get the location of this node
114  c_vector<double, DIM> node_j_location = p_node_j->rGetLocation();
115 
116  // Get the unit vector parallel to the line joining the two nodes (assuming no periodicities etc.)
117  c_vector<double, DIM> unit_vector = node_j_location - node_i_location;
118 
119  // Calculate the distance between the two nodes
120  double dij = norm_2(unit_vector);
121 
122  unit_vector /= dij;
123 
124  // Get the radius of the cell corresponding to this node
125  double radius_of_cell_j = p_node_j->GetRadius();
126 
127  if (dij < radius_of_cell_i + radius_of_cell_j)
128  {
129  // ...then compute the adhesion force and add it to the vector of forces...
130  double xij = 0.5*(radius_of_cell_i*radius_of_cell_i - radius_of_cell_j*radius_of_cell_j + dij*dij)/dij;
131 
132  Aij = M_PI*(radius_of_cell_i*radius_of_cell_i - xij*xij);
133 
134  // This is contribution from the sum term in (A7)
135  for (unsigned i=0; i<DIM; i++)
136  {
137  PetscMatTools::AddToElement(r_matrix, DIM*neighbour_node_index+i, DIM*neighbour_node_index+i, -damping_const*Aij);
138  PetscMatTools::AddToElement(r_matrix, DIM*node_index+i, DIM*node_index+i, damping_const*Aij);
139  }
140  }
141  }
142 
143  // This is the standard contribution (i.e. not in the sum) in (A7)
144  for (unsigned i=0; i<DIM; i++)
145  {
146  PetscMatTools::AddToElement(r_matrix, DIM*node_index+i, DIM*node_index+i, damping_const);
147  }
148 
149  // Add current positions to initial_conditions and RHS vector
150 
151  // Note that we define these vectors before setting them as otherwise the profiling build will break (see #2367)
152  c_vector<double, DIM> current_location;
153  c_vector<double, DIM> forces;
154  current_location = this->GetNode(global_node_index)->rGetLocation();
155  forces = this->GetNode(global_node_index)->rGetAppliedForce();
156 
157  for (unsigned i=0; i<DIM; i++)
158  {
159  PetscVecTools::SetElement(initial_condition, DIM*node_index+i, current_location(i));
160  PetscVecTools::SetElement(r_vector, DIM*node_index+i, forces(i));
161  }
162  }
163  PetscMatTools::Finalise(r_matrix);
164 
165  solver.SetInitialConditionVector(initial_condition);
166 
167  // Solve to get solution at next timestep
168  Vec soln_next_timestep = solver.SolveOneTimeStep();
169 
170  ReplicatableVector soln_next_timestep_repl(soln_next_timestep);
171 
172  // Iterate over all nodes associated with real cells to update the node locations
173  for (typename AbstractCellPopulation<DIM>::Iterator cell_iter = this->Begin();
174  cell_iter != this->End();
175  ++cell_iter)
176  {
177  // Get index of node associated with cell
178  unsigned global_node_index = this->GetLocationIndexUsingCell((*cell_iter));
179 
180  unsigned node_index = this->rGetMesh().SolveNodeMapping(global_node_index);
181 
182  c_vector<double, DIM> new_node_location;
183 
184  // Get new node location
185  for (unsigned i=0; i<DIM; i++)
186  {
187  new_node_location(i) = soln_next_timestep_repl[DIM*node_index+i];
188  }
189 
190  // Create ChastePoint for new node location
191  ChastePoint<DIM> new_point(new_node_location);
192 
193  // Move the node
194  this->SetNode(global_node_index, new_point);
195  }
196 
197  // Tidy up
198  PetscTools::Destroy(initial_condition);
199 }
200 
201 template<unsigned DIM>
203 {
204  // Currently no specific parameters to output all come from parent classes
205 
206  // Call method on direct parent class
208 }
209 
214 
215 // Serialization for Boost >= 1.36
void OutputCellPopulationParameters(out_stream &rParamsFile)
Definition: Node.hpp:58
static void SetElement(Vec vector, PetscInt row, double value)
static void AddToElement(Mat matrix, PetscInt row, PetscInt col, double value)
NodeBasedCellPopulationWithBuskeUpdate(NodesOnlyMesh< DIM > &rMesh, std::vector< CellPtr > &rCells, const std::vector< unsigned > locationIndices=std::vector< unsigned >(), bool deleteMesh=false)
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
static Vec CreateAndSetVec(int size, double value)
Definition: PetscTools.cpp:254
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:351
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition: Node.cpp:140
static void Finalise(Mat matrix)
double GetRadius()
Definition: Node.cpp:249
void SetInitialConditionVector(Vec initialConditionsVector)