diff --git a/postprocessing/Makefile b/postprocessing/Makefile
index 69cb856..3498a34 100644
--- a/postprocessing/Makefile
+++ b/postprocessing/Makefile
@@ -1,97 +1,97 @@
IDIR =./include
CXX=g++
#CPPFLAGS=-g -I$(IDIR) -I/home/samarth/codes/voronoi/include -std=c++11 -O3 -w -Wall
-#CPPFLAGS=-g -I$(IDIR) -O3
-CPPFLAGS=-I$(IDIR) -O3
-LDFLAGS=-O3
+CPPFLAGS=-g -I$(IDIR) -O3
+#CPPFLAGS=-I$(IDIR) -O3
+#LDFLAGS=-O3
POSTFIX=-L/home/samarth/codes/voronoi/lib/ -lvoro++
SRC_DIR = ./src
OBJ_DIR = ./obj
SRC_FILES = $(wildcard $(SRC_DIR)/*.cpp)
OBJ_FILES = $(patsubst $(SRC_DIR)/%.cpp,$(OBJ_DIR)/%.o,$(SRC_FILES))
BIN_SRC_DIR = ./bin_src
BIN_DIR = ./bin
all: unfold radius_of_gyration sq_vs_q pair_correlation percolation load_bearing_paths lr_iq lbp_brute_force percolation_invA chemical_distance lb_bonds_clusterwise_invA lbp_bfs
#all: psd_test
density_correlation: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o density_correlation.o target/density_correlation.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ density_correlation.o
unfold: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o unfold.o target/unfold.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ unfold.o
radius_of_gyration: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o radius_of_gyration.o target/radius_of_gyration.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ radius_of_gyration.o
sq_vs_q: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o sq_vs_q.o target/sq_vs_q.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ sq_vs_q.o
pair_correlation: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o pair_correlation.o target/pair_correlation.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ pair_correlation.o
psd_test: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o psd_test.o target/psd_test.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ psd_test.o
percolation: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o percolation.o target/percolation.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ percolation.o
load_bearing_paths: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o load_bearing_paths.o target/load_bearing_paths.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ load_bearing_paths.o
lbp_brute_force: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o lbp_brute_force.o target/lbp_brute_force.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ lbp_brute_force.o
percolation_invA: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o percolation_invA.o target/percolation_invA.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ percolation_invA.o
lb_bonds_invA: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o lb_bonds_invA.o target/lb_bonds_invA.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ lb_bonds_invA.o
lb_bonds_clusterwise_invA: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o lb_bonds_clusterwise_invA.o target/lb_bonds_clusterwise_invA.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ lb_bonds_clusterwise_invA.o
lb_bonds_via_cg: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o lb_bonds_via_cg.o target/lb_bonds_via_cg.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ lb_bonds_via_cg.o
chemical_distance: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o chemical_distance.o target/chemical_distance.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ chemical_distance.o
lr_iq: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o lr_iq.o target/lr_iq.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ lr_iq.o
lbp_bfs: $(OBJ_FILES) | $(BIN_DIR)
$(CXX) $(CPPFLAGS) -c -o lbp_bfs.o target/lbp_bfs.cpp
$(CXX) $(LDFLAGS) -o $(BIN_DIR)/$@ $^ lbp_bfs.o
$(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp | $(OBJ_DIR)
$(CXX) $(CPPFLAGS) -c -o $@ $<
$(BIN_DIR):
mkdir $(BIN_DIR)
$(OBJ_DIR):
mkdir $(OBJ_DIR)
clean:
rm *.o
rm -rf $(BIN_DIR)
rm -rf $(OBJ_DIR)
diff --git a/postprocessing/include/post.h b/postprocessing/include/post.h
index 5b8ce05..0ef9c35 100644
--- a/postprocessing/include/post.h
+++ b/postprocessing/include/post.h
@@ -1,407 +1,423 @@
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <sys/types.h>
#include <sys/stat.h>
#include <vector>
#include <sstream>
#include <fstream>
#include <random>
//#include <voro++/voro++.hh>
#include <map>
#include <algorithm>
#include <bits/stdc++.h>
#include <Eigen/Sparse>
#include <chrono>
#include <sys/resource.h>
#ifndef POSTPROCESSING_H
#define POSTPROCESSING_H
namespace post_p{
class postprocessing
{
public:
postprocessing(char *config_filename);
postprocessing(char *config_filename, bool build);
~postprocessing();
void read_params();
void read_config();
void memory_allocation_function();
std::vector<std::string> split_string_by_delimiter(const std::string& s, char delimiter);
void read_params_parser(char *config_filename);
void read_config_parser(char *config_filename);
void read_params_for_percolation_build(char *config_filename);
int dim() const;
int numParticles() const;
int maxAttachments() const;
int totalClusters() const;
double lo_hi(const int axis, const int limit) const;
double& lo_hi(const int axis, const int limit);
double L(const int axis) const;
double& L(const int axis);
double halfL(const int axis) const;
double& halfL(const int axis);
int periodic(const int axis) const;
int& periodic(const int axis);
double posDiff(const int axis) const;
double& posDiff(const int axis);
double pos(const int i, const int axis) const;
double& pos(const int i, const int axis);
int numAttachments(const int i) const;
int& numAttachments(const int i);
int original_seed(const int i) const;
int& original_seed(const int i);
int current_seed(const int i) const;
int& current_seed(const int i);
double diameter(const int i) const;
double& diameter(const int i);
double radius(const int i) const;
double& radius(const int i);
int unfoldedAttachments(const int i) const;
int& unfoldedAttachments(const int i);
int attachment(const int i, const int j) const;
int& attachment(const int i, const int j);
bool is_placed(const int i) const;
bool& is_placed(const int i);
bool attachments_placed(const int i) const;
bool& attachments_placed(const int i);
double unfolded_coords(const int i, const int axis) const;
double& unfolded_coords(const int i, const int axis);
double delta_coords(const int i, const int axis) const;
double& delta_coords(const int i, const int axis);
int delta_hist(const int i) const;
int& delta_hist(const int i);
int cluster_percolation(const int i, const int axis) const;
int& cluster_percolation(const int i, const int axis);
int load_bearing_paths(const int i, const int axis) const;
int& load_bearing_paths(const int i, const int axis);
double get_periodic_image(double x, const int axis);
void print_positions();
void calc_rij();
void dump_rij_file();
void dump_rij_file(char *filename);
void dump_rij_hist_file(double bin_size);
void dump_rij_hist_file(double bin_size, char *filename);
void calc_rij_hist(double bin_size);
void dump_unfolded_file();
void dump_unfolded_file(char *filename);
void dump_density_correlation(double bin_size);
void dump_density_correlation(double bin_size, char *filename);
void dump_percolation_file();
void dump_percolation_file(char *filename);
void dump_percolation_via_invA_file();
void dump_percolation_via_invA_file(char *filename);
void dump_load_bearing_paths_file(char *filename);
//void dump_load_bearing_paths_file(int origin, char *filename);
void dump_load_bearing_paths_file();
void dump_rog();
void dump_rog(char *filename);
void init_unfolding();
void init_unfolding_for_lbp();
void unfold(const int prev, const int next);
void unfold_for_lbp(const int prev, const int next, bool build_pb);
void calc_unfolded_rij();
void dump_unfolded_hist_file(double bin_size);
void calc_unfolded_hist(double bin_size);
void dump_scattering_function(double q_min, double q_max, int num_q);
void dump_scattering_function(double q_min, double q_max, int num_q, char *filename);
//void dump_lr_scattering_function(double q_min, double q_max, int num_q, double bin_size);
void dump_lr_scattering_function(double q_min, double q_max, int num_q, double bin_size, char *filename);
void calc_scattering_function(double q_min, double q_max, int num_q);
void calc_density_correlation(double bin_size);
void save_config();
void save_unfolded_config();
void save_unfolded_config(char *filename);
void get_headers();
bool check_percolation();
double calc_rog();
void psd(int total_iters);
void save_radius_distribution();
void save_radius_distribution(char *filename);
void init_unfolding_without_recursion();
void place_attachments(int i);
void build_bond_map();
void percolation_detection();
void reset_unfolding_params();
void switch_off_bonds(std::pair<int,int> bond);
void switch_on_bonds(std::pair<int,int> bond);
void create_subsystem();
void erase_subsystem();
void print_lbp();
void postprocess_lbp();
void switch_off_alt_lbp(int i);
void reset_bond_map(bool status);
void activate_path(std::vector<std::pair<int,int>> bonds, bool status);
void shred_path(int i);
void isolation_routine();
void martin_test(char *lb_file);
void lbp_brute_force();
void makeCombiUtil(int n, int left, int k);
void makeCombi(int n, int k);
void dump_percolation_via_inverse();
void percolation_via_inverse();
void modify_coords();
void modify_coords_after_minimization(int axis);
void build_A();
void build_A(int axis);
void build_b(int axis);
void calculate_bond_lengths_direction_wise(int axis);
void calculate_bond_lengths();
void determine_LB_bonds();
void copy_positions();
void copy_positions_for_cluster();
void unfold_for_clusterwise(int prev, int next);
void determine_LB_bonds_clusterwise(char *filename);
+ void determine_LB_bonds_network(char *filename);
void determine_LB_bonds_clusterwise_directionwise(char *filename);
void determine_LB_bonds_clusterwise_via_cg(char *filename);
bool check_if_particles_placed();
void modify_coords_for_cluster();
void modify_coords_after_minimization_for_cluster(int axis);
void build_A_for_cluster();
void build_A_for_cluster(int axis);
void build_b_for_cluster(int axis);
void build_A_for_cluster_for_cg();
void build_b_for_cluster_for_cg(int axis);
void calculate_bond_lengths_direction_wise_for_cluster(int axis);
void calculate_bond_lengths_for_cluster();
void calculate_bond_lengths_for_cluster(int axis);
void dump_lb_bonds_for_cluster_via_invA(char *filename);
void dump_lb_bonds_for_cluster_via_cg(char *filename);
void dump_chemical_distance(char *lb_filename, char *cd_filename);
void reset_chemical_distance_index();
void find_chemical_distance(int i, int j);
void find_shortest_path(int i, int j);
void find_paths_for_all_bonds(char *lb_file, char *cd_filename);
void find_all_paths(int source, int destination);
void find_all_paths_recursion(int source, int destination);
void distort_positions();
void increase_stiffness();
void copy_b();
void calculate_b_diff();
void put_particles_back_in_box(int axis);
void dump_xyz_stepwise(char *filename, int file_num);
void set_lbb_params(char *filename);
void init_lbb_unfolding();
void init_lbb_cluster_matrices();
void print_lbb_info();
void modify_A_matrix();
void switch_off_transverse_forces(int axis);
void print_minimized_coords(char *filename, int filenum);
void init_lbb_unfolding_without_recursion();
void unfold_for_clusterwise_without_recursion(int i);
void calc_total_bonds();
void print_mem_usage();
void print_cluster_vals(char *filename);
+ void unfold_for_lbp_without_recursion(int i);
+ void comb(int N, int K);
+
+ void set_lbb_params_for_network(char *filename);
+ void print_lbb_info_for_network(char *filename);
+ void dump_LB_bonds_for_network(char *filename);
+ void reset_unfolding_for_network();
+ void add_ghost_particle_to_network_cluster();
+ void fix_A_for_network_cluster();
+ void fix_b_for_network_cluster();
+ void build_A_for_network_cluster();
+ void build_b_for_network_cluster();
+ void init_lbb_cluster_matrices_for_network_cluster();
+ void calculate_bond_lengths_for_network_cluster();
+ bool check_sink_and_source();
private:
std::mt19937 generator;
std::uniform_real_distribution<double> dis;
int headers_;
int N_;
int D_;
int max_attachments_;
int folded_;
int lattice_;
int iters_;
int columns_;
double seed_mass_;
double alpha_;
double phi_;
double *lo_hi_;
double *L_;
double *halfL_;
double *diameter_;
double *radius_;
int *periodic_;
double *posDiff_;
int *posDiff_int_;
double *pos_;
int *num_attachments_;
int *attachment_;
int *original_seed_;
int *current_seed_;
double *r_ij_;
double *r_ij_hist_;
int *rho_hist_;
bool *is_placed_;
bool *attachments_placed_;
double *unfolded_coords_;
int temp_next;
double *delta_coords_;
double q_min;
double q_max;
double dq;
int num_q;
double *sq_;
double *lr_sq_;
int *unfolded_num_attachments_;
int *delta_hist_;
int max_delta_;
int *cluster_percolation_;
int *load_bearing_paths_;
int *temp_lbp_;
int totalClusters_;
int N_pairs_;
int r_ij_hist_bins_;
bool lattice_flag_ = false;
bool N_flag_ = false;
bool D_flag_ = false;
bool max_attachments_flag_ = false;
bool L_flag_ = false;
bool periodic_flag_ = false;
bool folded_flag_ = false;
bool iters_flag_ = false;
bool seed_mass_flag_ = false;
bool alpha_flag_ = false;
bool phi_flag_ = false;
bool columns_flag_ = false;
double *radius_dis;
double *centres;
int *type_dis;
int max_psd_iters;
int total_lbp;
std::map<std::pair<int,int>, int> bond_map_status;
std::vector<std::pair<int,int>> percolating_bonds;
std::pair<int,int> temp_pair;
int *divisors;
int *bond_bit_repr;
std::vector<int> path_components;
std::vector<std::pair<int,int>> current_path;
bool path_percolation;
std::vector<std::vector<std::pair<int,int>>> final_percolating_bonds;
std::vector<std::vector<std::pair<int,int>>> weak_links;
std::vector<bool> pb_status;
std::vector<std::pair<int,int>> unique_bonds;
std::vector<int> tmp;
std::vector<int> ans;
Eigen::SparseMatrix<double> A;
Eigen::VectorXd x;
Eigen::VectorXd b;
Eigen::VectorXd old_b;
Eigen::VectorXd b_diff;
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> solver;
Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower|Eigen::Upper> cg;
Eigen::MatrixXd modified_folded_x;
Eigen::VectorXd ref_pos;
Eigen::VectorXd bond_lengths;
Eigen::MatrixXd bond_lengths_direction_wise;
Eigen::VectorXd check_for_solver_mk;
int num_bonds;
int num_bonds_for_cluster;
int num_particles_for_cluster;
std::vector<std::pair<int,int>> unique_bonds_for_cluster;
std::map<int, int> particles_to_index;
std::map<int, int> index_to_particles;
int index_i;
int index_j;
int *chemical_distance_index;
std::vector<int> ref_nodes;
std::vector<int> neighbour_nodes;
std::vector<std::string> full_lb_list;
std::vector<std::string> lb_bonds_str;
bool *in_path;
std::vector<int> source_to_destination_path;
int current_path_index;
FILE *f_bfs;
double stiffness_coef;
double inv_stiffness_coef;
int *over_the_boundary;
double max_length;
int max_row;
int max_col;
double lbp_tolerance;
int nnz;
int a_i;
int a_j;
FILE *f_lbp;
/*std::chrono::steady_clock::time_point cp_1;
std::chrono::steady_clock::time_point cp_2;
std::chrono::steady_clock::time_point cp_3;
std::chrono::steady_clock::time_point cp_4;
std::chrono::steady_clock::time_point cp_5;*/
std::chrono::steady_clock::time_point cp_A1;
std::chrono::steady_clock::time_point cp_A2;
std::chrono::steady_clock::time_point cp_A3;
std::chrono::steady_clock::time_point cp_A4;
double time_A1 = 0.;
double time_A2 = 0.;
//Eigen::Triplet<double> T;
std::vector<Eigen::Triplet<double>> tripleList;
int n_lbb = 0;
int particle_counter;
std::vector<int> to_build_list;
int i_map;
int j_map;
int total_num_bonds;
std::chrono::steady_clock::time_point unf_cp_1;
std::chrono::steady_clock::time_point unf_cp_2;
std::chrono::steady_clock::time_point unf_cp_3;
std::chrono::steady_clock::time_point unf_cp_4;
std::chrono::steady_clock::time_point unf_cp_5;
std::chrono::steady_clock::time_point unf_cp_6;
std::chrono::steady_clock::time_point unf_cp_7;
double unf_time_1 = 0.;
double unf_time_2 = 0.;
double unf_time_3 = 0.;
double unf_time_4 = 0.;
double unf_time_5 = 0.;
double unf_time_6 = 0.;
int unf_checkpoint = 0;
struct rusage myusage;
long baseline;
};
}
#endif
diff --git a/postprocessing/src/lb_bonds_clusterwise.cpp b/postprocessing/src/lb_bonds_clusterwise.cpp
index ad0fff2..b9cfdf8 100644
--- a/postprocessing/src/lb_bonds_clusterwise.cpp
+++ b/postprocessing/src/lb_bonds_clusterwise.cpp
@@ -1,893 +1,862 @@
#include <post.h>
namespace post_p
{
void postprocessing::print_mem_usage()
{
getrusage(RUSAGE_SELF, &myusage);
std::cout<<" = "<<(1. * myusage.ru_maxrss)/1e3<<std::endl;
}
void postprocessing::dump_lb_bonds_for_cluster_via_invA(char *filename)
{
//std::cout<<"1 new"<<std::endl;
//getrusage(RUSAGE_SELF, &myusage);
//baseline = myusage.ru_maxrss;
//std::cout<<"on entering dump function = "<<(1. * baseline)/1e3<<std::endl;;
determine_LB_bonds_clusterwise(filename);
//std::cout<<"2"<<std::endl;
}
void postprocessing::calc_total_bonds()
{
total_num_bonds = 0;
for (int i = 0; i < numParticles(); i++)
total_num_bonds += numAttachments(i);
total_num_bonds /= 2;
}
bool postprocessing::check_if_particles_placed()
{
bool all_particles_placed = true;
for (int i = unf_checkpoint; i < numParticles(); i++) {
//all_particles_placed = all_particles_placed && is_placed(i);
if (is_placed(i) == false)
return false;
}
- //std::cout<<"bool values = "<<all_particles_placed<<std::endl;
-
return true;
}
void postprocessing::distort_positions()
{
for (int i = 0; i < numParticles(); i++){
for (int axis = 0; axis < dim(); axis++){
//std::cout<<"old position = "<<pos(i,axis)<<"\t";
pos(i,axis) += ((0.1 * dis(generator))-0.05);
//std::cout<<"new position = "<<pos(i,axis)<<"\t";
}
//std::cout<<"\n";
}
//std::cout<<"distorted"<<std::endl;
}
void postprocessing::set_lbb_params(char *filename)
{
reset_unfolding_params();
build_bond_map();
reset_bond_map(true);
total_lbp = 0;
ref_pos.resize(dim());
totalClusters_ = 0;
f_lbp = fopen(filename, "w");
fprintf(f_lbp,"clusterNumber,bond,");
for (int axis = 0; axis < dim(); axis++)
fprintf(f_lbp,"bond_component_x%d,",axis);
fprintf(f_lbp,"bond_length\n");
- /*if (!A.isCompressed()){
- A.makeCompressed();
- }*/
-
lbp_tolerance = 1e-6;
-
}
void postprocessing::init_lbb_unfolding()
{
num_bonds_for_cluster = 0;
num_particles_for_cluster = 0;
unique_bonds.clear();
index_to_particles.clear();
particles_to_index.clear();
for (int i = 0; i < numParticles(); i++){
if (is_placed(i) == false){
particles_to_index[i] = 0;
index_to_particles[0] = i;
num_particles_for_cluster += 1;
unfold_for_clusterwise(i,i);
break;
}
}
//std::cout<<"num_particles_for_cluster = "<<num_particles_for_cluster<<"\t num_bonds = "<<num_bonds_for_cluster<<std::endl;
}
void postprocessing::init_lbb_unfolding_without_recursion()
{
- //unf_cp_1 = std::chrono::steady_clock::now();
-
num_bonds_for_cluster = 0;
num_particles_for_cluster = 0;
unique_bonds.clear();
index_to_particles.clear();
particles_to_index.clear();
to_build_list.clear();
- //unf_cp_2 = std::chrono::steady_clock::now();
-
- //unf_time_1 += std::chrono::duration_cast<std::chrono::nanoseconds>(unf_cp_2 - unf_cp_1).count();
-
- //unf_cp_3 = std::chrono::steady_clock::now();
-
for (int i = unf_checkpoint; i < numParticles(); i++){
if (is_placed(i) == false){
particles_to_index[i] = num_particles_for_cluster;
index_to_particles[num_particles_for_cluster] = i;
to_build_list[num_particles_for_cluster] = i;
num_particles_for_cluster++;
unfold_for_clusterwise_without_recursion(i);
is_placed(i) = true;
for (int temp_index = 1; temp_index < num_particles_for_cluster; temp_index++){
unfold_for_clusterwise_without_recursion(to_build_list[temp_index]);
}
unf_checkpoint = i;
break;
}
}
-
- //unf_cp_4 = std::chrono::steady_clock::now();
-
- //unf_time_2 += std::chrono::duration_cast<std::chrono::nanoseconds>(unf_cp_4 - unf_cp_3).count();
-
- //std::cout<<"num_particles_for_cluster = "<<num_particles_for_cluster<<"\t num_bonds = "<<num_bonds_for_cluster<<std::endl;
-
}
void postprocessing::unfold_for_clusterwise_without_recursion(int i)
{
for (int att = 0; att < numAttachments(i); att++)
{
-
index_j = attachment(i, att);
- if (attachments_placed(index_j) == false)
+ if ((attachments_placed(index_j) == false) && (bond_map_status[{i,index_j}] == 1))
{
if (is_placed(index_j) == false){
particles_to_index[index_j] = num_particles_for_cluster;
index_to_particles[num_particles_for_cluster] = index_j;
is_placed(index_j) = true;
- //to_build_list.push_back(index_j);
to_build_list[num_particles_for_cluster] = index_j;
num_particles_for_cluster++;
-
- //std::cout<<"index_i = "<<i<<"\t index_j = "<<index_j<<std::endl;
}
i_map = particles_to_index[i];
j_map = particles_to_index[index_j];
- //unique_bonds.push_back({i_map,j_map});
unique_bonds[num_bonds_for_cluster] = {i_map,j_map};
num_bonds_for_cluster += 1;
- //std::cout<<"i = "<<i<<"\t j = "<<index_j<<std::endl;
-
}
}
attachments_placed(i) = true;
}
void postprocessing::unfold_for_clusterwise(const int prev, const int next)
{
for (int axis = 0; axis < dim(); axis++)
posDiff(axis) = get_periodic_image(pos(next,axis) - pos(prev,axis), axis);
for (int axis = 0; axis < dim(); axis++){
if (prev == next)
unfolded_coords(next,axis) = pos(prev,axis) + posDiff(axis);
else
unfolded_coords(next,axis) = unfolded_coords(prev,axis) + posDiff(axis);
}
is_placed(next) = true;
for (int att = 0; att < numAttachments(next); att++)
{
temp_next = attachment(next, att);
if (is_placed(temp_next) == false){
particles_to_index[temp_next] = num_particles_for_cluster;
index_to_particles[num_particles_for_cluster] = temp_next;
num_particles_for_cluster++;
/*if (particles_to_index.count(next) != 1){
std::cout<<"unfold 1"<<std::endl;
}
if (particles_to_index.count(temp_next) != 1){
std::cout<<"unfold 2"<<std::endl;
}*/
index_i = particles_to_index[next];
index_j = particles_to_index[temp_next];
if (std::find(unique_bonds.begin(), unique_bonds.end(), std::make_pair(index_j,index_i)) == unique_bonds.end()){
unique_bonds.push_back({index_i,index_j});
num_bonds_for_cluster += 1;
}
unfoldedAttachments(next) = unfoldedAttachments(next) + 1;
unfoldedAttachments(temp_next) = unfoldedAttachments(temp_next) + 1;
unfold_for_clusterwise(next, temp_next);
}
else {
index_i = particles_to_index[next];
index_j = particles_to_index[temp_next];
if (std::find(unique_bonds.begin(), unique_bonds.end(), std::make_pair(index_j,index_i)) == unique_bonds.end()){
unique_bonds.push_back({index_i,index_j});
num_bonds_for_cluster += 1;
}
for (int axis = 0; axis < dim(); axis++) {
posDiff(axis) = unfolded_coords(temp_next, axis) - unfolded_coords(next, axis);
if ((posDiff(axis) > halfL(axis)) || (posDiff(axis) < -halfL(axis))){
cluster_percolation(totalClusters_,axis) = 1;
load_bearing_paths_[axis]++;
}
}
}
}
}
void postprocessing::init_lbb_cluster_matrices()
{
modified_folded_x.resize(num_particles_for_cluster,dim());
copy_positions_for_cluster();
nnz = 2 * num_bonds_for_cluster + num_particles_for_cluster;
A.reserve(nnz);
A.resize((num_particles_for_cluster-1), (num_particles_for_cluster-1));
b.resize((num_particles_for_cluster-1));
old_b.resize((num_particles_for_cluster-1));
x.resize((num_particles_for_cluster-1));
bond_lengths_direction_wise.resize(num_bonds_for_cluster, dim());
bond_lengths.resize(num_bonds_for_cluster);
A.setZero();
build_A_for_cluster();
//A.makeCompressed();
}
void postprocessing::print_lbb_info()
{
max_length = bond_lengths.maxCoeff(&max_row);
if (max_length > lbp_tolerance){
switch_off_bonds({index_to_particles[unique_bonds[max_row].first], index_to_particles[unique_bonds[max_row].second]});
modify_A_matrix();
fprintf(f_lbp, "%d,%d-%d,",totalClusters_,index_to_particles[unique_bonds[max_row].first],index_to_particles[unique_bonds[max_row].second]);
for (int print_axis = 0; print_axis < dim(); print_axis++)
fprintf(f_lbp, "%lf,", sqrt(bond_lengths_direction_wise(max_row,print_axis)));
fprintf(f_lbp, "%lf\n", max_length);
n_lbb++;
//std::cout<<"num_c = "<<num_particles_for_cluster<<"\t n_lbb = "<<n_lbb<<std::endl;
}
}
void postprocessing::print_minimized_coords(char *filename, int filenum)
{
std::string file_test(filename);
/*std::vector<std::string> results;
results = split_string_by_delimiter(file_test, '.');
std::string full_filename = results[0] + "_" + std::to_string(filenum) + "." + results[1];
for (int i = 0; i < results.size(); i++)
std::cout<<"results "<<i<<" = "<<results[i]<<std::endl;*/
file_test = file_test.substr(0, file_test.size()-4);
std::string full_filename = file_test + "_" + std::to_string(filenum) + ".csv";
std::ofstream f_coords;
f_coords.open(full_filename);
f_coords << "id,";
for (int axis = 0; axis < dim(); axis++){
if (axis == (dim() - 1)){
f_coords<<"x"<<axis<<"\n";
}
else{
f_coords<<"x"<<axis<<",";
}
}
int p_index;
for (int i = 0; i < num_particles_for_cluster; i++){
//p_index = particles_to_index[i];
p_index = index_to_particles[i];
//std::cout<<"i = "<<i<<"\t p_index = "<<std::endl;
f_coords<<p_index<<",";
for (int axis = 0; axis < dim(); axis++){
if (axis == (dim()-1)){
//f_coords<<modified_folded_x(p_index,axis)<<"\n";
f_coords<<modified_folded_x(i,axis)<<"\n";
}
else{
//f_coords<<modified_folded_x(p_index,axis)<<",";
f_coords<<modified_folded_x(i,axis)<<",";
}
}
}
f_coords.close();
}
void postprocessing::print_cluster_vals(char *filename)
{
std::string file_test(filename);
file_test = file_test.substr(0, file_test.size()-4);
std::string full_filename = file_test + "_clusters.csv";
std::ofstream f_coords;
f_coords.open(full_filename, std::ios_base::app);
for (int i = 0; i < num_particles_for_cluster; i++)
f_coords<<to_build_list[i]<<","<<totalClusters_<<"\n";
f_coords.close();
}
void postprocessing::determine_LB_bonds_clusterwise(char *filename)
{
std::chrono::steady_clock::time_point cp_1;
std::chrono::steady_clock::time_point cp_2;
std::chrono::steady_clock::time_point cp_3;
std::chrono::steady_clock::time_point cp_4;
std::chrono::steady_clock::time_point cp_5;
std::chrono::steady_clock::time_point cp_6;
std::chrono::steady_clock::time_point cp_7;
std::chrono::steady_clock::time_point cp_8;
- std::string file_test(filename);
+ /*std::string file_test(filename);
file_test = file_test.substr(0, file_test.size()-4);
std::string full_filename = file_test + "_clusters.csv";
std::ofstream f_coords;
f_coords.open(full_filename);
f_coords << "id,clusterNumber\n";
- f_coords.close();
+ f_coords.close();*/
//std::cout<<"on entering determine function";print_mem_usage();
double unf_time = 0.;
double A_time = 0.;
double invA_time = 0.;
double while_loop_time = 0.;
int min_cluster_size = 3;
set_lbb_params(filename);
int num_A_constructions = 0;
int num_invAb_ops = 0;
//std::cout<<"before to_build_list";print_mem_usage();
to_build_list.resize(numParticles());
//std::cout<<"after to_build_list";print_mem_usage();
calc_total_bonds();
//std::cout<<"before unique bonds";print_mem_usage();
unique_bonds.resize(total_num_bonds);
//std::cout<<"after unique bonds";print_mem_usage();
//std::cout<<"after init bullshit";print_mem_usage();
totalClusters_ = 0;
while (!check_if_particles_placed()) {
//cp_7 = std::chrono::steady_clock::now();
//cp_1 = std::chrono::steady_clock::now();
//std::cout<<"before unfolding";print_mem_usage();
init_lbb_unfolding_without_recursion();
//std::cout<<"after unfolding";print_mem_usage();
//init_lbb_unfolding();
//cp_2 = std::chrono::steady_clock::now();
//unf_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_2 - cp_1).count();
if (num_particles_for_cluster >= min_cluster_size) {
//std::cout<<"cluster size = "<<num_particles_for_cluster<<"\t num bonds = "<<num_bonds_for_cluster<<"\t time = "<<(unf_time * 1e-9)<<std::endl;
//cp_3 = std::chrono::steady_clock::now();
//std::cout<<"before init lbb matrices";print_mem_usage();
init_lbb_cluster_matrices();
//std::cout<<"after init lbb matrices";print_mem_usage();
//num_A_constructions++;
//cp_4 = std::chrono::steady_clock::now();
//A_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_4 - cp_3).count();
//cp_5 = std::chrono::steady_clock::now();
do {
//std::cout<<"before modify coords";print_mem_usage();
modify_coords_for_cluster();
//std::cout<<"after modify coords";print_mem_usage();
if (!A.isCompressed()){
A.makeCompressed();
//std::cout<<"test"<<std::endl;
}
//std::cout<<"before analyze pattern";print_mem_usage();
solver.analyzePattern(A);
solver.factorize(A);
//std::cout<<"after analyze pattern";print_mem_usage();
//std::cout<<"lbb = "<<n_lbb<<std::endl;
for (int axis = 0; axis < dim(); axis++){
b.setZero();
build_b_for_cluster(axis);
x = solver.solve(b);
num_invAb_ops++;
modify_coords_after_minimization_for_cluster(axis);
}
//std::cout<<"after solving pattern";print_mem_usage();
calculate_bond_lengths_for_cluster();
//print_minimized_coords(filename, n_lbb);
print_lbb_info();
} while(max_length > lbp_tolerance);
//cp_6 = std::chrono::steady_clock::now();
//invA_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_6 - cp_5).count();
}
- print_cluster_vals(filename);
+ //print_cluster_vals(filename);
totalClusters_++;
//cp_8 = std::chrono::steady_clock::now();
//while_loop_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_8 - cp_7).count();
//std::cout<<"total clusters = "<<totalClusters_<<std::endl;
}
//std::cout<<"while loop time = "<<(while_loop_time * 1e-9)<<std::endl;
fclose(f_lbp);
//std::cout<<"unfolding time = "<<unf_time * 1e-9<<std::endl;
//std::cout<<"A time = "<<A_time * 1e-9<<std::endl;
//std::cout<<"invA time = "<<invA_time * 1e-9<<std::endl;
//std::cout<<"n_lbb = "<<n_lbb<<std::endl;
//std::cout<<"num_A_constructions = "<<num_A_constructions<<std::endl;
//std::cout<<"num_invAb_ops = "<<num_invAb_ops<<std::endl;
//std::cout<<"unf 1 = "<<unf_time_1 * 1e-9<<std::endl;
//std::cout<<"unf 2 = "<<unf_time_2 * 1e-9<<std::endl;
//std::cout<<"A1 module = "<<time_A1<<std::endl;
//std::cout<<"A2 module = "<<time_A2<<std::endl;
}
-
-
void postprocessing::copy_positions_for_cluster()
{
int index;
for (int i = 0; i < num_particles_for_cluster; i++){
for (int axis = 0; axis < dim(); axis++){
index = index_to_particles[i];
modified_folded_x(i,axis) = pos(index,axis);
modified_folded_x(i,axis) -= L(axis) * round(modified_folded_x(i,axis)/L(axis));
}
}
}
void postprocessing::modify_coords_for_cluster()
{
int index;
for (int axis = 0; axis < dim(); axis++){
ref_pos(axis) = modified_folded_x(num_particles_for_cluster-1, axis);
}
for (int i = 0; i < num_particles_for_cluster; i++){
for (int axis = 0; axis < dim(); axis++){
modified_folded_x(i,axis) = modified_folded_x(i,axis) - ref_pos(axis);
modified_folded_x(i,axis) -= L(axis) * round(modified_folded_x(i,axis)/L(axis));
//std::cout<<pos(i,axis)<<" ";
}
//std::cout<<"\n";
}
}
void postprocessing::modify_coords_after_minimization_for_cluster(int axis)
{
for (int i = 0; i < num_particles_for_cluster-1; i++)
modified_folded_x(i,axis) = (L(axis) * x(i)) + ref_pos(axis);
modified_folded_x(num_particles_for_cluster-1,axis) = ref_pos(axis);
for (int i = 0; i < num_particles_for_cluster; i++)
modified_folded_x(i,axis) -= L(axis) * round(modified_folded_x(i,axis)/L(axis));
}
void postprocessing::put_particles_back_in_box(int axis)
{
for (int i = 0; i < num_particles_for_cluster; i++)
modified_folded_x(i,axis) -= L(axis) * round(modified_folded_x(i,axis)/L(axis));
}
void postprocessing::build_A_for_cluster()
{
int i,j;
//int index_i, index_j;
for (int index = 0; index < num_bonds_for_cluster; index++){
//cp_A1 = std::chrono::steady_clock::now();
i = unique_bonds[index].first;
j = unique_bonds[index].second;
//cp_A2 = std::chrono::steady_clock::now();
//time_A1 += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_A2 - cp_A1).count();
//index_i = index_to_particles[i];
//index_j = index_to_particles[j];
/*if (bond_map_status.count({index_i,index_j}) == 0){
std::cout<<"problem with A"<<std::endl;
}*/
//if (bond_map_status[{index_i,index_j}] == 1){
cp_A3 = std::chrono::steady_clock::now();
if (i == (num_particles_for_cluster-1)){
tripleList.push_back(Eigen::Triplet<double>(j,j,-1));
}
else if (j == (num_particles_for_cluster-1)){
tripleList.push_back(Eigen::Triplet<double>(i,i,-1));
}
else{
tripleList.push_back(Eigen::Triplet<double>(i,j,1));
tripleList.push_back(Eigen::Triplet<double>(j,i,1));
tripleList.push_back(Eigen::Triplet<double>(i,i,-1));
tripleList.push_back(Eigen::Triplet<double>(j,j,-1));
}
//cp_A4 = std::chrono::steady_clock::now();
//time_A2 += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_A4 - cp_A3).count();
//}
}
A.setFromTriplets(tripleList.begin(), tripleList.end());
tripleList.clear();
}
void postprocessing::build_b_for_cluster(int axis)
{
int i,j;
int index_i, index_j;
double dx;
for (int index = 0; index < num_bonds_for_cluster; index++){
i = unique_bonds[index].first;
j = unique_bonds[index].second;
index_i = index_to_particles[i];
index_j = index_to_particles[j];
/*if (bond_map_status.count({index_i,index_j}) == 0){
std::cout<<"problem with b"<<std::endl;
}*/
if (bond_map_status[{index_i,index_j}] == 1){
dx = round((modified_folded_x(i,axis) - modified_folded_x(j,axis))/L(axis));
if (i == (num_particles_for_cluster-1)){
b(j) += dx;
}
else if (j == (num_particles_for_cluster-1)){
b(i) -= dx;
}
else {
b(i) -= dx;
b(j) += dx;
}
}
}
}
void postprocessing::modify_A_matrix()
{
a_i = unique_bonds[max_row].first;
a_j = unique_bonds[max_row].second;
if (a_i == (num_particles_for_cluster-1)){
A.coeffRef(a_j,a_j) += 1;
}
else if (a_j == (num_particles_for_cluster-1)){
A.coeffRef(a_i, a_i) += 1;
}
else{
A.coeffRef(a_i,a_j) = 0;
A.coeffRef(a_j,a_i) = 0;
A.coeffRef(a_i,a_i) += 1;
A.coeffRef(a_j,a_j) += 1;
}
}
void postprocessing::copy_b()
{
for (int i = 0; i < b.size(); i++)
old_b(i) = b(i);
}
void postprocessing::calculate_b_diff()
{
b_diff = b - old_b;
if (b_diff.norm() > 1e-4)
std::cout<<"b vector changes"<<std::endl;
}
void postprocessing::increase_stiffness()
{
for (int i = 0; i < b.size(); i++){
b(i) = 1. * b(i);
}
}
void postprocessing::calculate_bond_lengths_direction_wise_for_cluster(int axis)
{
//bond_distance.setZero();
int i,j;
int index_i, index_j;
for (int index = 0; index < num_bonds_for_cluster; index++){
i = unique_bonds[index].first;
j = unique_bonds[index].second;
index_i = index_to_particles[i];
index_j = index_to_particles[j];
if (bond_map_status[{index_i,index_j}] == 1){
bond_lengths_direction_wise(index, axis) = (modified_folded_x(i, axis) - modified_folded_x(j, axis));
bond_lengths_direction_wise(index, axis) -= (L(axis) * round(bond_lengths_direction_wise(index, axis)/L(axis)));
bond_lengths_direction_wise(index, axis) *= bond_lengths_direction_wise(index, axis);
}
else {
bond_lengths_direction_wise(index, axis) = 0.;
}
}
}
void postprocessing::calculate_bond_lengths_for_cluster()
{
bond_lengths_direction_wise.setZero();
bond_lengths.setZero();
for (int axis = 0; axis < dim(); axis++){
calculate_bond_lengths_direction_wise_for_cluster(axis);
}
for (int i = 0; i < num_bonds_for_cluster; i++){
for (int axis = 0; axis < dim(); axis++) {
bond_lengths(i) += bond_lengths_direction_wise(i,axis);
}
bond_lengths(i) = sqrt(bond_lengths(i));
}
}
void postprocessing::calculate_bond_lengths_for_cluster(int axis)
{
//bond_distance.setZero();
int i,j;
int index_i, index_j;
for (int index = 0; index < num_bonds_for_cluster; index++){
i = unique_bonds[index].first;
j = unique_bonds[index].second;
index_i = index_to_particles[i];
index_j = index_to_particles[j];
if (bond_map_status[{index_i,index_j}] == 1){
bond_lengths(index) = (modified_folded_x(i, axis) - modified_folded_x(j, axis));
bond_lengths(index) -= (L(axis) * round(bond_lengths(index)/L(axis)));
bond_lengths(index) = abs(bond_lengths(index));
}
else {
bond_lengths(index) = 0.;
}
//std::cout<<index_i<<"-"<<index_j<<" = "<<bond_lengths[index]<<std::endl;
}
//std::cout<<"-----------------------------------------------------"<<std::endl;
}
void postprocessing::dump_xyz_stepwise(char *filename, int file_num)
{
std::string ext_filename;
ext_filename = filename;
ext_filename += "_" + std::to_string(file_num) + ".csv";
std::ofstream test_f;
test_f.open(ext_filename);
for (int axis = 0; axis < dim(); axis++){
if (axis != (dim()-1))
test_f <<"x"+std::to_string(axis)+",";
else
test_f <<"x"+std::to_string(axis)+"\n";
}
for (int i = 0; i < numParticles(); i++){
for (int axis = 0; axis < dim(); axis++){
if (axis != (dim()-1))
test_f <<std::to_string(modified_folded_x(i,axis))+",";
else
test_f <<std::to_string(modified_folded_x(i,axis))+"\n";
}
}
test_f.close();
}
}
\ No newline at end of file
diff --git a/postprocessing/src/lb_bonds_network.cpp b/postprocessing/src/lb_bonds_network.cpp
new file mode 100644
index 0000000..8a816b7
--- /dev/null
+++ b/postprocessing/src/lb_bonds_network.cpp
@@ -0,0 +1,237 @@
+#include <post.h>
+
+namespace post_p
+{
+
+void postprocessing::dump_LB_bonds_for_network(char *filename)
+{
+ determine_LB_bonds_network(filename);
+}
+
+void postprocessing::reset_unfolding_for_network()
+{
+ unf_checkpoint = 0;
+
+ for (int i = 0; i < numParticles(); i++){
+ is_placed(i) = false;
+ attachments_placed(i) = false;
+ }
+
+}
+
+void postprocessing::set_lbb_params_for_network(char *filename)
+{
+
+ FILE *f_network;
+
+ reset_unfolding_params();
+ build_bond_map();
+ reset_bond_map(true);
+
+ f_network = fopen(filename, "w");
+ fprintf(f_network,"bond,bond_length\n");
+ fclose(f_network);
+
+ lbp_tolerance = 1e-6;
+
+}
+
+void postprocessing::print_lbb_info_for_network(char *filename)
+{
+
+ FILE *f_network;
+
+ f_network = fopen(filename, "a");
+ max_length = bond_lengths.maxCoeff(&max_row);
+
+ if (max_length > lbp_tolerance){
+
+ switch_off_bonds({index_to_particles[unique_bonds[max_row].first], index_to_particles[unique_bonds[max_row].second]});
+ fprintf(f_network, "%d-%d,%lf\n",index_to_particles[unique_bonds[max_row].first],index_to_particles[unique_bonds[max_row].second],max_length);
+ n_lbb++;
+
+ }
+
+ fclose(f_network);
+
+}
+
+void postprocessing::build_A_for_network_cluster()
+{
+
+ int i,j;
+
+ for (int index = 0; index < num_bonds_for_cluster; index++){
+
+ i = unique_bonds[index].first;
+ j = unique_bonds[index].second;
+
+ tripleList.push_back(Eigen::Triplet<double>(i,j,1));
+ tripleList.push_back(Eigen::Triplet<double>(j,i,1));
+ tripleList.push_back(Eigen::Triplet<double>(i,i,-1));
+ tripleList.push_back(Eigen::Triplet<double>(j,j,-1));
+
+ }
+
+ A.setFromTriplets(tripleList.begin(), tripleList.end());
+ tripleList.clear();
+
+}
+
+void postprocessing::build_b_for_network_cluster()
+{
+ for (int i = 0; i < num_particles_for_cluster; i++)
+ {
+ b(i) = current_seed(index_to_particles[i]);
+ }
+}
+
+bool postprocessing::check_sink_and_source()
+{
+
+ bool sink = false;
+ bool source = false;
+
+ for (int i = 0; i < num_particles_for_cluster; i++){
+
+ if (current_seed(index_to_particles[i]) == 1)
+ sink = true;
+
+ if (current_seed(index_to_particles[i]) == -1)
+ source = true;
+
+ }
+
+ if (sink && source)
+ return true;
+
+ else
+ return false;
+
+}
+
+void postprocessing::fix_A_for_network_cluster()
+{
+
+ int i, j;
+
+ for (int i = 0; i < num_particles_for_cluster; i++){
+
+ if (current_seed(index_to_particles[i]) != 0)
+ A.coeffRef(i,i) -= 1;
+
+ }
+
+ //A.setFromTriplets(tripleList.begin(), tripleList.end());
+ //tripleList.clear();
+
+}
+
+void postprocessing::init_lbb_cluster_matrices_for_network_cluster()
+{
+ nnz = (2 * num_bonds_for_cluster) + num_particles_for_cluster;
+ A.reserve(nnz);
+ A.resize(num_particles_for_cluster, num_particles_for_cluster);
+
+ b.resize(num_particles_for_cluster);
+ x.resize(num_particles_for_cluster);
+ bond_lengths.resize(num_bonds_for_cluster);
+
+ A.setZero();
+ b.setZero();
+ x.setZero();
+ bond_lengths.setZero();
+
+ build_A_for_network_cluster();
+ build_b_for_network_cluster();
+
+ fix_A_for_network_cluster();
+
+ double temp;
+
+ /*for (int i = 0; i < num_particles_for_cluster; i++){
+ for (int j = 0; j < num_particles_for_cluster; j++){
+ std::cout<<A.coeff(i,j)<<"\t";
+ }
+ std::cout<<"\n";
+ }*/
+
+ //for (int i = 0; i < num_particles_for_cluster; i++)
+ //std::cout<<b.coeff(i)<<"\n";
+
+ //std::cout<<num_particles_for_cluster<<std::endl;
+ //std::cout<<num_bonds_for_cluster<<std::endl;
+
+
+}
+
+void postprocessing::calculate_bond_lengths_for_network_cluster()
+{
+
+ int i, j;
+
+ for (int index = 0; index < num_bonds_for_cluster; index++){
+
+ i = unique_bonds[index].first;
+ j = unique_bonds[index].second;
+
+ bond_lengths[index] = (x(i) - x(j)) * (x(i) - x(j));
+
+ }
+
+}
+
+void postprocessing::determine_LB_bonds_network(char *filename)
+{
+
+ set_lbb_params_for_network(filename);
+
+ to_build_list.resize(numParticles());
+
+ calc_total_bonds();
+ unique_bonds.resize(total_num_bonds);
+
+ totalClusters_ = 0;
+
+ while (!check_if_particles_placed()) {
+
+ init_lbb_unfolding_without_recursion();
+
+ //std::cout<<"num particles = "<<num_particles_for_cluster<<"\t num_bonds = "<<num_bonds_for_cluster<<"\n";
+
+ if (check_sink_and_source()) {
+
+ init_lbb_cluster_matrices_for_network_cluster();
+
+ do {
+
+ if (!A.isCompressed()){
+ A.makeCompressed();
+ }
+
+ solver.analyzePattern(A);
+ solver.factorize(A);
+
+ x = solver.solve(b);
+
+ calculate_bond_lengths_for_network_cluster();
+ print_lbb_info_for_network(filename);
+
+ if (max_length > lbp_tolerance){
+ reset_unfolding_for_network();
+ break;
+ }
+
+
+ } while(max_length > lbp_tolerance);
+
+
+ }
+
+ }
+
+}
+
+
+
+}
\ No newline at end of file
diff --git a/postprocessing/src/percolation_algo.cpp b/postprocessing/src/percolation_algo.cpp
index 1571b08..afc9775 100644
--- a/postprocessing/src/percolation_algo.cpp
+++ b/postprocessing/src/percolation_algo.cpp
@@ -1,879 +1,449 @@
#include <post.h>
namespace post_p
{
void postprocessing::build_bond_map()
{
for (int i = 0; i < numParticles(); i++){
for (int j = 0; j < numAttachments(i); j++)
bond_map_status.insert({{i,attachment(i,j)}, 1});
}
}
void postprocessing::switch_off_bonds(std::pair<int,int> bond)
{
bond_map_status[bond] = 0;
temp_pair.first = bond.second;
temp_pair.second = bond.first;
bond_map_status[temp_pair] = 0;
}
void postprocessing::switch_on_bonds(std::pair<int,int> bond)
{
bond_map_status[bond] = 1;
temp_pair.first = bond.second;
temp_pair.second = bond.first;
bond_map_status[temp_pair] = 1;
}
-void postprocessing::activate_path(std::vector<std::pair<int,int>> bonds, bool status)
-{
-
- for (int i = 0; i < bonds.size(); i++){
-
- if (status)
- switch_on_bonds(bonds[i]);
-
- else
- switch_off_bonds(bonds[i]);
-
-
- }
-
-}
-
void postprocessing::reset_bond_map(bool status)
{
if (status) {
for (int i = 0; i < numParticles(); i++){
for (int j = 0; j < numAttachments(i); j++) {
switch_on_bonds({i, attachment(i,j)});
}
}
}
else {
for (int i = 0; i < numParticles(); i++){
for (int j = 0; j < numAttachments(i); j++) {
switch_off_bonds({i, attachment(i,j)});
}
}
}
-}
-
-void postprocessing::create_subsystem()
-{
-
- reset_bond_map(false);
-
- for (int i = 0; i < current_path.size(); i++)
- switch_on_bonds(current_path[i]);
-
-
-}
-
-void postprocessing::erase_subsystem()
-{
-
- reset_bond_map(true);
-
- for (int i = 0; i < final_percolating_bonds.size(); i++){
-
- if (pb_status[i]) {
-
- for (int j = 0; j < final_percolating_bonds[i].size(); j++) {
-
- switch_off_bonds(final_percolating_bonds[i][j]);
-
- }
-
- }
- }
-
-
-}
-
-void postprocessing::switch_off_alt_lbp(int index)
-{
-
- for (int i = 0; i < final_percolating_bonds.size(); i++){
-
- if ((i != index) && (pb_status[i])) {
- //if (i != index) {
-
- for (int j = 0; j < final_percolating_bonds[i].size(); j++)
- switch_off_bonds(final_percolating_bonds[i][j]);
-
- }
-
- }
-
-
}
void postprocessing::reset_unfolding_params()
{
for (int j = 0; j < numParticles(); j++){
is_placed(j) = false;
attachments_placed(j) = false;
for (int axis = 0; axis < dim(); axis++)
unfolded_coords(j,axis) = 0.;
}
for (int axis = 0; axis < dim(); axis++)
temp_lbp_[axis] = 0;
}
-void postprocessing::print_lbp()
-{
-
- /*for (int i = 0; i < final_percolating_bonds.size(); i++){
-
- std::cout<<"lbp "<<i<<" = \t";
-
- for (int j = 0; j < final_percolating_bonds[i].size(); j++) {
-
- std::cout<<final_percolating_bonds[i][j].first<<"-"<<final_percolating_bonds[i][j].second<<"\t";
-
- }
-
- std::cout<<"\n";
-
- }*/
-
- for (int i = 0; i < weak_links.size(); i++){
-
- std::cout<<"lbp "<<i<<" = \t";
-
- for (int j = 0; j < weak_links[i].size(); j++) {
-
- std::cout<<weak_links[i][j].first<<"-"<<weak_links[i][j].second<<"\t";
-
- }
-
- std::cout<<"\n";
-
- }
-
-}
-
-void postprocessing::isolation_routine()
-{
-
- reset_unfolding_params();
-
- for (int i = 0; i < numParticles(); i++){
-
- if (is_placed(i) == false){
-
- path_percolation = false;
- path_components.clear();
- current_path.clear();
- unfold_for_lbp(i,i,true);
-
- if (path_percolation == true){
- create_subsystem();
- percolation_detection();
- path_percolation = true;
- total_lbp++;
- erase_subsystem();
- break;
- }
-
- }
-
-
- }
-
-}
-
void postprocessing::init_unfolding_for_lbp()
{
build_bond_map();
- //martin_test();
-
- /*bool build_pb = false;
- int ref_axis = 0;
- total_lbp = 0;
- path_percolation = false;
-
- for (int i = 0; i < numParticles(); i++){
- for (int axis = 0; axis < dim(); axis++){
- cluster_percolation(i,axis) = 0;
- }
- }
-
- do {
-
- isolation_routine();
-
- } while (path_percolation);
-
- postprocess_lbp();*/
- //print_lbp();
-
-}
-
-void postprocessing::percolation_detection()
-{
- percolating_bonds.clear();
- int size;
-
- for (int index = 0; index < current_path.size(); index++){
-
- reset_unfolding_params();
- switch_off_bonds(current_path[index]);
-
- path_percolation = false;
-
- for (int i = 0; i < path_components.size(); i++){
-
- if (is_placed(path_components[i]) == false){
-
- unfold_for_lbp(path_components[i],path_components[i],false);
-
- if (path_percolation == true)
- break;
-
-
- }
-
- }
-
- if (path_percolation == false)
- percolating_bonds.push_back(current_path[index]);
-
- switch_on_bonds(current_path[index]);
-
- }
-
- //pb_status[final_percolating_bonds.size()] = 1;
- //std::cout<<"size = "<<final_percolating_bonds.size()<<std::endl;
- final_percolating_bonds.push_back(percolating_bonds);
- pb_status.push_back(true);
-
-
-}
-
-void postprocessing::shred_path(int index)
-{
-
- reset_bond_map(true);
- switch_off_alt_lbp(index);
- switch_off_bonds(final_percolating_bonds[index][0]);
- isolation_routine();
-
- for (int i = 0; i < final_percolating_bonds[index].size(); i++)
- std::cout<<final_percolating_bonds[index][i].first<<"-"<<final_percolating_bonds[index][i].second<<",";
-
- std::cout<<"\n";
-
- index = final_percolating_bonds.size() - 1;
-
- for (int i = 0; i < final_percolating_bonds[index].size(); i++)
- std::cout<<final_percolating_bonds[index][i].first<<"-"<<final_percolating_bonds[index][i].second<<",";
-
- std::cout<<"\n";
-
- do {
-
- isolation_routine();
-
- } while (path_percolation);
-
-
-}
-
-void postprocessing::postprocess_lbp()
-{
-
- bool active_flag;
-
- for (int i = 0; i < final_percolating_bonds.size(); i++){
-
- //std::cout<<"size = "<<final_percolating_bonds.size()<<std::endl;
-
- reset_bond_map(true);
- percolating_bonds.clear();
- switch_off_alt_lbp(i);
- active_flag = true;
-
- for (int j = 0; j < final_percolating_bonds[i].size(); j++) {
-
- switch_off_bonds(final_percolating_bonds[i][j]);
- reset_unfolding_params();
- path_percolation = false;
-
- for (int k = 0; k < numParticles(); k++){
-
- if (is_placed(k) == false){
-
- if (path_percolation == true)
- break;
-
-
- }
-
- }
-
- if (path_percolation == false) {
- percolating_bonds.push_back(final_percolating_bonds[i][j]);
- active_flag = false;
- }
-
- switch_on_bonds(final_percolating_bonds[i][j]);
-
-
- }
-
- if (active_flag) {
- pb_status[i] = false;
- shred_path(i);
- }
-
- else {
- weak_links.push_back(percolating_bonds);
- }
-
- }
-
-
-}
-
-void postprocessing::unfold_for_lbp(const int prev, const int next, bool build_pb)
-{
-
- for (int axis = 0; axis < dim(); axis++)
- posDiff(axis) = get_periodic_image(pos(next,axis) - pos(prev,axis), axis);
-
- for (int axis = 0; axis < dim(); axis++){
- if (prev == next)
- unfolded_coords(next,axis) = pos(prev,axis) + posDiff(axis);
- else
- unfolded_coords(next,axis) = unfolded_coords(prev,axis) + posDiff(axis);
- }
-
- is_placed(next) = true;
-
- if (build_pb)
- path_components.push_back(next);
-
- for (int att = 0; att < numAttachments(next); att++)
- {
-
- //if (path_percolation)
- //break;
-
- temp_next = attachment(next, att);
-
- if (bond_map_status[{next, temp_next}] == 1){
-
- if ((build_pb == true)) {
-
- if(std::find(current_path.begin(), current_path.end(), std::make_pair(temp_next, next)) == current_path.end()) {
- current_path.push_back({next, temp_next});
- }
-
- }
-
- if (is_placed(temp_next) == false){
- unfoldedAttachments(next) = unfoldedAttachments(next) + 1;
- unfoldedAttachments(temp_next) = unfoldedAttachments(temp_next) + 1;
- unfold_for_lbp(next, temp_next, build_pb);
- }
-
- else {
-
- for (int axis = 0; axis < dim(); axis++) {
-
- posDiff(axis) = unfolded_coords(temp_next, axis) - unfolded_coords(next, axis);
-
- if ((posDiff(axis) > halfL(axis)) || (posDiff(axis) < -halfL(axis))){
- cluster_percolation(total_lbp,axis) = 1;
- temp_lbp_[axis]++;
-
- //if (axis == per_axis){
- path_percolation = true;
- //std::cout<<"axis = "<<axis<<"per_axis = "<<per_axis<<std::endl;
- //}
- }
-
- }
-
- }
-
- }
-
- }
}
void postprocessing::martin_test(char *lb_file)
{
reset_bond_map(true);
std::ifstream parser(lb_file, std::ifstream::in);
std::string str;
std::vector<std::string> results;
std::vector<std::pair<int,int>> lb_bonds;
getline(parser,str);
//std::vector<std::pair<int,int>> test_bonds;
while (getline(parser,str)){
results = split_string_by_delimiter(str, ',');
str = results[0];
//std::cout<<str<<std::endl;
results = split_string_by_delimiter(str, '-');
std::cout<<stoi(results[0])<<"\t"<<(results[1])<<std::endl;
lb_bonds_str.push_back(str);
lb_bonds.push_back({stoi(results[0]), stoi(results[1])});
}
parser.close();
for (int i = 0; i < lb_bonds.size(); i++)
switch_off_bonds(lb_bonds[i]);
path_percolation = false;
reset_unfolding_params();
int test_cluster_number = 0;
for (int i = 0; i < numParticles(); i++){
if (is_placed(i) == false){
//std::cout<<"i="<<i<<std::endl;
unfold_for_lbp(i, i, false);
test_cluster_number++;
}
}
std::cout<<"total clusters = "<<test_cluster_number<<std::endl;
if (path_percolation)
std::cout<<"still percolating"<<std::endl;
if (!path_percolation)
std::cout<<"no percolation at all"<<std::endl;
}
void postprocessing::makeCombiUtil(int n, int left, int k)
{
// Pushing this vector to a vector of vector
//std::cout<<"axis = "<<axis<<std::endl;
if (k == 0) {
if (path_percolation){
reset_bond_map(true);
for (int index = 0; index < tmp.size(); index++)
switch_off_bonds(unique_bonds[tmp[index]]);
path_percolation = false;
reset_unfolding_params();
for (int index = 0; index < numParticles(); index++){
if (is_placed(index) == false){
- unfold_for_lbp(index, index, false);
+ //to_build_list.clear();
+ num_particles_for_cluster = 0;
+ to_build_list[num_particles_for_cluster] = index;
+ num_particles_for_cluster++;
+ is_placed(index) = true;
+
+ //unfold_for_lbp(index, index, false);
+ unfold_for_lbp_without_recursion(index);
+
+ for (int iter_index = 1; iter_index < num_particles_for_cluster; iter_index++)
+ unfold_for_lbp_without_recursion(to_build_list[iter_index]);
+
}
+ if (path_percolation)
+ break;
+
}
if (!path_percolation){
for (int index = 0; index < tmp.size(); index++)
ans.push_back(tmp[index]);
}
}
return;
}
// i iterates from left to n. First time
// left will be 1
for (int i = left; i < n; ++i)
{
if (!path_percolation)
return;
tmp.push_back(i);
makeCombiUtil(n, i + 1, k - 1);
// Popping out last inserted element
// from the vector
tmp.pop_back();
}
}
// Prints all combinations of size k of numbers
// from 1 to n.
void postprocessing::makeCombi(int n, int k)
{
makeCombiUtil(n, 0, k);
}
void postprocessing::lbp_brute_force()
{
//makeCombi(5,3);
reset_bond_map(true);
+ to_build_list.resize(numParticles());
int num_bonds = 0;
int bond_sum = 0;
for (int i = 0; i < numParticles(); i++){
for (int j = 0; j < numAttachments(i); j++){
if (std::find(unique_bonds.begin(), unique_bonds.end(), std::make_pair(attachment(i,j),i)) == unique_bonds.end()) {
unique_bonds.push_back({i, attachment(i,j)});
num_bonds++;
}
}
}
path_percolation = true;
std::cout<<"bond"<<std::endl;
//for (int axis = 0; axis < dim(); axis++){
reset_bond_map(true);
path_percolation = true;
ans.clear();
- for (int i = 1; i < num_bonds; i++) {
+ for (int i = 0; i < num_bonds; i++) {
- makeCombi(num_bonds, i);
+ //std::cout<<"i = "<<i<<std::endl;
+ //makeCombi(num_bonds, i);
+ comb(num_bonds, i);
if (!path_percolation){
for (int k = 0; k < ans.size(); k++)
std::cout<<unique_bonds[ans[k]].first<<"-"<<unique_bonds[ans[k]].second<<"\n";
break;
}
}
//}
-
-
-}
-
}
-
-/*void postprocessing::percolation_detection()
+void postprocessing::comb(int N, int K)
{
+ std::string bitmask(K, 1); // K leading 1's
+ bitmask.resize(N, 0); // N-K trailing 0's
+
+ // print integers and permute bitmask
+ do {
- int lbp_min;
- int i_min;
- int bond_sum;
-
- for (int axis = 0; axis < dim(); axis++){
-
- lbp_min = load_bearing_paths(10,axis);
- i_min=10;
+ reset_bond_map(true);
+ ans.clear();
- /*for (int i = 0; i < numParticles(); i++){
-
- if (load_bearing_paths(i,axis) < lbp_min){
- lbp_min = load_bearing_paths(i,axis);
- i_min = i;
+ for (int i = 0; i < N; ++i) // [0..N-1] integers
+ {
+ if (bitmask[i]){
+ switch_off_bonds(unique_bonds[i]);
+ ans.push_back(i);
}
-
}
- if (lbp_min > 0) {
-
- for (int j = 0; j < numParticles(); j++){
- is_placed(j) = false;
- unfoldedAttachments(j) = 0;
- for (int axis_2 = 0; axis_2 < dim(); axis_2++)
- unfolded_coords(j,axis_2) = 0.;
- }
-
- for (int axis_2 = 0; axis_2 < dim(); axis_2++)
- temp_lbp_[axis_2] = 0;
-
- percolating_bonds.clear();
- unfold_for_lbp(i_min,i_min, true, axis);
-
- lbp_min = percolating_bonds.size();
-
- std::cout<<"lbp_min = "<<lbp_min<<std::endl;
-
- divisors = (int*)malloc(sizeof(int) * lbp_min);
- bond_bit_repr = (int*)malloc(sizeof(int) * lbp_min);
-
- for (int i = 0; i < lbp_min; i++)
- divisors[i] = std::pow(2, (lbp_min-1-i));
-
-
- //while (temp_lbp_[axis] != 0){
- for (int i = (std::pow(2,lbp_min)-1); i >= 0; i--) {
-
- for (int j = 0; j < lbp_min; j++)
- bond_bit_repr[j] = (i/divisors[j])%2;
-
- /*std::cout<<"i = "<<i<<"\t";
-
- for (int j = 0; j < lbp_min; j++)
- std::cout<<bond_bit_repr[j];
-
- std::cout<<"\n";
-
- for (int j = 0; j < lbp_min; j++) {
-
- bond_map_status[{percolating_bonds[j].first, percolating_bonds[j].second}] = bond_bit_repr[j];
- bond_map_status[{percolating_bonds[j].second, percolating_bonds[j].first}] = bond_bit_repr[j];
-
- }
-
- for (int j = 0; j < numParticles(); j++){
- is_placed(j) = false;
- unfoldedAttachments(j) = 0;
- for (int axis_2 = 0; axis_2 < dim(); axis_2++)
- unfolded_coords(j,axis_2) = 0.;
- }
+ path_percolation = false;
+ reset_unfolding_params();
- for (int axis_2 = 0; axis_2 < dim(); axis_2++)
- temp_lbp_[axis_2] = 0;
+ for (int index = 0; index < numParticles(); index++){
- unfold_for_lbp(i_min, i_min, false, axis);
+ if (is_placed(index) == false){
+ //to_build_list.clear();
+ num_particles_for_cluster = 0;
+ to_build_list[num_particles_for_cluster] = index;
+ num_particles_for_cluster++;
+ is_placed(index) = true;
+
+ //unfold_for_lbp(index, index, false);
+ unfold_for_lbp_without_recursion(index);
- /*for (int j = 0; j < numParticles(); j++){
+ for (int iter_index = 1; iter_index < num_particles_for_cluster; iter_index++)
+ unfold_for_lbp_without_recursion(to_build_list[iter_index]);
- if (is_placed(j) == false)
- std::cout<<"problem here"<<std::endl;
+ }
- }
+ if (path_percolation)
+ break;
+ }
- if (temp_lbp_[axis] == 0){
+ if (!path_percolation){
+ break;
+ }
- bond_sum = 0;
+ //std::cout << std::endl;
+ } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
- for (int j = 0; j < lbp_min; j++)
- bond_sum += bond_bit_repr[j];
+}
- std::cout<<"lbp = "<<(lbp_min-bond_sum)<<std::endl;
+void postprocessing::unfold_for_lbp_without_recursion(int i)
+{
+
+ for (int att = 0; att < numAttachments(i); att++)
+ {
- for (int j = 0; j < lbp_min; j++) {
+ index_j = attachment(i, att);
- bond_map_status[{percolating_bonds[j].first, percolating_bonds[j].second}] = 1;
- bond_map_status[{percolating_bonds[j].second, percolating_bonds[j].first}] = 1;
+ if (bond_map_status[{i, index_j}] == 1)
+ {
+ if (is_placed(index_j) == false)
+ {
+ is_placed(index_j) = true;
+ for (int axis = 0; axis < dim(); axis++){
+ posDiff(axis) = get_periodic_image(pos(index_j,axis) - pos(i,axis), axis);
+ unfolded_coords(index_j, axis) = unfolded_coords(i,axis) + posDiff(axis);
}
- for (int j = 0; j < percolating_bonds.size(); j++)
- std::cout<<"p = "<<percolating_bonds[j].first<<"\t"<<percolating_bonds[j].second<<std::endl;
-
+ to_build_list[num_particles_for_cluster] = index_j;
+ num_particles_for_cluster++;
+ }
- break;
+ else
+ {
- }
+ for (int axis = 0; axis < dim(); axis++){
- for (int j = 0; j < lbp_min; j++) {
+ posDiff(axis) = unfolded_coords(index_j, axis) - unfolded_coords(i,axis);
- bond_map_status[{percolating_bonds[j].first, percolating_bonds[j].second}] = 1;
- bond_map_status[{percolating_bonds[j].second, percolating_bonds[j].first}] = 1;
+ if ((posDiff(axis) > halfL(axis)) || (posDiff(axis) < -halfL(axis)))
+ path_percolation = true;
}
- /*std::cout<<"i = "<<i<<"\t";
-
- for (int j = 0; j < lbp_min; j++)
- std::cout<<bond_bit_repr[j];
-
- std::cout<<"\t axis = "<<axis<<"\t lbp = "<<temp_lbp_[axis]<<std::endl;
-
-
-
}
-
}
}
}
-
-void postprocessing::percolation_detection()
+void postprocessing::unfold_for_lbp(const int prev, const int next, bool build_pb)
{
- int lbp_min = 0;
- int i_min;
- int bond_sum;
- int lbp_min_axis;
-
- percolating_bonds.clear();
+ for (int axis = 0; axis < dim(); axis++)
+ posDiff(axis) = get_periodic_image(pos(next,axis) - pos(prev,axis), axis);
for (int axis = 0; axis < dim(); axis++){
-
- lbp_min = load_bearing_paths(0,axis);
- i_min=0;
-
- for (int i = 0; i < numParticles(); i++){
-
- if (load_bearing_paths(i,axis) < lbp_min){
- lbp_min = load_bearing_paths(i,axis);
- i_min = i;
- }
-
- }
-
- for (int j = 0; j < numParticles(); j++){
- is_placed(j) = false;
- unfoldedAttachments(j) = 0;
- for (int axis_2 = 0; axis_2 < dim(); axis_2++)
- unfolded_coords(j,axis_2) = 0.;
- }
-
- for (int axis_2 = 0; axis_2 < dim(); axis_2++)
- temp_lbp_[axis_2] = 0;
-
-
- unfold_for_lbp(i_min,i_min, true, axis);
-
+ if (prev == next)
+ unfolded_coords(next,axis) = pos(prev,axis) + posDiff(axis);
+ else
+ unfolded_coords(next,axis) = unfolded_coords(prev,axis) + posDiff(axis);
}
- lbp_min = percolating_bonds.size();
-
- //std::cout<<"lbp_min="<<lbp_min<<std::endl;
-
- if (lbp_min > 0) {
-
-
- divisors = (int*)malloc(sizeof(int) * lbp_min);
- bond_bit_repr = (int*)malloc(sizeof(int) * lbp_min);
-
- for (int i = 0; i < lbp_min; i++)
- divisors[i] = std::pow(2, (lbp_min-1-i));
-
-
- //while (temp_lbp_[axis] != 0){
- for (int i = (std::pow(2,lbp_min)-1); i >= 0; i--) {
+ is_placed(next) = true;
- for (int j = 0; j < lbp_min; j++)
- bond_bit_repr[j] = (i/divisors[j])%2;
+ if (build_pb)
+ path_components.push_back(next);
- /*std::cout<<"i = "<<i<<"\t";
+ for (int att = 0; att < numAttachments(next); att++)
+ {
- for (int j = 0; j < lbp_min; j++)
- std::cout<<bond_bit_repr[j];
+ //if (path_percolation)
+ //break;
- std::cout<<"\n";
+ temp_next = attachment(next, att);
- for (int j = 0; j < lbp_min; j++) {
+ if (bond_map_status[{next, temp_next}] == 1){
- bond_map_status[{percolating_bonds[j].first, percolating_bonds[j].second}] = bond_bit_repr[j];
- bond_map_status[{percolating_bonds[j].second, percolating_bonds[j].first}] = bond_bit_repr[j];
+ if ((build_pb == true)) {
- }
+ if(std::find(current_path.begin(), current_path.end(), std::make_pair(temp_next, next)) == current_path.end()) {
+ current_path.push_back({next, temp_next});
+ }
- for (int j = 0; j < numParticles(); j++){
- is_placed(j) = false;
- unfoldedAttachments(j) = 0;
- for (int axis_2 = 0; axis_2 < dim(); axis_2++)
- unfolded_coords(j,axis_2) = 0.;
}
-
- for (int axis_2 = 0; axis_2 < dim(); axis_2++)
- temp_lbp_[axis_2] = 0;
-
- unfold_for_lbp(i_min, i_min, false, 0);
-
- for (int j = 0; j < numParticles(); j++){
-
- if (is_placed(j) == false)
- std::cout<<"problem here"<<std::endl;
-
+
+ if (is_placed(temp_next) == false){
+ unfoldedAttachments(next) = unfoldedAttachments(next) + 1;
+ unfoldedAttachments(temp_next) = unfoldedAttachments(temp_next) + 1;
+ unfold_for_lbp(next, temp_next, build_pb);
}
- lbp_min_axis = 0;
-
- for (int axis = 0; axis < dim(); axis++)
- lbp_min_axis += temp_lbp_[axis];
-
- if (lbp_min_axis == 0){
-
- bond_sum = 0;
-
- for (int j = 0; j < lbp_min; j++)
- bond_sum += bond_bit_repr[j];
+ else {
+
+ for (int axis = 0; axis < dim(); axis++) {
- std::cout<<"lbp = "<<(lbp_min-bond_sum)<<std::endl;
+ posDiff(axis) = unfolded_coords(temp_next, axis) - unfolded_coords(next, axis);
- for (int j = 0; j < lbp_min; j++) {
+ if ((posDiff(axis) > halfL(axis)) || (posDiff(axis) < -halfL(axis))){
+ cluster_percolation(total_lbp,axis) = 1;
+ temp_lbp_[axis]++;
- bond_map_status[{percolating_bonds[j].first, percolating_bonds[j].second}] = 1;
- bond_map_status[{percolating_bonds[j].second, percolating_bonds[j].first}] = 1;
+ //if (axis == per_axis){
+ path_percolation = true;
+ //std::cout<<"axis = "<<axis<<"per_axis = "<<per_axis<<std::endl;
+ //}
+ }
}
- break;
-
}
-
}
}
-
-
-}*/
+}
+}
\ No newline at end of file
diff --git a/postprocessing/target/lb_bonds_clusterwise_invA.cpp b/postprocessing/target/lb_bonds_clusterwise_invA.cpp
index 8858726..3bfd5a9 100644
--- a/postprocessing/target/lb_bonds_clusterwise_invA.cpp
+++ b/postprocessing/target/lb_bonds_clusterwise_invA.cpp
@@ -1,40 +1,42 @@
#include <post.h>
#include <ctime>
int main(int argc, char *argv[])
{
/*std::chrono::steady_clock::time_point cp_1;
std::chrono::steady_clock::time_point cp_2;
std::chrono::steady_clock::time_point cp_3;
std::chrono::steady_clock::time_point cp_4;
double total_time = 0.;
double read_time = 0.;
cp_1 = std::chrono::steady_clock::now();
cp_3 = std::chrono::steady_clock::now();*/
post_p::postprocessing *test = new post_p::postprocessing(argv[1]);
//cp_4 = std::chrono::steady_clock::now();
//read_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_4 - cp_3).count();
//std::cout<<"read time = "<<read_time * 1e-9<<std::endl;
FILE *f_time;
if (argc < 3){
std::cout<<"please provide config_filename results_filename"<<std::endl;
return 0;
}
- else
- test->dump_lb_bonds_for_cluster_via_invA(argv[2]);
+ else{
+ //test->dump_lb_bonds_for_cluster_via_invA(argv[2]);
+ test->dump_LB_bonds_for_network(argv[2]);
+ }
//cp_2 = std::chrono::steady_clock::now();
//total_time += std::chrono::duration_cast<std::chrono::nanoseconds>(cp_2 - cp_1).count();
//std::cout<<"total time = "<<total_time * 1e-9<<std::endl;
delete test;
return 0;
}
\ No newline at end of file
diff --git a/scripts/setup_simulation_rsp_2d.py b/scripts/setup_simulation_rsp_2d.py
index 7891c7c..d6e06bf 100755
--- a/scripts/setup_simulation_rsp_2d.py
+++ b/scripts/setup_simulation_rsp_2d.py
@@ -1,175 +1,175 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 4 13:56:09 2021
@author: samarth
"""
import pandas as pd
import os
import sys
from pathlib import Path
from joblib import Parallel, delayed
from matplotlib import pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from random import random, seed
def create_dirs(phi, N, seed_pct, dim):
global D
D = dim
global params_dir
phi_str = format(phi, '.8f')
params_dir = './'+str(D)+'d/params/phi='+phi_str+'/seed_pct='+str(seed_pct)+'/L='+str(N)+'/'
path = Path(params_dir)
try:
path.mkdir(parents=True, exist_ok=True)
except OSError:
pass
else:
pass
global results_dir
results_dir = './'+str(D)+'d/results/phi='+phi_str+'/seed_pct='+str(seed_pct)+'/L='+str(N)+'/'
path = Path(results_dir)
try:
path.mkdir(parents=True, exist_ok=True)
except OSError:
pass
else:
pass
global config_files_dir
config_files_dir = './'+str(D)+'d/config_files/phi='+phi_str+'/seed_pct='+str(seed_pct)+'/L='+str(N)+'/'
path = Path(config_files_dir)
try:
path.mkdir(parents=True, exist_ok=True)
except OSError:
pass
else:
pass
def run_simulation(phi, L, seed_pct, rng_seed):
lattice=1
params_file = params_dir+'params_'+str(rng_seed)+'.csv'
file = open(params_file, "w")
file.write("system=random_site_percolation\n")
file.write("lattice="+str(lattice)+"\n")
file.write("D="+str(D)+"\n")
file.write("phi="+str(phi)+"\n")
file.write("seed_pct="+str(seed_pct)+"\n")
for axis in range(D):
file.write("x"+str(axis)+"_bc=periodic\n")
for axis in range(D):
file.write("x"+str(axis)+"_L="+str(L)+"\n")
file.write("rng_seed="+str(rng_seed)+"\n")
file.close()
print("phi = {}, L={},rng_seed={} simulation".format(phi,L,rng_seed))
config_filename=config_files_dir+str(rng_seed)+'.csv'
command = 'simulation '+params_file+' '+config_filename
os.system(command)
#print(command)
print("phi = {}, L={},rng_seed={} lbp".format(phi,L,rng_seed))
rog_filename=results_dir+'lbp_'+str(rng_seed)+'.csv'
command = 'lb_bonds_clusterwise_invA ' + config_filename +' ' + rog_filename
os.system(command)
command = 'rm ' + config_filename
os.system(command)
def create_dirs_for_mk_table(dim):
global D
D = dim
global params_dir
params_dir = './mk_table/params/'
path = Path(params_dir)
try:
path.mkdir(parents=True, exist_ok=True)
except OSError:
pass
else:
pass
global results_dir
results_dir = './mk_table/results/'
path = Path(results_dir)
try:
path.mkdir(parents=True, exist_ok=True)
except OSError:
pass
else:
pass
global config_files_dir
config_files_dir = './mk_table/config_files/'
path = Path(config_files_dir)
try:
path.mkdir(parents=True, exist_ok=True)
except OSError:
pass
else:
pass
#phi=0.15
seeds=100
seed_pct=100
dim = 2
L_range = [50, 75, 100, 200, 500, 1000, 2000]
phi_min = 0.005
phi_max = 0.05
phi_c = 0.5927462
phi_min += phi_c
phi_max += phi_c
total_steps = 15
d_phi = (np.log(phi_max/phi_min))/(total_steps)
phi_range = []
for i in range(total_steps):
temp_phi = phi_min * np.exp(i * (d_phi))
phi_range.append(temp_phi)
-for i in phi_range:
- for L in L_range:
+for L in L_range:
+ for i in phi_range:
create_dirs(i, L, seed_pct, dim)
Parallel(n_jobs=10)(delayed(run_simulation)(i, L, seed_pct, rng_seed) for rng_seed in range(seeds))
diff --git a/simulation_files/bin/simulation b/simulation_files/bin/simulation
new file mode 100755
index 0000000..e36369f
Binary files /dev/null and b/simulation_files/bin/simulation differ
diff --git a/simulation_files/obj/brownian_movement.o b/simulation_files/obj/brownian_movement.o
new file mode 100644
index 0000000..cf29019
Binary files /dev/null and b/simulation_files/obj/brownian_movement.o differ
diff --git a/simulation_files/obj/brownian_movement_offlattice.o b/simulation_files/obj/brownian_movement_offlattice.o
new file mode 100644
index 0000000..e53b643
Binary files /dev/null and b/simulation_files/obj/brownian_movement_offlattice.o differ
diff --git a/simulation_files/obj/check_aggregation.o b/simulation_files/obj/check_aggregation.o
new file mode 100644
index 0000000..0e36cf2
Binary files /dev/null and b/simulation_files/obj/check_aggregation.o differ
diff --git a/simulation_files/obj/check_aggregation_offlattice.o b/simulation_files/obj/check_aggregation_offlattice.o
new file mode 100644
index 0000000..20205a3
Binary files /dev/null and b/simulation_files/obj/check_aggregation_offlattice.o differ
diff --git a/simulation_files/obj/cluster.o b/simulation_files/obj/cluster.o
new file mode 100644
index 0000000..4b86a08
Binary files /dev/null and b/simulation_files/obj/cluster.o differ
diff --git a/simulation_files/obj/dlma_iterator.o b/simulation_files/obj/dlma_iterator.o
new file mode 100644
index 0000000..d3ebabf
Binary files /dev/null and b/simulation_files/obj/dlma_iterator.o differ
diff --git a/simulation_files/obj/dlma_save_config.o b/simulation_files/obj/dlma_save_config.o
new file mode 100644
index 0000000..2de7174
Binary files /dev/null and b/simulation_files/obj/dlma_save_config.o differ
diff --git a/simulation_files/obj/dlma_system.o b/simulation_files/obj/dlma_system.o
new file mode 100644
index 0000000..e4ce78a
Binary files /dev/null and b/simulation_files/obj/dlma_system.o differ
diff --git a/simulation_files/obj/dlma_system_offlattice.o b/simulation_files/obj/dlma_system_offlattice.o
new file mode 100644
index 0000000..eea821a
Binary files /dev/null and b/simulation_files/obj/dlma_system_offlattice.o differ
diff --git a/simulation_files/obj/dlma_system_onlattice.o b/simulation_files/obj/dlma_system_onlattice.o
new file mode 100644
index 0000000..ac17b55
Binary files /dev/null and b/simulation_files/obj/dlma_system_onlattice.o differ
diff --git a/simulation_files/obj/iterator_factory.o b/simulation_files/obj/iterator_factory.o
new file mode 100644
index 0000000..921fc8e
Binary files /dev/null and b/simulation_files/obj/iterator_factory.o differ
diff --git a/simulation_files/obj/mass_aggregation_condition.o b/simulation_files/obj/mass_aggregation_condition.o
new file mode 100644
index 0000000..55aed14
Binary files /dev/null and b/simulation_files/obj/mass_aggregation_condition.o differ
diff --git a/simulation_files/obj/normal_bind.o b/simulation_files/obj/normal_bind.o
new file mode 100644
index 0000000..eb99c32
Binary files /dev/null and b/simulation_files/obj/normal_bind.o differ
diff --git a/simulation_files/obj/off_lattice.o b/simulation_files/obj/off_lattice.o
new file mode 100644
index 0000000..0c599e5
Binary files /dev/null and b/simulation_files/obj/off_lattice.o differ
diff --git a/simulation_files/obj/on_lattice.o b/simulation_files/obj/on_lattice.o
new file mode 100644
index 0000000..21bfb20
Binary files /dev/null and b/simulation_files/obj/on_lattice.o differ
diff --git a/simulation_files/obj/particle.o b/simulation_files/obj/particle.o
new file mode 100644
index 0000000..9d493b2
Binary files /dev/null and b/simulation_files/obj/particle.o differ
diff --git a/simulation_files/obj/periodic.o b/simulation_files/obj/periodic.o
new file mode 100644
index 0000000..f5e6766
Binary files /dev/null and b/simulation_files/obj/periodic.o differ
diff --git a/simulation_files/obj/split_string.o b/simulation_files/obj/split_string.o
new file mode 100644
index 0000000..71b66c7
Binary files /dev/null and b/simulation_files/obj/split_string.o differ
diff --git a/simulation_files/obj/system_factory.o b/simulation_files/obj/system_factory.o
new file mode 100644
index 0000000..ffd89e5
Binary files /dev/null and b/simulation_files/obj/system_factory.o differ
diff --git a/simulation_files/simulation.o b/simulation_files/simulation.o
new file mode 100644
index 0000000..d8510c5
Binary files /dev/null and b/simulation_files/simulation.o differ