ESyS-Particle  4.0.1
SphereNeighbours.hpp
00001 
00002 //                                                         //
00003 // Copyright (c) 2003-2011 by The University of Queensland //
00004 // Earth Systems Science Computational Centre (ESSCC)      //
00005 // http://www.uq.edu.au/esscc                              //
00006 //                                                         //
00007 // Primary Business: Brisbane, Queensland, Australia       //
00008 // Licensed under the Open Software License version 3.0    //
00009 // http://www.opensource.org/licenses/osl-3.0.php          //
00010 //                                                         //
00012 
00013 
00014 #include "Foundation/StringUtil.h"
00015 #include <functional>
00016 #include <limits>
00017 
00018 namespace esys
00019 {
00020   namespace lsm
00021   {
00022     template <typename TmplSphere, typename TmplIdPair>
00023     SphereNeighbours<TmplSphere,TmplIdPair>::SphereNeighbours(
00024       double maxDist,
00025       const BoundingBox &bBox,
00026       const BoolVector &circDimensions
00027     )
00028       : m_connectionPoolPtr(new IdPairPool(4096)),
00029         m_connectionSet(),
00030         m_nTablePtr(),
00031         m_minRadius(std::numeric_limits<double>::max()),
00032         m_maxRadius(-std::numeric_limits<double>::max()),
00033         m_maxDist(maxDist),
00034         m_minPt(bBox.getMinPt()),
00035         m_maxPt(bBox.getMaxPt())
00036     {
00037       const double gridSize =
00038         (
00039           max(
00040             bBox.getSizes()[0],
00041             max(
00042               bBox.getSizes()[1],
00043               bBox.getSizes()[2]
00044             )
00045           )
00046         )/5.0;
00047       m_nTablePtr =
00048         NTablePtr(
00049           new NTable(
00050             bBox,
00051             gridSize,
00052             circDimensions,
00053             2*gridSize
00054           )
00055       );
00056     }
00057 
00058     template <typename TmplSphere, typename TmplIdPair>
00059     SphereNeighbours<TmplSphere,TmplIdPair>::~SphereNeighbours()
00060     {
00061     }
00062 
00063     template <typename TmplSphere, typename TmplIdPair>
00064     int SphereNeighbours<TmplSphere,TmplIdPair>::getNumSpheres() const
00065     {
00066       return m_nTablePtr->size();
00067     }
00068 
00069     template <typename TmplSphere, typename TmplIdPair>
00070     int SphereNeighbours<TmplSphere,TmplIdPair>::getNumIdPairs() const
00071     {
00072       return m_connectionSet.size();
00073     }
00074 
00075     template <typename TmplSphere, typename TmplIdPair>
00076     double SphereNeighbours<TmplSphere,TmplIdPair>::getMinRadius() const
00077     {
00078       return m_minRadius;
00079     }
00080 
00081     template <typename TmplSphere, typename TmplIdPair>
00082     double SphereNeighbours<TmplSphere,TmplIdPair>::getMaxRadius() const
00083     {
00084       return m_maxRadius;
00085     }
00086 
00087     template <typename TmplSphere, typename TmplIdPair>
00088     typename SphereNeighbours<TmplSphere,TmplIdPair>::SphereConstIterator
00089     SphereNeighbours<TmplSphere,TmplIdPair>::getSphereIterator() const
00090     {
00091       return m_nTablePtr->getIterator();
00092     }
00093 
00094     template <typename TmplSphere, typename TmplIdPair>
00095     const typename SphereNeighbours<TmplSphere,TmplIdPair>::IdPair &
00096     SphereNeighbours<TmplSphere,TmplIdPair>::createIdPair(
00097       const Sphere &p1,
00098       const Sphere &p2
00099     )
00100     {
00101       return 
00102         **(m_connectionSet.insert(
00103           m_connectionPoolPtr->construct(p1.getId(), p2.getId())
00104         ).first);
00105     }
00106 
00107     template <typename TmplSphere>
00108     class CmpSphereId
00109     {
00110     public:
00111       bool operator()(const TmplSphere &p1, const TmplSphere &p2) const
00112       {
00113         return (p1.getId() < p2.getId());
00114       }
00115 
00116       bool operator()(const TmplSphere *p1, const TmplSphere *p2) const
00117       {
00118         return (p1->getId() < p2->getId());
00119       }
00120     };
00121 
00122     template <typename TmplType>
00123     class Deref
00124     {
00125     public:
00126       typedef const TmplType& result_type;
00127       typedef const TmplType* argument_type;
00128   
00129       result_type
00130       operator()(argument_type a) const
00131       { return *a; }
00132     };
00133     
00134     template <typename TmplSphere, typename TmplIdPair>
00135     template <typename TmplSphereIterator>
00136     typename SphereNeighbours<TmplSphere,TmplIdPair>::IdPairVector
00137     SphereNeighbours<TmplSphere,TmplIdPair>::getNeighbours(
00138       TmplSphereIterator it
00139     )
00140     {
00141       // Put the spheres in the neighbour table
00142       typedef std::set<Sphere *, CmpSphereId<Sphere> > SphereSet;
00143       SphereSet pSet;
00144       while (it.hasNext())
00145       {
00146         Sphere &p = it.next();
00147         insert(p);
00148         pSet.insert(&p);
00149       }
00150       ConstIdPairSet idPairSet;
00151       // Resize ntable according to min and max sphere radii.
00152       m_nTablePtr->resize(
00153         getSphereBBox(),
00154         4.1*getMinRadius(),
00155         2.1*getMaxRadius()
00156       );
00157 
00158       // For each sphere in the iterator it, determine it's
00159       // neighours within m_maxDist distance.
00160       for (
00161         typename SphereSet::const_iterator pIt = pSet.begin();
00162         pIt != pSet.end();
00163         pIt++
00164       )
00165       {
00166         // get the vector of spheres which are contained in nearby cells.
00167         typename NTable::ParticleVector nVector =
00168           m_nTablePtr->getNeighbourVector(
00169             (*pIt)->getPos(),
00170             (*pIt)->getRad() + m_maxDist
00171           );
00172         // for each of these cell-spheres, determine if they
00173         // are closer than m_maxDist.
00174         for (
00175           typename NTable::ParticleVector::const_iterator nIt = nVector.begin();
00176           nIt != nVector.end();
00177           nIt++
00178         )
00179         {
00180           Sphere *p1 = (*pIt);
00181           Sphere *p2 = (*nIt);
00182 
00183           if (
00184             (
00185               (p1->getId() < p2->getId())
00186               &&
00187               (pSet.find(p1) != pSet.end())
00188               &&
00189               (pSet.find(p2) != pSet.end())
00190             )
00191             ||
00192             (
00193               ((pSet.find(p1)==pSet.end()) && (pSet.find(p2)!= pSet.end()))
00194               ||
00195               ((pSet.find(p1)!=pSet.end()) && (pSet.find(p2)== pSet.end()))
00196             )
00197           )
00198           {
00199             p1 =
00200               ((*pIt)->getId() < (*nIt)->getId())
00201               ?
00202               (*pIt)
00203               :
00204               (*nIt);
00205             p2 =
00206               ((*pIt)->getId() < (*nIt)->getId())
00207               ?
00208               (*nIt)
00209               :
00210               (*pIt);
00211             const double radiusSumPlusTol =
00212               m_maxDist + p1->getRad() + p2->getRad();
00213             const double radiusSumPlusTolSqrd =
00214               radiusSumPlusTol*radiusSumPlusTol;
00215   
00216             if (
00217               (p1->getPos() - p2->getPos()).norm2()
00218               <=
00219               (radiusSumPlusTolSqrd)
00220             )
00221             {
00222               idPairSet.insert(&createIdPair(*p1, *p2));
00223             }
00224           }
00225         }
00226       }
00227       IdPairVector idPairVector;
00228       idPairVector.reserve(idPairSet.size());
00229       std::transform(
00230         idPairSet.begin(),
00231         idPairSet.end(),
00232         std::back_insert_iterator<IdPairVector>(idPairVector),
00233         Deref<IdPair>()
00234       );
00235       return idPairVector;
00236     }
00237 
00238     template <typename TmplSphere, typename TmplIdPair>
00239     void
00240     SphereNeighbours<TmplSphere,TmplIdPair>::insert(Sphere &p)
00241     {
00242       if (p.getRad() < m_minRadius)
00243       {
00244         m_minRadius = p.getRad();
00245       }
00246       if (p.getRad() > m_maxRadius)
00247       {
00248         m_maxRadius = p.getRad();
00249       }
00250 
00251       m_nTablePtr->insert(p);
00252 
00253       for (int i = 0; i < 3; i++)
00254       {
00255         if (!(m_nTablePtr->getPeriodicDimensions()[i]))
00256         {
00257           if (p.getPos()[i]-p.getRad() < m_minPt[i])
00258           {
00259             m_minPt[i] = p.getPos()[i]-p.getRad();
00260           }
00261           if (p.getPos()[i]+p.getRad() > m_maxPt[i])
00262           {
00263             m_maxPt[i] = p.getPos()[i]+p.getRad();
00264           }
00265         }
00266       }
00267     }
00268 
00269     template <typename TmplSphere, typename TmplIdPair>
00270     BoundingBox
00271     SphereNeighbours<TmplSphere,TmplIdPair>::getSphereBBox() const
00272     {
00273       return BoundingBox(m_minPt, m_maxPt);
00274     }
00275   }
00276 }