Line |
Branch |
Exec |
Source |
1 |
|
|
#include "analyzers.hpp" |
2 |
|
|
#include "interaction.hpp" |
3 |
|
|
#include "particle.hpp" |
4 |
|
|
#include "potentials.hpp" |
5 |
|
|
#include "types.hpp" |
6 |
|
|
#include <algorithm> |
7 |
|
|
#include <iterator> |
8 |
|
|
#include <map> |
9 |
|
|
#include <numeric> |
10 |
|
|
#include <iostream> |
11 |
|
|
|
12 |
|
✗ |
Clusterer::Clusterer(std::vector<Particle> &particles, const Vec &L, double cutoff) |
13 |
|
✗ |
: particles(particles), L(L), cutoff(cutoff), cutoffsq(cutoff * cutoff) { |
14 |
|
|
} |
15 |
|
|
|
16 |
|
✗ |
void Clusterer::reset() { |
17 |
|
✗ |
int n = particles.size(); |
18 |
|
✗ |
particle2Coordination.resize(n); |
19 |
|
✗ |
std::fill(particle2Coordination.begin(), particle2Coordination.end(), 0); |
20 |
|
✗ |
particle2Cluster.resize(n); |
21 |
|
✗ |
cluster2Particles.clear(); |
22 |
|
✗ |
for (int i = 0; i < n; ++i) { |
23 |
|
✗ |
particle2Cluster[i] = i; |
24 |
|
✗ |
cluster2Particles[i] = {i}; |
25 |
|
|
} |
26 |
|
|
} |
27 |
|
|
|
28 |
|
✗ |
double Clusterer::operator()(Particle &p1, Particle &p2) { |
29 |
|
✗ |
int i = std::distance(&particles.front(), &p1); |
30 |
|
✗ |
int j = std::distance(&particles.front(), &p2); |
31 |
|
✗ |
double distsq = r_ij(L, particles[i].r, particles[j].r).squaredNorm(); |
32 |
|
✗ |
if (distsq < cutoffsq && particle2Cluster[i] != particle2Cluster[j]) { |
33 |
|
✗ |
particle2Coordination[i]++; |
34 |
|
✗ |
particle2Coordination[j]++; |
35 |
|
✗ |
int oldCluster = std::max(particle2Cluster[i], particle2Cluster[j]); |
36 |
|
✗ |
int newCluster = std::min(particle2Cluster[i], particle2Cluster[j]); |
37 |
|
✗ |
for (int i : cluster2Particles[oldCluster]) { |
38 |
|
✗ |
particle2Cluster[i] = newCluster; |
39 |
|
|
} |
40 |
|
✗ |
cluster2Particles[newCluster].splice(cluster2Particles[newCluster].end(), |
41 |
|
✗ |
cluster2Particles[oldCluster]); |
42 |
|
|
} |
43 |
|
✗ |
return 0.; |
44 |
|
|
} |
45 |
|
|
|
46 |
|
✗ |
void Clusterer::sortClusters() { |
47 |
|
|
// At this point, particles of the same cluster have the same id in particle2Cluster, which can |
48 |
|
|
// still be in the range [0, total number of particles) |
49 |
|
|
// Now change data structures such that cluster ids go from 0 to numClusters - 1 and are ordered |
50 |
|
|
// by cluster size |
51 |
|
✗ |
int n = particles.size(); |
52 |
|
✗ |
std::set<int> clusterIdSet(particle2Cluster.begin(), |
53 |
|
✗ |
particle2Cluster.end()); // get all unique ids |
54 |
|
✗ |
std::vector<int> clusterIdVector(clusterIdSet.begin(), |
55 |
|
✗ |
clusterIdSet.end()); // transform to vector for accessible index |
56 |
|
✗ |
int numClusters = clusterIdVector.size(); |
57 |
|
|
// Construct map oldId -> newId where oldId in [0, numParticles) and newId in [0, numClusters) |
58 |
|
✗ |
std::map<int, int> clusterIdMap; |
59 |
|
✗ |
for (int newId = 0; newId < numClusters; ++newId) { |
60 |
|
✗ |
int oldId = clusterIdVector[newId]; |
61 |
|
✗ |
clusterIdMap[oldId] = newId; |
62 |
|
|
} |
63 |
|
|
// Change cluster ids of data structures to new restricted range |
64 |
|
✗ |
for (int &clusterId : particle2Cluster) { |
65 |
|
✗ |
clusterId = clusterIdMap[clusterId]; |
66 |
|
|
} |
67 |
|
✗ |
std::map<int, std::list<int>> cluster2ParticlesRestricted; |
68 |
|
✗ |
for (int i = 0; i < n; ++i) { |
69 |
|
✗ |
int newId = clusterIdMap[i]; |
70 |
|
✗ |
cluster2ParticlesRestricted[newId].splice(cluster2ParticlesRestricted[newId].end(), |
71 |
|
✗ |
cluster2Particles[i]); |
72 |
|
|
} |
73 |
|
|
// At this point, cluster ids are in [0, numClusters), but they are not yet sorted by the cluster |
74 |
|
|
// size |
75 |
|
✗ |
std::vector<int> clusterSizes(numClusters, 0); |
76 |
|
✗ |
for (int clusterId : particle2Cluster) { |
77 |
|
✗ |
clusterSizes[clusterId]++; |
78 |
|
|
} |
79 |
|
✗ |
std::vector<int> permutation(numClusters); |
80 |
|
✗ |
std::iota(permutation.begin(), permutation.end(), 0); |
81 |
|
✗ |
std::sort(permutation.begin(), permutation.end(), |
82 |
|
✗ |
[&clusterSizes](int a, int b) { return clusterSizes[a] > clusterSizes[b]; }); |
83 |
|
✗ |
std::vector<int> inversePermutation(numClusters); |
84 |
|
✗ |
for (int i = 0; i < numClusters; ++i) { |
85 |
|
✗ |
inversePermutation[permutation[i]] = i; |
86 |
|
|
} |
87 |
|
✗ |
for (int &clusterId : particle2Cluster) { |
88 |
|
✗ |
clusterId = inversePermutation[clusterId]; |
89 |
|
|
} |
90 |
|
✗ |
cluster2Particles.clear(); |
91 |
|
✗ |
for (int i = 0; i < numClusters; ++i) { |
92 |
|
✗ |
cluster2Particles[i].splice(cluster2Particles[i].end(), cluster2ParticlesRestricted[permutation[i]]); |
93 |
|
|
} |
94 |
|
|
} |
95 |
|
|
|
96 |
|
✗ |
PercolationAnalyzer::PercolationAnalyzer(System *system) |
97 |
|
✗ |
: system(system), tLast(std::numeric_limits<double>::quiet_NaN()), |
98 |
|
✗ |
clusterer(system->particles, system->L, system->interaction->internal->range) {} |
99 |
|
|
|
100 |
|
✗ |
void PercolationAnalyzer::update() { |
101 |
|
✗ |
if (system->t == tLast) { |
102 |
|
|
return; |
103 |
|
|
} |
104 |
|
✗ |
this->tLast = system->t; |
105 |
|
✗ |
int n = system->particles.size(); |
106 |
|
|
// Do the clustering algorithm |
107 |
|
✗ |
clusterer.reset(); |
108 |
|
✗ |
system->interaction->iteratePairs(clusterer, true, false); |
109 |
|
✗ |
clusterer.sortClusters(); |
110 |
|
|
// Do the coordination statistics |
111 |
|
|
// Just go through all particles and make a histogram of the coordinations |
112 |
|
✗ |
coordinationCount.clear(); |
113 |
|
✗ |
coordinationCount.resize(n + 1); |
114 |
|
✗ |
clusterSizeCount.clear(); |
115 |
|
✗ |
clusterSizeCount.resize(n + 1); |
116 |
|
✗ |
for (int i = 0; i < n; ++i) { |
117 |
|
✗ |
coordinationCount[i] = 0; |
118 |
|
✗ |
clusterSizeCount[i] = 0; |
119 |
|
|
} |
120 |
|
✗ |
for (int coord : clusterer.particle2Coordination) { |
121 |
|
✗ |
coordinationCount[coord]++; |
122 |
|
|
} |
123 |
|
✗ |
numClusters = clusterer.cluster2Particles.size(); |
124 |
|
✗ |
clusterSizeMax = numClusters == 0 ? 0 : clusterer.cluster2Particles[0].size(); |
125 |
|
✗ |
for (int i = 0; i < numClusters; ++i) { |
126 |
|
✗ |
clusterSizeCount[i] = clusterer.cluster2Particles[i].size(); |
127 |
|
|
} |
128 |
|
|
} |
129 |
|
|
|
130 |
|
✗ |
Analyzer *AnalyzerFactory::create(System *system, const std::string &key) { |
131 |
|
✗ |
if (key == "percolation") { |
132 |
|
✗ |
return new PercolationAnalyzer(system); |
133 |
|
|
} else { |
134 |
|
✗ |
throw std::invalid_argument("No " + key + " analyzer found"); |
135 |
|
|
} |
136 |
|
|
} |
137 |
|
|
|