AbstractTwoBodyInteractionForce.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 "AbstractTwoBodyInteractionForce.hpp"
00030 
00031 
00032 template<unsigned DIM>
00033 AbstractTwoBodyInteractionForce<DIM>::AbstractTwoBodyInteractionForce()
00034    : AbstractForce<DIM>(),
00035      mUseCutOffLength(false),
00036      mMechanicsCutOffLength(DBL_MAX)
00037 {
00038 }
00039 
00040 template<unsigned DIM>
00041 bool AbstractTwoBodyInteractionForce<DIM>::GetUseCutOffLength()
00042 {
00043     return mUseCutOffLength;
00044 }
00045 
00046 template<unsigned DIM>
00047 void AbstractTwoBodyInteractionForce<DIM>::SetCutOffLength(double cutOffLength)
00048 {
00049     assert(cutOffLength > 0.0);
00050     mUseCutOffLength = true;
00051     mMechanicsCutOffLength = cutOffLength;
00052 }
00053 
00054 
00055 template<unsigned DIM>
00056 double AbstractTwoBodyInteractionForce<DIM>::GetCutOffLength()
00057 {
00058     return mMechanicsCutOffLength;
00059 }
00060 
00061 
00062 template<unsigned DIM>
00063 void AbstractTwoBodyInteractionForce<DIM>::AddForceContribution(std::vector<c_vector<double, DIM> >& rForces,
00064                                                                 AbstractCellPopulation<DIM>& rCellPopulation)
00065 {
00066     // Throw an exception message if not using a subclass of AbstractCentreBasedCellPopulation
00067     if (dynamic_cast<AbstractCentreBasedCellPopulation<DIM>*>(&rCellPopulation) == NULL)
00068     {
00069         EXCEPTION("Subclasses of AbstractTwoBodyInteractionForce are to be used with subclasses of AbstractCentreBasedCellPopulation only");
00070     }
00071 
00072     if (rCellPopulation.IsMeshBasedCellPopulation())
00073     {
00074         MeshBasedCellPopulation<DIM>* p_static_cast_cell_population = static_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation);
00075 
00076         // Iterate over all springs and add force contributions
00077         for (typename MeshBasedCellPopulation<DIM>::SpringIterator spring_iterator = p_static_cast_cell_population->SpringsBegin();
00078              spring_iterator != p_static_cast_cell_population->SpringsEnd();
00079              ++spring_iterator)
00080         {
00081             unsigned nodeA_global_index = spring_iterator.GetNodeA()->GetIndex();
00082             unsigned nodeB_global_index = spring_iterator.GetNodeB()->GetIndex();
00083 
00084             // Calculate the force between nodes
00085             c_vector<double, DIM> force = CalculateForceBetweenNodes(nodeA_global_index, nodeB_global_index, rCellPopulation);
00086 
00087             // Add the force contribution to each node
00088             rForces[nodeB_global_index] -= force;
00089             rForces[nodeA_global_index] += force;
00090         }
00091     }
00092     else
00093     {
00094         std::set< std::pair<Node<DIM>*, Node<DIM>* > >& r_node_pairs = (static_cast<NodeBasedCellPopulation<DIM>*>(&rCellPopulation))->rGetNodePairs();
00095 
00096 //        assert(DIM==2); // 3d boxes not implemented yet - if fails nightly uncomment the double node loop below
00097                         // and use that for the 3d case
00098         for (typename std::set< std::pair<Node<DIM>*, Node<DIM>* > >::iterator iter = r_node_pairs.begin();
00099             iter != r_node_pairs.end();
00100             iter++)
00101         {
00102             std::pair<Node<DIM>*, Node<DIM>* > pair = *iter;
00103 
00104             unsigned node_a_index = pair.first->GetIndex();
00105             unsigned node_b_index = pair.second->GetIndex();
00106 
00107             // Calculate the force between nodes
00108             c_vector<double, DIM> force = CalculateForceBetweenNodes(node_a_index, node_b_index, rCellPopulation);
00109             for (unsigned j=0; j<DIM; j++)
00110             {
00111                 assert(!std::isnan(force[j]));
00112             }
00113 
00114             // Add the force contribution to each node
00115             rForces[node_a_index] += force;
00116             rForces[node_b_index] -= force;
00117         }
00118 
00119 //        // Iterate over nodes
00120 //        for (unsigned node_a_index=0; node_a_index<rCellPopulation.GetNumNodes(); node_a_index++)
00121 //        {
00122 //            // Iterate over nodes
00123 //            for (unsigned node_b_index=node_a_index+1; node_b_index<rCellPopulation.GetNumNodes(); node_b_index++)
00124 //            {
00125 //                // Calculate the force between nodes
00126 //                c_vector<double, DIM> force = CalculateForceBetweenNodes(node_a_index, node_b_index, rCellPopulation);
00127 //                for (unsigned j=0; j<DIM; j++)
00128 //                {
00129 //                    assert(!std::isnan(force[j]));
00130 //                }
00131 //
00132 //                // Add the force contribution to each node
00133 //                rForces[node_a_index] += force;
00134 //                rForces[node_b_index] -= force;
00135 //            }
00136 //        }
00137     }
00138 }
00139 
00140 template<unsigned DIM>
00141 void AbstractTwoBodyInteractionForce<DIM>::OutputForceParameters(out_stream& rParamsFile)
00142 {
00143     *rParamsFile << "\t\t\t<UseCutOffLength>" << mUseCutOffLength << "</UseCutOffLength> \n";
00144     *rParamsFile << "\t\t\t<CutOffLength>" << mMechanicsCutOffLength << "</CutOffLength> \n";
00145 
00146     // Call method on direct parent class
00147     AbstractForce<DIM>::OutputForceParameters(rParamsFile);
00148 }
00149 
00151 // Explicit instantiation
00153 
00154 template class AbstractTwoBodyInteractionForce<1>;
00155 template class AbstractTwoBodyInteractionForce<2>;
00156 template class AbstractTwoBodyInteractionForce<3>;

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