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 <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 }