GeneralisedLinearSpringForce.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
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 "GeneralisedLinearSpringForce.hpp"
00030 
00031 
00032 template<unsigned DIM>
00033 GeneralisedLinearSpringForce<DIM>::GeneralisedLinearSpringForce()
00034    : AbstractTwoBodyInteractionForce<DIM>()
00035 {
00036 }
00037 
00038 template<unsigned DIM>
00039 double GeneralisedLinearSpringForce<DIM>::VariableSpringConstantMultiplicationFactor(unsigned nodeAGlobalIndex,
00040                                                                                      unsigned nodeBGlobalIndex,
00041                                                                                      AbstractTissue<DIM>& rTissue,
00042                                                                                      bool isCloserThanRestLength)
00043 {
00044     return 1.0;
00045 }
00046 
00047 template<unsigned DIM>
00048 GeneralisedLinearSpringForce<DIM>::~GeneralisedLinearSpringForce()
00049 {
00050 }
00051 
00052 template<unsigned DIM>
00053 c_vector<double, DIM> GeneralisedLinearSpringForce<DIM>::CalculateForceBetweenNodes(unsigned nodeAGlobalIndex,
00054                                                                                     unsigned nodeBGlobalIndex,
00055                                                                                     AbstractTissue<DIM>& rTissue)
00056 {
00057     // We should only ever calculate the force between two distinct nodes
00058     assert(nodeAGlobalIndex!=nodeBGlobalIndex);
00059 
00060     // Get the node locations
00061     c_vector<double, DIM> node_a_location = rTissue.GetNode(nodeAGlobalIndex)->rGetLocation();
00062     c_vector<double, DIM> node_b_location = rTissue.GetNode(nodeBGlobalIndex)->rGetLocation();
00063 
00064     // Get the unit vector parallel to the line joining the two nodes
00065     c_vector<double, DIM> unit_difference;
00066     if (rTissue.HasMesh())
00067     {
00068         // We use the mesh method GetVectorFromAtoB() to compute the direction of the unit vector
00069         // along the line joining the two nodes, rather than simply subtract their positions,
00070         // because this method can be overloaded, e.g. to enforce a periodic boundary in Cylindrical2dMesh
00071         unit_difference = (static_cast<MeshBasedTissue<DIM>*>(&rTissue))->rGetMesh().GetVectorFromAtoB(node_a_location, node_b_location);
00072     }
00073     else
00074     {
00075         unit_difference = node_b_location - node_a_location;
00076     }
00077 
00078     // Calculate the distance between the two nodes
00079     double distance_between_nodes = norm_2(unit_difference);
00080     assert(distance_between_nodes > 0);
00081     assert(!std::isnan(distance_between_nodes));
00082 
00083     unit_difference /= distance_between_nodes;
00084 
00085     // If mUseCutoffPoint has been set, then there is zero force between
00086     // two nodes located a distance apart greater than mUseCutoffPoint
00087     if (this->mUseCutoffPoint)
00088     {
00089         if (distance_between_nodes >= TissueConfig::Instance()->GetMechanicsCutOffLength())
00090         {
00091             return zero_vector<double>(DIM); // c_vector<double,DIM>() is not guaranteed to be fresh memory
00092         }
00093     }
00094 
00095     // Calculate the rest length of the spring connecting the two nodes
00096 
00097     double rest_length = 1.0;
00098 
00099     assert( !(rTissue.IsGhostNode(nodeAGlobalIndex)) && !(rTissue.IsGhostNode(nodeBGlobalIndex)) );
00100 
00101     double ageA = rTissue.rGetCellUsingLocationIndex(nodeAGlobalIndex).GetAge();
00102     double ageB = rTissue.rGetCellUsingLocationIndex(nodeBGlobalIndex).GetAge();
00103 
00104     assert(!std::isnan(ageA));
00105     assert(!std::isnan(ageB));
00106 
00107     TissueCell& r_cell_A = rTissue.rGetCellUsingLocationIndex(nodeAGlobalIndex);
00108     TissueCell& r_cell_B = rTissue.rGetCellUsingLocationIndex(nodeBGlobalIndex);
00109 
00110     // If the cells are both newly divided, then the rest length of the spring
00111     // connecting them grows linearly with time, until 1 hour after division
00112     if ( ageA<TissueConfig::Instance()->GetMDuration() && ageB<TissueConfig::Instance()->GetMDuration() )
00113     {
00114         if (rTissue.HasMesh())
00115         {
00116             if ( (static_cast<MeshBasedTissue<DIM>*>(&rTissue))->IsMarkedSpring(r_cell_A, r_cell_B) )
00117             {
00118                 // Spring rest length increases from ???? to normal rest length, 1.0, over 1 hour
00119                 double lambda = TissueConfig::Instance()->GetDivisionRestingSpringLength();
00120                 rest_length = lambda + (1.0-lambda)*(ageA/(TissueConfig::Instance()->GetMDuration()));
00121             }
00122             if (ageA+SimulationTime::Instance()->GetTimeStep() >= TissueConfig::Instance()->GetMDuration())
00123             {
00124                 // This spring is about to go out of scope
00125                 (static_cast<MeshBasedTissue<DIM>*>(&rTissue))->UnmarkSpring(r_cell_A, r_cell_B);
00126             }
00127         }
00128         else
00129         {
00130             // Spring rest length increases from mDivisionRestingSpringLength to normal rest length, 1.0, over 1 hour
00131             double lambda = TissueConfig::Instance()->GetDivisionRestingSpringLength();
00132             rest_length = lambda + (1.0-lambda)*(ageA/(TissueConfig::Instance()->GetMDuration()));
00133         }
00134     }
00135 
00136     double a_rest_length = rest_length*0.5;
00137     double b_rest_length = a_rest_length;
00138 
00139     // If either of the cells has begun apoptosis, then the length of the spring
00140     // connecting them decreases linearly with time
00141     if (rTissue.rGetCellUsingLocationIndex(nodeAGlobalIndex).HasApoptosisBegun())
00142     {
00143         double time_until_death_a = rTissue.rGetCellUsingLocationIndex(nodeAGlobalIndex).TimeUntilDeath();
00144         a_rest_length = a_rest_length*(time_until_death_a)/(TissueConfig::Instance()->GetApoptosisTime());
00145     }
00146     if (rTissue.rGetCellUsingLocationIndex(nodeBGlobalIndex).HasApoptosisBegun())
00147     {
00148         double time_until_death_b = rTissue.rGetCellUsingLocationIndex(nodeBGlobalIndex).TimeUntilDeath();
00149         b_rest_length = b_rest_length*(time_until_death_b)/(TissueConfig::Instance()->GetApoptosisTime());
00150     }
00151 
00152     rest_length = a_rest_length + b_rest_length;
00153 
00154     assert(rest_length<=1.0+1e-12);
00155 
00156     bool is_closer_than_rest_length = true;
00157 
00158     if (distance_between_nodes - rest_length >0)
00159     {
00160         is_closer_than_rest_length = false;
00161     }
00162 
00163     // Although in this class the 'spring constant' is a constant parameter, in
00164     // subclasses it can depend on properties of each of the cells
00165     double multiplication_factor = 1.0;
00166     multiplication_factor *= VariableSpringConstantMultiplicationFactor(nodeAGlobalIndex, nodeBGlobalIndex, rTissue, is_closer_than_rest_length);
00167 
00168     if (rTissue.HasMesh())
00169     {
00170         return multiplication_factor * TissueConfig::Instance()->GetSpringStiffness() * unit_difference * (distance_between_nodes - rest_length);
00171     }
00172     else
00173     {
00174         // A reasonably stable simple force law
00175         if (distance_between_nodes > rest_length)
00176         {
00177             double alpha = 5;
00178             c_vector<double, DIM> temp = TissueConfig::Instance()->GetSpringStiffness() * unit_difference * (distance_between_nodes - rest_length)*exp(-alpha*(distance_between_nodes-rest_length));
00179             for (unsigned i=0; i<DIM; i++)
00180             {
00181                 assert(!std::isnan(temp[i]));
00182             }
00183             return temp;
00184         }
00185         else
00186         {
00187             c_vector<double, DIM> temp = TissueConfig::Instance()->GetSpringStiffness() * unit_difference * log(1 + distance_between_nodes - rest_length);
00188             for (unsigned i=0; i<DIM; i++)
00189             {
00190                 assert(!std::isnan(temp[i]));
00191             }
00192             return temp;
00193         }
00194     }
00195 }
00196 
00197 
00199 // Explicit instantiation
00201 
00202 template class GeneralisedLinearSpringForce<1>;
00203 template class GeneralisedLinearSpringForce<2>;
00204 template class GeneralisedLinearSpringForce<3>;

Generated on Tue Aug 4 16:10:21 2009 for Chaste by  doxygen 1.5.5