| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #include "percolation.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 | |||
| 10 | ✗ | Clusterer::Clusterer(std::vector<Particle> &particles, const Vec &L, double cutoff, | |
| 11 | std::vector<int> &particle2Coordination, std::vector<int> &particle2Cluster, | ||
| 12 | ✗ | std::map<int, std::list<int>> &cluster2Particles) | |
| 13 | ✗ | : particles(particles), L(L), cutoff(cutoff), cutoffsq(cutoff * cutoff), | |
| 14 | ✗ | particle2Coordination(particle2Coordination), particle2Cluster(particle2Cluster), | |
| 15 | ✗ | cluster2Particles(cluster2Particles) {} | |
| 16 | |||
| 17 | ✗ | double Clusterer::operator()(Particle &p1, Particle &p2) { | |
| 18 | ✗ | int i = std::distance(&particles.front(), &p1); | |
| 19 | ✗ | int j = std::distance(&particles.front(), &p2); | |
| 20 | ✗ | double distsq = r_ij(L, particles[i].r, particles[j].r).squaredNorm(); | |
| 21 | ✗ | if (distsq < cutoffsq && particle2Cluster[i] != particle2Cluster[j]) { | |
| 22 | ✗ | particle2Coordination[i]++; | |
| 23 | ✗ | particle2Coordination[j]++; | |
| 24 | ✗ | int oldCluster = std::max(particle2Cluster[i], particle2Cluster[j]); | |
| 25 | ✗ | int newCluster = std::min(particle2Cluster[i], particle2Cluster[j]); | |
| 26 | ✗ | for (int i : cluster2Particles[oldCluster]) { | |
| 27 | ✗ | particle2Cluster[i] = newCluster; | |
| 28 | } | ||
| 29 | ✗ | cluster2Particles[newCluster].splice(cluster2Particles[newCluster].end(), | |
| 30 | ✗ | cluster2Particles[oldCluster]); | |
| 31 | } | ||
| 32 | ✗ | return 0.; | |
| 33 | } | ||
| 34 | |||
| 35 | 75 | Percolation::Percolation(System *system) | |
| 36 | 75 | : system(system), tLast(std::numeric_limits<double>::quiet_NaN()) {} | |
| 37 | |||
| 38 | ✗ | void Percolation::update() { | |
| 39 | ✗ | if (system->t == tLast) { | |
| 40 | ✗ | return; | |
| 41 | } | ||
| 42 | // Initialize data structures | ||
| 43 | ✗ | this->tLast = system->t; | |
| 44 | ✗ | int n = system->particles.size(); | |
| 45 | ✗ | numClusters = 0; | |
| 46 | ✗ | coordinationCount.clear(); | |
| 47 | ✗ | coordinationCount.resize(n + 1); | |
| 48 | ✗ | clusterSizeCount.clear(); | |
| 49 | ✗ | clusterSizeCount.resize(n + 1); | |
| 50 | ✗ | particle2Coordination.resize(n); | |
| 51 | ✗ | std::fill(particle2Coordination.begin(), particle2Coordination.end(), 0); | |
| 52 | ✗ | particle2Cluster.resize(n); | |
| 53 | ✗ | for (int i = 0; i < n; ++i) { | |
| 54 | ✗ | particle2Cluster[i] = i; | |
| 55 | ✗ | cluster2Particles[i] = {i}; | |
| 56 | } | ||
| 57 | // Do the clustering algorithm | ||
| 58 | ✗ | Clusterer clusterer(system->particles, system->L, system->interaction->internal->range, | |
| 59 | ✗ | particle2Coordination, particle2Cluster, cluster2Particles); | |
| 60 | ✗ | system->interaction->iteratePairs(clusterer, true, false); | |
| 61 | // Do the coordination statistics | ||
| 62 | // Just go through all particles and make a histogram of the coordinations | ||
| 63 | ✗ | for (int coord : particle2Coordination) { | |
| 64 | ✗ | coordinationCount[coord]++; | |
| 65 | } | ||
| 66 | // Do the cluster sizes | ||
| 67 | // Go through all cluster lists and get their sizes | ||
| 68 | ✗ | for (auto &[_, particlesInCluster] : cluster2Particles) { | |
| 69 | ✗ | int clusterSize = particlesInCluster.size(); | |
| 70 | ✗ | if (clusterSize > 0) { | |
| 71 | ✗ | numClusters++; | |
| 72 | ✗ | clusterSizeCount[clusterSize] += clusterSize; | |
| 73 | } | ||
| 74 | } | ||
| 75 | // coordinationCount and clusterSizeCount are normed to the number of particles | ||
| 76 | } | ||
| 77 |