From 2f3979fbca10725eb120a21a51afab6ffda2b0d4 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Tue, 10 Dec 2024 15:11:17 +0100 Subject: Removed percolation-specific data from graph class. --- percolation.cpp | 118 ++++++++++++++++++++++++++------------------------------ 1 file changed, 55 insertions(+), 63 deletions(-) diff --git a/percolation.cpp b/percolation.cpp index 08c6970..a5c2a36 100644 --- a/percolation.cpp +++ b/percolation.cpp @@ -14,128 +14,128 @@ using Vector = Eigen::Matrix; using Matrix = Eigen::Matrix; class Vertex; -class Edge; class HalfEdge { public: - Edge& parent; + unsigned index; Vertex& neighbor; - HalfEdge(Edge& e, Vertex& v) : parent(e), neighbor(v) {}; + HalfEdge(unsigned index, Vertex& v) : index(index), neighbor(v) {}; }; class Vertex { public: std::vector bonds; unsigned index; - unsigned cluster; }; -class Edge { -public: - bool active; -}; +class Cluster : public std::list> {}; -class Cluster : public std::list> {}; -class Clusters { -public: - std::vector lists; - std::vector indices; -}; +unsigned uPow(unsigned x, unsigned a) { + unsigned result = 1; + while (a) { + if (a & 1) { + result *= x; + } + a >>= 1; + x *= x; + } + return result; +} class Graph { public: unsigned D; unsigned L; std::vector vertices; - std::vector edges; - std::vector clusters; + std::vector multiplicities; unsigned Nv() const { return vertices.size(); } unsigned Ne() const { - return edges.size(); + return D * Nv(); } unsigned squaredDistance(unsigned vi, unsigned vj) const { unsigned sd = 0; for (unsigned i = 0; i < D; i++) { - int x1 = (unsigned)fmod((unsigned)(vi / (unsigned)pow(L, i)), pow(L, i + 1)); - int x2 = (unsigned)fmod((unsigned)(vj / (unsigned)pow(L, i)), pow(L, i + 1)); + int x1 = (vi / uPow(L, i)) % uPow(L, i + 1); + int x2 = (vj / uPow(L, i)) % uPow(L, i + 1); unsigned Δx = abs(x1 - x2); if (Δx > L / 2) { Δx = L - Δx; } - sd += pow(Δx, 2); + sd += Δx * Δx; } return sd; } /* Initialize a square lattice */ - Graph(unsigned D, unsigned L) : D(D), L(L), vertices(pow(L, D)), edges(D * pow(L, D)) { + Graph(unsigned D, unsigned L) : D(D), L(L), vertices(uPow(L, D)), multiplicities(D * L * L / 4) { for (unsigned i = 0; i < Nv(); i++) { vertices[i].index = i; vertices[i].bonds.reserve(2 * D); for (unsigned j = 0; j < D; j++) { - unsigned n1 = (unsigned)pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) + - (unsigned)fmod(i + pow(L, j), pow(L, j + 1)); - unsigned n2 = (unsigned)pow(L, j + 1) * (i / ((unsigned)pow(L, j + 1))) + - (unsigned)fmod(pow(L, j + 1) + i - pow(L, j), pow(L, j + 1)); + unsigned n1 = uPow(L, j + 1) * (i / uPow(L, j + 1)) + (i + uPow(L, j)) % uPow(L, j + 1); + unsigned n2 = uPow(L, j + 1) * (i / uPow(L, j + 1)) + (uPow(L, j + 1) + i - uPow(L, j)) % uPow(L, j + 1); unsigned e1 = j * Nv() + i; unsigned e2 = j * Nv() + n2; - HalfEdge he1(edges[e1], vertices[n1]); - HalfEdge he2(edges[e2], vertices[n2]); + HalfEdge he1(e1, vertices[n1]); + HalfEdge he2(e2, vertices[n2]); vertices[i].bonds.push_back(he1); vertices[i].bonds.push_back(he2); } } - } - void dilute(Real p, Rng& r) { - for (Edge& e : edges) { - e.active = p > r.variate(0.0, 1.0); + for (const Vertex& v1 : vertices) { + for (const Vertex& v2 : vertices) { + unsigned dist = squaredDistance(v1.index, v2.index); + if (dist > 0) { + multiplicities[dist - 1]++; + } + } } } - void markClusters() { - clusters = {}; - for (Vertex& v : vertices) { - v.cluster = 0; - }; + std::vector markClusters(const std::vector& activeEdges) const { + std::vector clusters; + std::vector indicies(Nv()); unsigned nClusters = 0; - for (Vertex& v : vertices) { - if (v.cluster == 0) { + for (const Vertex& v : vertices) { + if (indicies[v.index] == 0) { nClusters++; clusters.push_back({}); - std::stack> inCluster; + std::stack> inCluster; inCluster.push(v); while (!inCluster.empty()) { - Vertex& vn = inCluster.top(); + const Vertex& vn = inCluster.top(); inCluster.pop(); - if (vn.cluster == 0) { - vn.cluster = nClusters; + if (indicies[vn.index] == 0) { + indicies[vn.index] = nClusters; clusters[nClusters - 1].push_back(vn); } - for (HalfEdge& e : vn.bonds) { - if (e.parent.active && e.neighbor.cluster == 0) { + for (const HalfEdge& e : vn.bonds) { + if (activeEdges[e.index] && indicies[e.neighbor.index] == 0) { inCluster.push(e.neighbor); } } } } } + + return clusters; } - Matrix laplacian() const { + Matrix laplacian(const std::vector& activeEdges) const { Matrix L = Matrix::Zero(Nv(), Nv()); for (const Vertex& v : vertices) { for (const HalfEdge& e : v.bonds) { - if (e.parent.active) { + if (activeEdges[e.index]) { L(v.index, v.index) += 1; L(v.index, e.neighbor.index) -= 1; } @@ -175,27 +175,19 @@ int main(int argc, char* argv[]) { } } - std::vector multiplicities(D * L * L / 4); - Graph G(D, L); - for (const Vertex& v1 : G.vertices) { - for (const Vertex& v2 : G.vertices) { - unsigned dist = G.squaredDistance(v1.index, v2.index); - if (dist > 0) { - multiplicities[dist - 1]++; - } - } - } - for (unsigned trials = 0; trials < n; trials++) { - G.dilute(p, r); - G.markClusters(); - Matrix laplacian = G.laplacian(); + std::vector activeEdges(G.Ne()); + for (short unsigned& edge : activeEdges) { + edge = r.uniform(0.0, 1.0) < p; + } + std::vector clusters = G.markClusters(activeEdges); + Matrix laplacian = G.laplacian(activeEdges); std::vector conductivities(L * L * D / 4); - for (const Cluster& c : G.clusters) { + for (const Cluster& c : clusters) { std::vector inds; inds.reserve(c.size()); for (const Vertex& v : c){ @@ -221,8 +213,8 @@ int main(int argc, char* argv[]) { } for (unsigned i = 0; i < conductivities.size(); i++) { - if (multiplicities[i] != 0) { - std::cout << conductivities[i] / multiplicities[i] << " "; + if (G.multiplicities[i] != 0) { + std::cout << conductivities[i] / G.multiplicities[i] << " "; } else { std::cout << 0 << " "; } -- cgit v1.2.3-70-g09d2