OpenLB 1.7
Loading...
Searching...
No Matches
particleUtilities.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2021-2023 Nicolas Hafen, Jan E. Marquardt, Mathias J. Krause
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24
25
26#ifndef PARTICLE_UTILITIES_H
27#define PARTICLE_UTILITIES_H
28
29#include <cassert>
30#include <unordered_set>
31
33
34namespace olb {
35
36namespace particles {
37
38
39
40
41//Get points on spherical hull
42//TODO: has to be reworked thoroughly
43// - separate position
44// - pre-check necessary check by bounding box first
45// - ordered by nested looping: x,y,z
46template<typename T>
47std::array<Vector<T,3>,26> discretePointsOnSphericalHull(
48 Vector<T,3> position, T radius )
49{
50 //TODO: check correctness
51 T distA = radius;
52 T distB = (1/std::sqrt(2))*distA;
53 T distC = (1/std::sqrt(2))*distB;
54
55 //Point on spherical hull
56 std::array<Vector<T,3>,26> pointsOnSphericalHull = {
57 Vector<T,3>(position[0]-distC, position[1]-distC, position[2]-distC),
58 Vector<T,3>(position[0], position[1]-distB, position[2]-distB),
59 Vector<T,3>(position[0]+distC, position[1]-distC, position[2]-distC),
60
61 Vector<T,3>(position[0]-distB, position[1], position[2]-distB),
62 Vector<T,3>(position[0], position[1], position[2]-distA), //bottom
63 Vector<T,3>(position[0]+distB, position[1], position[2]-distB),
64
65 Vector<T,3>(position[0]-distC, position[1]+distC, position[2]-distC),
66 Vector<T,3>(position[0], position[1]+distB, position[2]-distB),
67 Vector<T,3>(position[0]+distC, position[1]+distC, position[2]-distC),
68
69
70 Vector<T,3>(position[0]-distB, position[1]-distB, position[2] ),
71 Vector<T,3>(position[0], position[1]-distA, position[2] ), //right
72 Vector<T,3>(position[0]+distB, position[1]-distB, position[2] ),
73
74 Vector<T,3>(position[0]-distA, position[1], position[2] ), //back
75 // Centre excluded here!
76 Vector<T,3>(position[0]+distA, position[1], position[2] ), //front
77
78 Vector<T,3>(position[0]-distB, position[1]+distB, position[2] ),
79 Vector<T,3>(position[0], position[1]+distA, position[2] ), //left
80 Vector<T,3>(position[0]+distB, position[1]+distB, position[2] ),
81
82
83 Vector<T,3>(position[0]-distC, position[1]-distC, position[2]+distC),
84 Vector<T,3>(position[0], position[1]-distB, position[2]+distB),
85 Vector<T,3>(position[0]+distC, position[1]-distC, position[2]+distC),
86
87 Vector<T,3>(position[0]-distB, position[1], position[2]+distB),
88 Vector<T,3>(position[0], position[1], position[2]+distA), //top
89 Vector<T,3>(position[0]+distB, position[1], position[2]+distB),
90
91 Vector<T,3>(position[0]-distC, position[1]+distC, position[2]+distC),
92 Vector<T,3>(position[0], position[1]+distB, position[2]+distB),
93 Vector<T,3>(position[0]+distC, position[1]+distC, position[2]+distC)
94 };
95
96 return pointsOnSphericalHull;
97}
98
99template<typename T>
100std::array<Vector<T,2>,8> discretePointsOnSphericalHull(
101 Vector<T,2> position, T radius )
102{
103 //TODO: check correctness
104 T distA = radius;
105 T distB = (1/std::sqrt(2))*distA;
106
107 //Point on spherical hull
108 std::array<Vector<T,2>,8> pointsOnSphericalHull = {
109 Vector<T,2>(position[0]-distB, position[1]-distB),
110 Vector<T,2>(position[0], position[1]-distA), //right
111 Vector<T,2>(position[0]+distB, position[1]-distB),
112
113 Vector<T,2>(position[0]-distA, position[1]), //back
114 // Centre excluded here!
115 Vector<T,2>(position[0]+distA, position[1]), //front
116
117 Vector<T,2>(position[0]-distB, position[1]+distB),
118 Vector<T,2>(position[0], position[1]+distA), //left
119 Vector<T,2>(position[0]+distB, position[1]+distB)
120 };
121
122 return pointsOnSphericalHull;
123}
124
126 template<typename T, unsigned D>
127 static constexpr auto calculate( Vector<T,D> position, T radius )
128 {
129 if constexpr(D==3 || D==2){
130 return discretePointsOnSphericalHull( position, radius );
131 } else {
132 std::cerr << "ERROR: Only 2D and 3D supported!" << std::endl;
133 return false;
134 }
135 }
136 template<typename T, typename PARTICLETYPE>
137 static constexpr auto calculate( Particle<T,PARTICLETYPE>& particle )
138 {
139 auto position = access::getPosition( particle );
140 T radius = access::getRadius( particle );
141 return calculate( position, radius );
142 }
143};
144
145template<typename T, typename PARTICLETYPE>
147 XParticleSystem<T,PARTICLETYPE>& xParticleSystem )
148{
149 using PCONDITION = conditions::invalid_particles;
150 if constexpr (access::providesParallelization<PARTICLETYPE>()){
151 auto& sParticleSystem = xParticleSystem;
152 //Iterate over particle systems
154 [&](ParticleSystem<T,PARTICLETYPE>& particleSystem, int iC, int globiC){
155 //Iterate over particles
156 std::size_t iP=0;
157 while(iP < particleSystem.size()){
158 auto particle = particleSystem.get(iP);
159 ++iP;
160 //Execute F when particle meets condition
161 doWhenMeetingCondition<T,PARTICLETYPE,PCONDITION>( particle,
162 [&](Particle<T,PARTICLETYPE> particle){
163 //Erase particle
164 auto& container = particleSystem.get();
165 container.erase( particle.getId() );
166 --iP; //reset counter
167 },globiC);
168 }
169 });
170 } else {
171 auto& particleSystem = xParticleSystem;
172 std::size_t iP=0;
173 while(iP < particleSystem.size()){
174 auto particle = particleSystem.get(iP);
175 ++iP;
176 //Execute F when particle meets condition
177 doWhenMeetingCondition<T,PARTICLETYPE,PCONDITION>( particle,
178 [&](Particle<T,PARTICLETYPE> particle){
179 //Erase particle
180 auto& container = particleSystem.get();
181 container.erase( particle.getId() );
182 --iP; //reset counter
183 });
184 }
185 }
186}
187
188template<typename T, typename PARTICLETYPE, std::size_t selectedID, typename F>
190 const bool hasFieldValid = access::providesValid<PARTICLETYPE>();
191 using PCOND = std::conditional_t<
192 hasFieldValid,
193 conditions::template valid_particle_matching_ID<selectedID>,
194 conditions::template particle_matching_ID<selectedID>
195 >;
196 //Loop over particle and find the one matching iD
197 forParticlesInXParticleSystem<T,PARTICLETYPE,PCOND>( xParticleSystem, f );
198}
199
200//Search particle in ParticleSystem by globalParticleID and return localParticleID if found
201template<typename T, typename PARTICLETYPE, typename PCONDITION=conditions::valid_particles>
203 std::size_t globalIDrequested, std::size_t& localParticleID )
204{
205 using namespace descriptors;
206 static_assert(PARTICLETYPE::template providesNested<PARALLELIZATION,ID>(), "Field PARALLELIZATION:ID has to be provided");
207 //Initialize return quantities
208 localParticleID = 0;
209 bool found = false;
210 //Loop over particles
211 for (std::size_t iP=0; iP<particleSystem.size(); ++iP) {
212 if (!found){
213 auto particle = particleSystem.get(iP);
214 //Execute F when particle meets condition
215 doWhenMeetingCondition<T,PARTICLETYPE,PCONDITION>( particle,
216 [&](Particle<T,PARTICLETYPE> particle){
217 auto globalID = particle.template getField<PARALLELIZATION,ID>();
218 if (globalID==globalIDrequested){
219 localParticleID = particle.getId();
220 found=true;
221 }
222 });
223 } else {
224 break;
225 }
226 }
227 return found;
228}
229
230
231
232namespace communication {
233
234
235
239template<typename T, typename PARTICLETYPE, typename PCONDITION=conditions::valid_particle_centres, bool sync=true>
241 std::size_t globalIDrequested, std::size_t& localParticleID, int& globiC )
242{
243 using namespace descriptors;
244 static_assert(PARTICLETYPE::template providesNested<PARALLELIZATION,ID>(), "Field PARALLELIZATION:ID has to be provided");
245 //Retrieve load balancer
246 auto& loadBalancer = sParticleSystem.getSuperStructure().getLoadBalancer();
247 //Initialize return quantities
248 localParticleID = 0;
249 globiC = 0;
250 bool found = false;
251 //Loop over particle systems
252 for (int iC=0; iC<loadBalancer.size(); ++iC){
253 if (!found){
254 int globiC = loadBalancer.glob(iC);
255 //Retrieve container
256 auto bParticleSystems = sParticleSystem.getBlockParticleSystems();
257 auto& particleSystem = *bParticleSystems[iC];
258 //Search particle on particleSystem
259 //NOTE: Can not be substituted with searchParticle( ParticleSystem ..) due to globiC
260 for (std::size_t iP=0; iP<particleSystem.size(); ++iP) {
261 if (!found){
262 auto particle = particleSystem.get(iP);
263 //Execute F when particle meets condition
264 doWhenMeetingCondition<T,PARTICLETYPE,PCONDITION>( particle,
265 [&](Particle<T,PARTICLETYPE> particle){
266 auto globalID = particle.template getField<PARALLELIZATION,ID>();
267 if (globalID==globalIDrequested){
268 localParticleID = particle.getId();
269 found=true;
270 }
271 }, globiC);
272 } else {
273 break;
274 }
275 } //for (std::size_t iP=0; iP<particleSystem.size(); ++iP)
276 } else {
277 break;
278 }
279 } //for (int iC=0; iC<loadBalancer.size(); ++iC)
280#ifdef PARALLEL_MODE_MPI
281 if constexpr(sync){
282 singleton::mpi().reduceAndBcast(found, MPI_LOR);
283 if (found){
284 singleton::mpi().reduceAndBcast(globiC, MPI_MAX);
285 singleton::mpi().reduceAndBcast(localParticleID, MPI_MAX);
286 }
287 }
288#endif
289 return found;
290}
291
292
297template<typename T, unsigned D, typename F=std::function<void(int)>,
298 bool domainWarning=false, bool checkDiscretePoints=false>
299std::unordered_set<int> getSurfaceTouchingICs(
300 CuboidGeometry<T,D>& cuboidGeometry,
301 Vector<T,D> position, T circumRadius,
302 const Vector<bool, D>& periodicity,
303 int globiC, F f=[](int){} )
304{
305 // This option is prone to errors
306 if constexpr(checkDiscretePoints) {
307 //Set of destination iCs preventing sending to same iC multiple times
308 std::unordered_set<int> iCs;
309 //Retrieve points on hull
310 auto pointsOnHull = discrete_points_on_hull::calculate( position, circumRadius );
311 //Iterate over points on hull
312 for (const PhysR<T,D>& point : pointsOnHull){
313 //Retrieve residence iC
314 int globiConHull=0;
315 const bool cuboidFound = getCuboid(cuboidGeometry,
316 periodicity, point, globiConHull);
317 if constexpr(domainWarning){
318 if (!cuboidFound){
319 std::cerr << "Warning [during getSurfaceTouchingICs]: Cuboid not found for point "
320 << point << std::endl;
321 }
322 }
323 //Check whether point has left responsible ic
324 if (cuboidFound && globiConHull!=globiC){
325 //Add destination ic to list (to avoid resending)
326 if(iCs.insert(globiConHull).second) {
327 f(globiConHull);
328 }
329 } //if (globiConHull!=globiC)
330 } //for (auto point : pointsOnHull)
331 return iCs;
332 }
333 else {
334 const T physDeltaX = cuboidGeometry.getMotherCuboid().getDeltaR();
335 std::vector<int> tmp = communication::getNeighborCuboids<T,D,false>(
336 cuboidGeometry, util::ceil(circumRadius/physDeltaX), globiC);
337 std::unordered_set<int> iCs;
338 for(const int entry : tmp) {
339 if( globiC != entry ) {
340 iCs.insert(entry);
341 f(entry);
342 }
343 }
344 return iCs;
345 }
346 __builtin_unreachable();
347}
348
349template<typename T, typename PARTICLETYPE, typename F=std::function<void(int)>>
350std::unordered_set<int> getSurfaceTouchingICs(
352 Vector<T,PARTICLETYPE::d> position, T circumRadius,
353 const Vector<bool, PARTICLETYPE::d>& periodicity,
354 int globiC, F f=[](int){} )
355{
356 std::unordered_set<int> neighbors = sParticleSystem.getCuboidNeighborhood()[globiC];
357 for(const int neighbor: neighbors) {
358 f(neighbor);
359 }
360 return neighbors;
361}
362
363
364//Get list of surface touching iCs (that are not globiC)
365template<typename T, unsigned D, typename F=std::function<void(int)>,
366 bool domainWarning=false>
368 CuboidGeometry<T,D>& cuboidGeometry,
369 Vector<T,D> position, T circumRadius,
370 const Vector<bool, D>& periodicity,
371 int globiC, F f=[](int){} )
372{
373 const std::unordered_set<int> touchedICs{getSurfaceTouchingICs<T,D,F,domainWarning>(
374 cuboidGeometry, position, circumRadius, periodicity, globiC, f)};
375 std::vector<int> listOfICs;
376 listOfICs.insert(listOfICs.end(), touchedICs.begin(), touchedICs.end());
377 return listOfICs;
378}
379
380
381} //namespace communication
382
383
384namespace sorting {
385
386
387//Sort particleSystem by provided nested fields
388//TODO: include derived field function
389//TODO: for now only boolean, numeric comparison should be provided as well
390//TODO: possibly generalize to work with all containers or objects, such as blockLattice
391// - this should be possible at least, when particleSystem and Lattice have the same
392// structure/base type
393template<typename T, typename PARTICLETYPE, typename ...NESTED_FIELDS>
395 using namespace descriptors;
396 //Find first occurence of element belonging to partition B
397 std::size_t iPfirstB=particleSystem.size(); //Fallback (simple to check)
398 bool valB = true;
399 for (std::size_t iP=0; iP<particleSystem.size(); ++iP){
400 auto particle = particleSystem.get(iP);
401 bool active = particle.template getField<NESTED_FIELDS...>();
402 if (active==valB){
403 iPfirstB=iP;
404 break;
405 }
406 }
407 //Find succeeding elements belonging to partition A and move them to othe beginning
408 for (std::size_t iP=iPfirstB; iP<particleSystem.size(); ++iP){
409 auto particle = particleSystem.get(iP);
410 bool active = particle.template getField<NESTED_FIELDS...>();
411 if (active!=valB){
412 particleSystem.swapParticles(iP,iPfirstB);
413 ++iPfirstB;
414 }
415 }
416 //Return splitpoint between partitions
417 return iPfirstB;
418}
419
420} //namespace sorting
421
422} //namespace particles
423
424} //namespace olb
425
426
427#endif
LoadBalancer< T > & getLoadBalancer()
Read and write access to the load balancer.
Plain old scalar vector.
Definition vector.h:47
auto & get()
Expose container.
void swapParticles(std::size_t iP, std::size_t jP)
Swap particles by index.
constexpr std::size_t size()
Size of ParticleSystem.
std::size_t getId() const
Return memory ID of the currently represented particle.
Definition particle.hh:67
SuperStructure< T, PARTICLETYPE::d > & getSuperStructure()
const std::vector< std::unordered_set< int > > & getCuboidNeighborhood()
std::vector< ParticleSystem< T, PARTICLETYPE > * > & getBlockParticleSystems()
T getRadius(Particle< T, PARTICLETYPE > &particle)
Vector< T, PARTICLETYPE::d > getPosition(Particle< T, PARTICLETYPE > particle)
void forSystemsInSuperParticleSystem(SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, F f)
std::unordered_set< int > getSurfaceTouchingICs(CuboidGeometry< T, D > &cuboidGeometry, Vector< T, D > position, T circumRadius, const Vector< bool, D > &periodicity, int globiC, F f=[](int){})
Get a set of surface touching iCs (that are not globiC) Allows to run an optional function per unique...
std::vector< int > getVectorOfSurfaceTouchingICs(CuboidGeometry< T, D > &cuboidGeometry, Vector< T, D > position, T circumRadius, const Vector< bool, D > &periodicity, int globiC, F f=[](int){})
bool searchParticleGlobally(SuperParticleSystem< T, PARTICLETYPE > &sParticleSystem, std::size_t globalIDrequested, std::size_t &localParticleID, int &globiC)
Search particle in SuperParticleSystem by globalParticleID and return globiC if found.
std::size_t partitionParticleSystem(ParticleSystem< T, PARTICLETYPE > &particleSystem)
std::conditional_t< PARTICLETYPE::template providesNested< descriptors::PARALLELIZATION >(), SuperParticleSystem< T, PARTICLETYPE >, ParticleSystem< T, PARTICLETYPE > > XParticleSystem
void purgeInvalidParticles(XParticleSystem< T, PARTICLETYPE > &xParticleSystem)
std::array< Vector< T, 3 >, 26 > discretePointsOnSphericalHull(Vector< T, 3 > position, T radius)
bool searchParticleLocally(ParticleSystem< T, PARTICLETYPE > &particleSystem, std::size_t globalIDrequested, std::size_t &localParticleID)
void doForParticleMatchingID(XParticleSystem< T, PARTICLETYPE > &xParticleSystem, F f)
Top level namespace for all of OpenLB.
std::conditional_t< D==2, CuboidGeometry2D< T >, CuboidGeometry3D< T > > CuboidGeometry
Definition aliases.h:47
static constexpr auto calculate(Particle< T, PARTICLETYPE > &particle)
static constexpr auto calculate(Vector< T, D > position, T radius)