ESyS-Particle
4.0.1
|
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 }