DiscreteSystemForceCalculator.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 "DiscreteSystemForceCalculator.hpp"
00029 
00030 DiscreteSystemForceCalculator::DiscreteSystemForceCalculator(MeshBasedCellPopulation<2>& rCellPopulation,
00031                                                              std::vector<AbstractTwoBodyInteractionForce<2>*> forceCollection)
00032     : mrCellPopulation(rCellPopulation),
00033       mForceCollection(forceCollection),
00034       mEpsilon(0.01)
00035 {
00036 }
00037 
00038 
00039 std::vector< std::vector<double> > DiscreteSystemForceCalculator::CalculateExtremalNormalForces()
00040 {
00041     unsigned num_nodes = mrCellPopulation.GetNumNodes();
00042 
00043     std::vector< std::vector<double> > extremal_normal_forces;
00044     std::vector<double> minimum_normal_forces(num_nodes);
00045     std::vector<double> maximum_normal_forces(num_nodes);
00046 
00047     for (unsigned i=0; i<num_nodes; i++)
00048     {
00049         std::vector<double> sampling_angles = GetSamplingAngles(i);
00050         std::vector<double> extremal_angles = GetExtremalAngles(i, sampling_angles);
00051 
00052         double minimum_normal_force_for_node_i = DBL_MAX;
00053         double maximum_normal_force_for_node_i = -DBL_MAX;
00054 
00055         for (unsigned j=0; j<extremal_angles.size(); j++)
00056         {
00057             double current_normal_force = CalculateFtAndFn(i, extremal_angles[j])[1];
00058 
00059             if (current_normal_force > maximum_normal_force_for_node_i)
00060             {
00061                 maximum_normal_force_for_node_i = current_normal_force;
00062             }
00063             if (current_normal_force < minimum_normal_force_for_node_i)
00064             {
00065                 minimum_normal_force_for_node_i = current_normal_force;
00066             }
00067         }
00068 
00069         assert( minimum_normal_force_for_node_i <= maximum_normal_force_for_node_i);
00070 
00071         minimum_normal_forces[i] = minimum_normal_force_for_node_i;
00072         maximum_normal_forces[i] = maximum_normal_force_for_node_i;
00073     }
00074 
00075     extremal_normal_forces.push_back(minimum_normal_forces);
00076     extremal_normal_forces.push_back(maximum_normal_forces);
00077 
00078     return extremal_normal_forces;
00079 }
00080 
00081 
00082 void DiscreteSystemForceCalculator::WriteResultsToFile(std::string simulationOutputDirectory)
00083 {
00084     double time = SimulationTime::Instance()->GetTime();
00085     std::ostringstream time_string;
00086     time_string << time;
00087     std::string results_directory = simulationOutputDirectory + "/results_from_time_" + time_string.str();
00088 
00089     OutputFileHandler output_file_handler2(results_directory+"/", false);
00090     mpVizStressResultsFile = output_file_handler2.OpenOutputFile("results.vizstress");
00091 
00092     (*mpVizStressResultsFile) <<  time << "\t";
00093 
00094     double global_index;
00095     double x;
00096     double y;
00097     double minimum;
00098     double maximum;
00099 
00100     TetrahedralMesh<2,2>& r_mesh = mrCellPopulation.rGetMesh();
00101 
00102     std::vector< std::vector<double> > extremal_normal_forces = CalculateExtremalNormalForces();
00103 
00104     for (unsigned i=0; i<r_mesh.GetNumNodes(); i++)
00105     {
00106         global_index = (double) i;
00107 
00108         x = r_mesh.GetNode(i)->rGetLocation()[0];
00109         y = r_mesh.GetNode(i)->rGetLocation()[1];
00110 
00111         minimum = extremal_normal_forces[0][i];
00112         maximum = extremal_normal_forces[1][i];
00113 
00114         (*mpVizStressResultsFile) << global_index << " " << x << " " << y << " " << minimum << " " << maximum << " ";
00115     }
00116 
00117     (*mpVizStressResultsFile) << "\n";
00118     mpVizStressResultsFile->close();
00119 }
00120 
00121 
00122 std::set<unsigned> DiscreteSystemForceCalculator::GetNeighbouringNodeIndices(unsigned index)
00123 {
00124     TetrahedralMesh<2,2>& r_mesh = mrCellPopulation.rGetMesh();
00125 
00126     Node<2>* p_node = r_mesh.GetNode(index);
00127 
00128     std::set<unsigned> neighbouring_node_indices;
00129 
00130     for (Node<2>::ContainingElementIterator it = p_node->ContainingElementsBegin();
00131          it != p_node->ContainingElementsEnd();
00132          ++it)
00133     {
00134         Element<2,2>* p_element = r_mesh.GetElement(*it);
00135         for (unsigned i=0; i<p_element->GetNumNodes(); i++)
00136         {
00137             unsigned node_index = p_element->GetNodeGlobalIndex(i);
00138             if (node_index!=index)
00139             {
00140                 neighbouring_node_indices.insert(node_index);
00141             }
00142         }
00143     }
00144     return neighbouring_node_indices;
00145 }
00146 
00147 
00148 std::vector<double> DiscreteSystemForceCalculator::CalculateFtAndFn(unsigned index, double theta)
00149 {
00150     TetrahedralMesh<2,2>& r_mesh = mrCellPopulation.rGetMesh();
00151 
00152     std::set<unsigned> neighbouring_node_indices = GetNeighbouringNodeIndices(index);
00153 
00154     double tangential_force = 0.0;
00155     double normal_force = 0.0;
00156     double alpha;
00157 
00158     c_vector<double,2> unit_vec_between_nodes(2);
00159 
00160     for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
00161          iter != neighbouring_node_indices.end();
00162          ++iter)
00163     {
00164         // The method GetAngleBetweenNodes() returns an angle in the range (-pi,pi]
00165         alpha = r_mesh.GetAngleBetweenNodes(index, *iter);
00166 
00167         assert(alpha <= M_PI);
00168         assert(alpha > -M_PI);
00169 
00170         if ( sin(alpha-theta) > DBL_EPSILON )
00171         {
00172             // Initialise a zero force vector between neighbouring nodes
00173             c_vector<double,2> force_between_nodes = zero_vector<double>(2);
00174 
00175             // Iterate over vector of forces present and add up forces between nodes
00176             for (std::vector<AbstractTwoBodyInteractionForce<2>*>::iterator force_iter = mForceCollection.begin();
00177                  force_iter != mForceCollection.end();
00178                  ++force_iter)
00179             {
00180                force_between_nodes += (*force_iter)->CalculateForceBetweenNodes(index, *iter, mrCellPopulation);
00181             }
00182 
00183             unit_vec_between_nodes[0] = cos(alpha);
00184             unit_vec_between_nodes[1] = sin(alpha);
00185 
00186             double plusminus_norm_force = inner_prod(force_between_nodes,unit_vec_between_nodes);
00187             tangential_force += plusminus_norm_force * cos(alpha-theta);
00188             normal_force += plusminus_norm_force * sin(alpha-theta);
00189         }
00190     }
00191 
00192     std::vector<double> ret(2);
00193     ret[0] = tangential_force;
00194     ret[1] = normal_force;
00195 
00196     return ret;
00197 }
00198 
00199 
00200 std::vector<double> DiscreteSystemForceCalculator::GetSamplingAngles(unsigned index)
00201 {
00202     TetrahedralMesh<2,2>& r_mesh = mrCellPopulation.rGetMesh();
00203 
00204     std::set<unsigned> neighbouring_node_indices = GetNeighbouringNodeIndices(index);
00205 
00206     std::vector<double> sampling_angles(4*neighbouring_node_indices.size());
00207 
00208     unsigned i=0;
00209 
00210     for (std::set<unsigned>::iterator iter = neighbouring_node_indices.begin();
00211          iter != neighbouring_node_indices.end();
00212          ++iter)
00213     {
00214         // The method GetAngleBetweenNodes() returns an angle in the range (-pi,pi]
00215         double alpha = r_mesh.GetAngleBetweenNodes(index, *iter);
00216 
00217         double alpha_minus_epsilon = alpha - mEpsilon;
00218         double alpha_plus_epsilon = alpha + mEpsilon;
00219         double alpha_plus_pi_minus_epsilon = alpha + M_PI - mEpsilon;
00220         double alpha_plus_pi_plus_epsilon = alpha + M_PI + mEpsilon;
00221 
00222         // Calculate sampling angles in the range (-pi,pi]
00223 
00224         #define COVERAGE_IGNORE
00225         if (alpha_minus_epsilon <= -M_PI)
00226         {
00227             alpha_minus_epsilon += 2*M_PI;
00228         }
00229         #undef COVERAGE_IGNORE
00230         sampling_angles[i] = alpha_minus_epsilon;
00231 
00232         assert(sampling_angles[i] <= M_PI);
00233         assert(sampling_angles[i] > -M_PI);
00234         i++;
00235 
00236         if (alpha_plus_epsilon > M_PI)
00237         {
00238             alpha_plus_epsilon -= 2*M_PI;
00239         }
00240         sampling_angles[i] = alpha_plus_epsilon;
00241 
00242         assert(sampling_angles[i] <= M_PI);
00243         assert(sampling_angles[i] > -M_PI);
00244         i++;
00245 
00246         if (alpha_plus_pi_minus_epsilon > M_PI)
00247         {
00248             alpha_plus_pi_minus_epsilon -= 2*M_PI;
00249         }
00250         sampling_angles[i] = alpha_plus_pi_minus_epsilon;
00251 
00252         assert(sampling_angles[i] <= M_PI);
00253         assert(sampling_angles[i] > -M_PI);
00254         i++;
00255 
00256         if (alpha_plus_pi_plus_epsilon > M_PI)
00257         {
00258             alpha_plus_pi_plus_epsilon -= 2*M_PI;
00259         }
00260         sampling_angles[i] = alpha_plus_pi_plus_epsilon;
00261 
00262         assert(sampling_angles[i] <= M_PI);
00263         assert(sampling_angles[i] > -M_PI);
00264         i++;
00265     }
00266 
00267     sort(sampling_angles.begin(), sampling_angles.end());
00268     return sampling_angles;
00269 }
00270 
00271 
00272 double DiscreteSystemForceCalculator::GetLocalExtremum(unsigned index, double angle1, double angle2)
00273 {
00274     // We always pass in angle1 and angle2 such that angle1<angle2,
00275     // but note that angle1 may be <M_PI
00276     assert(angle1 < angle2);
00277 
00278     double tolerance = 1e-5;
00279     unsigned counter = 0;
00280 
00281     double previous_angle;
00282     double current_error;
00283     double current_angle = angle1;
00284 
00285     current_error = angle2 - angle1;
00286     std::vector<double> current_ft_and_fn(2);
00287 
00288     while (current_error > tolerance)
00289     {
00290         previous_angle = current_angle;
00291         current_ft_and_fn = CalculateFtAndFn(index, current_angle);
00292         current_angle -= current_ft_and_fn[0]/current_ft_and_fn[1];
00293         current_error = fabs(current_angle - previous_angle);
00294         counter++;
00295     }
00296 
00297     assert(current_angle>angle1 && current_angle<angle2);
00298     assert(current_error < tolerance);
00299 
00300     return current_angle;
00301 }
00302 
00303 
00304 std::vector<double> DiscreteSystemForceCalculator::GetExtremalAngles(unsigned index, std::vector<double> samplingAngles)
00305 {
00306     std::vector<double> extremal_angles;
00307     std::vector<double> ft_and_fn(2);
00308     std::vector<double> tangential_force(samplingAngles.size());
00309 
00310     for (unsigned i=0; i<samplingAngles.size(); i++)
00311     {
00312         ft_and_fn = CalculateFtAndFn(index, samplingAngles[i]);
00313         tangential_force[i] = ft_and_fn[0];
00314     }
00315 
00316     unsigned n = samplingAngles.size()-1;
00317 
00318     for (unsigned i=0; i<n; i++)
00319     {
00320         if ( ( tangential_force[i%n]>0 && tangential_force[(i+1)%n]<0 ) ||
00321              ( tangential_force[i%n]<0 && tangential_force[(i+1)%n]>0 ) )
00322         {
00323             double next_extremal_angle;
00324 
00325             // If we are in the interval that crosses the branch line at pi,
00326             // then subtract 2*pi from the positive angle
00327             if (i==n-1)
00328             {
00329                 samplingAngles[i%n] -= 2*M_PI;
00330             }
00331 
00332             if (samplingAngles[(i+1)%n] - samplingAngles[i%n] < 2*mEpsilon + 1e-6 )
00333             {
00334                 // If we find a jump through zero, then the local extremum is
00335                 // simply at the mid-point of the interval
00336                 next_extremal_angle = 0.5*(samplingAngles[(i+1)%n] + samplingAngles[i%n]);
00337             }
00338             else
00339             {
00340                 // Otherwise we need to find it using Newton's method
00341                 next_extremal_angle = GetLocalExtremum(index, samplingAngles[i%n], samplingAngles[(i+1)%n]);
00342             }
00343 
00344             if (next_extremal_angle <= -M_PI)
00345             {
00346                 next_extremal_angle += 2*M_PI;
00347             }
00348             assert(next_extremal_angle>-M_PI && next_extremal_angle<=M_PI);
00349             extremal_angles.push_back(next_extremal_angle);
00350         }
00351     }
00352 
00353     return extremal_angles;
00354 }

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