ESyS-Particle  4.0.1
DistConnections.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 <limits>
00016 
00017 namespace esys
00018 {
00019   namespace lsm
00020   {
00021     template <typename TmplParticle, typename TmplConnection>
00022     DistConnections<TmplParticle,TmplConnection>::DistConnections(
00023       double maxDist,
00024       Tag defaultTag,
00025       const BoundingBox &bBox,
00026       const BoolVector &circDimensions
00027     )
00028       : m_connectionPoolPtr(new ConnectionPool(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         m_defaultTag(defaultTag)
00037     {
00038       const double gridSize =
00039         (
00040           max(
00041             bBox.getSizes()[0],
00042             max(
00043               bBox.getSizes()[1],
00044               bBox.getSizes()[2]
00045             )
00046           )
00047         )/5.0;
00048       m_nTablePtr =
00049         NTablePtr(
00050           new NTable(
00051             bBox,
00052             gridSize,
00053             circDimensions,
00054             2*gridSize
00055           )
00056       );
00057     }
00058 
00059     template <typename TmplParticle, typename TmplConnection>
00060     DistConnections<TmplParticle,TmplConnection>::~DistConnections()
00061     {
00062     }
00063 
00064     template <typename TmplParticle, typename TmplConnection>
00065     int DistConnections<TmplParticle,TmplConnection>::getNumParticles() const
00066     {
00067       return m_nTablePtr->size();
00068     }
00069 
00070     template <typename TmplParticle, typename TmplConnection>
00071     int DistConnections<TmplParticle,TmplConnection>::getNumConnections() const
00072     {
00073       return m_connectionSet.size();
00074     }
00075 
00076     template <typename TmplParticle, typename TmplConnection>
00077     double DistConnections<TmplParticle,TmplConnection>::getMinRadius() const
00078     {
00079       return m_minRadius;
00080     }
00081 
00082     template <typename TmplParticle, typename TmplConnection>
00083     double DistConnections<TmplParticle,TmplConnection>::getMaxRadius() const
00084     {
00085       return m_maxRadius;
00086     }
00087 
00088     template <typename TmplParticle, typename TmplConnection>
00089     typename DistConnections<TmplParticle,TmplConnection>::ParticleConstIterator
00090     DistConnections<TmplParticle,TmplConnection>::getParticleIterator() const
00091     {
00092       return m_nTablePtr->getIterator();
00093     }
00094 
00095     template <typename TmplParticle, typename TmplConnection>
00096     void
00097     DistConnections<TmplParticle,TmplConnection>::createConnection(
00098       const Particle &p1,
00099       const Particle &p2,
00100       Tag tag
00101     )
00102     {
00103       m_connectionSet.insert(
00104         m_connectionPoolPtr->construct(p1.getId(), p2.getId(), tag)
00105       );
00106     }
00107     template <typename TmplParticle>
00108     class CmpParticleId
00109     {
00110     public:
00111       bool operator()(const TmplParticle &p1, const TmplParticle &p2) const
00112       {
00113         return (p1.getId() < p2.getId());
00114       }
00115 
00116       bool operator()(const TmplParticle *p1, const TmplParticle *p2) const
00117       {
00118         return (p1->getId() < p2->getId());
00119       }
00120     };
00121 
00122     template <typename TmplParticle, typename TmplConnection>
00123     template <typename TmplParticleIterator>
00124     void
00125     DistConnections<TmplParticle,TmplConnection>::create(
00126       TmplParticleIterator it,
00127       Tag tag
00128     )
00129     {
00130       typedef std::set<Particle *, CmpParticleId<Particle> > ParticleSet;
00131       ParticleSet pSet;
00132       while (it.hasNext())
00133       {
00134         Particle &p = it.next();
00135         insert(p);
00136         pSet.insert(&p);
00137       }
00138       m_nTablePtr->resize(getParticleBBox(), 4.1*getMinRadius(), 2.1*getMaxRadius());
00139 
00140       for (
00141         typename ParticleSet::const_iterator pIt = pSet.begin();
00142         pIt != pSet.end();
00143         pIt++
00144       )
00145       {
00146         typename NTable::ParticleVector nVector =
00147           m_nTablePtr->getNeighbourVector(
00148             (*pIt)->getPos(),
00149             (*pIt)->getRad() + m_maxDist
00150           );
00151         for (
00152           typename NTable::ParticleVector::const_iterator nIt = nVector.begin();
00153           nIt != nVector.end();
00154           nIt++
00155         )
00156         {
00157           Particle *p1 = (*pIt);
00158           Particle *p2 = (*nIt);
00159 
00160           if (
00161             (
00162               (pSet.find(p1) != pSet.end())
00163               &&
00164               (pSet.find(p2) != pSet.end())
00165               &&
00166               (p1->getId() < p2->getId())
00167             )
00168             ||
00169             (
00170               ((pSet.find(p1)==pSet.end()) && (pSet.find(p2)!= pSet.end()))
00171               ||
00172               ((pSet.find(p1)!=pSet.end()) && (pSet.find(p2)== pSet.end()))
00173             )
00174           )
00175           {
00176             p1 =
00177               ((*pIt)->getId() < (*nIt)->getId())
00178               ?
00179               (*pIt)
00180               :
00181               (*nIt);
00182             p2 =
00183               ((*pIt)->getId() < (*nIt)->getId())
00184               ?
00185               (*nIt)
00186               :
00187               (*pIt);
00188             const double radiusSum =
00189               m_maxDist + p1->getRad() + p2->getRad();
00190             const double radiusSumSqrd = radiusSum*radiusSum;
00191   
00192             if (
00193               (p1->getPos() - p2->getPos()).norm2()
00194               <=
00195               (radiusSumSqrd)
00196             )
00197             {
00198 #if 0
00199               console.Debug()
00200                 << "creating connection: \n"
00201                 << StringUtil::toString(*p1)
00202                 << "->"
00203                 << StringUtil::toString(*p2) << "\n";
00204 #endif
00205               createConnection(*p1, *p2, tag);
00206             }
00207           }
00208         }
00209       }
00210     }
00211 
00212     template <typename TmplParticle, typename TmplConnection>
00213     template <typename TmplParticleIterator>
00214     void
00215     DistConnections<TmplParticle,TmplConnection>::create(
00216       TmplParticleIterator it
00217     )
00218     {
00219       create(it, getDefaultTag());
00220     }
00221     
00222     template <typename TmplParticle, typename TmplConnection>
00223     void
00224     DistConnections<TmplParticle,TmplConnection>::insert(Particle &p)
00225     {
00226       if (p.getRad() < m_minRadius)
00227       {
00228         m_minRadius = p.getRad();
00229       }
00230       if (p.getRad() > m_maxRadius)
00231       {
00232         m_maxRadius = p.getRad();
00233       }
00234 
00235       m_nTablePtr->insert(p);
00236 
00237       for (int i = 0; i < 3; i++)
00238       {
00239         if (!(m_nTablePtr->getPeriodicDimensions()[i]))
00240         {
00241           if (p.getPos()[i]-p.getRad() < m_minPt[i])
00242           {
00243             m_minPt[i] = p.getPos()[i]-p.getRad();
00244           }
00245           if (p.getPos()[i]+p.getRad() > m_maxPt[i])
00246           {
00247             m_maxPt[i] = p.getPos()[i]+p.getRad();
00248           }
00249         }
00250       }
00251     }
00252 
00253     template <typename TmplParticle, typename TmplConnection>
00254     typename DistConnections<TmplParticle,TmplConnection>::Tag
00255     DistConnections<TmplParticle,TmplConnection>::getDefaultTag() const
00256     {
00257       return m_defaultTag;
00258     }
00259 
00260     template <typename TmplParticle, typename TmplConnection>
00261     void
00262     DistConnections<TmplParticle,TmplConnection>::setDefaultTag(Tag defaultTag)
00263     {
00264       m_defaultTag = defaultTag;
00265     }
00266 
00267     template <typename TmplParticle, typename TmplConnection>
00268     BoundingBox
00269     DistConnections<TmplParticle,TmplConnection>::getParticleBBox() const
00270     {
00271       return BoundingBox(m_minPt, m_maxPt);
00272     }
00273 
00274   }
00275 }