Skip to content
Snippets Groups Projects

Mrk br0

Closed Mrk br0
mrk_br0 into main
Closed Masudur Rahaman Kazi requested to merge mrk_br0 into main
Compare
and
6 files
+ 426
− 118
Compare changes
  • Side-by-side
  • Inline
Files
6
ex02/src/cell_particle.cpp
+ 198
− 50
// ConsoleApplication1.cpp : This file contains the 'main' function. Program execution begins and ends there.
//
#include "particle_class.h"
#include <iostream>
#include <fstream>
#include <string>
@@ -9,22 +9,105 @@
#include <cmath>
using namespace std;
//class Cell {
//public:
// int value;
// Cell* next;
//};
class Particle {
class Cell {
public:
double v[3] = {0.0, 0.0, 0.0}; //velocity set with gaussian distribution
double a[3] = {0.0, 0.0, 0.0};
double p[3] = {0.0, 0.0, 0.0}; //position of the nodes of the 3D matrix
int cell_ptr = -2;//stores the particle-wise value for the linked list of particles belonging to same cell
int particle_ptr = -1;
int grid_pos[3] = {0, 0, 0};
vector<int> neighbours;
vector<int> findCellNeighbour(const int grid_pos[3], int cell_idx, int grid_edge) {
vector<int> neighbours;
if (grid_pos[0] > 0)
neighbours.push_back(cell_idx - 1);
else
neighbours.push_back(cell_idx+grid_edge-1);
if (grid_pos[0] < grid_edge - 1)
neighbours.push_back(cell_idx + 1);
else
neighbours.push_back(cell_idx-(grid_edge-1));
for (int i = neighbours.size() - 1; i >= 0; i--) {
if (grid_pos[1] > 0)
neighbours.push_back(neighbours[i] - grid_edge);
else
neighbours.push_back(neighbours[i] + grid_edge*(grid_edge-1));
if (grid_pos[1] < grid_edge - 1)
neighbours.push_back(neighbours[i] + grid_edge);
else
neighbours.push_back(neighbours[i] - grid_edge*(grid_edge-1));
}
if (grid_pos[1] > 0)
neighbours.push_back(cell_idx - grid_edge);
else
neighbours.push_back(cell_idx + grid_edge*(grid_edge-1));
if (grid_pos[1] < grid_edge - 1)
neighbours.push_back(cell_idx + grid_edge);
else
neighbours.push_back(cell_idx - grid_edge*(grid_edge-1));
for (int i = neighbours.size() - 1; i >= 0; i--) {
if (grid_pos[2] > 0)
neighbours.push_back(neighbours[i] - grid_edge * grid_edge);
else
neighbours.push_back(neighbours[i] + grid_edge*grid_edge*(grid_edge-1));
if (grid_pos[2] < grid_edge - 1)
neighbours.push_back(neighbours[i] + grid_edge * grid_edge);
else
neighbours.push_back(neighbours[i] - grid_edge*grid_edge*(grid_edge-1));
}
if (grid_pos[2] > 0)
neighbours.push_back(cell_idx - grid_edge * grid_edge);
else
neighbours.push_back(cell_idx + grid_edge*grid_edge*(grid_edge-1));
if (grid_pos[2] < grid_edge - 1)
neighbours.push_back(cell_idx + grid_edge * grid_edge);
else
neighbours.push_back(cell_idx - grid_edge*grid_edge*(grid_edge-1));
return neighbours;
}
void addParticle(int cellIdx,int particleIdx, Particle* particle){
particle->setCellPtr(this->particle_ptr);
particle->setCellIdx(cellIdx);
this->particle_ptr = particleIdx;
}
void removeParticle(int cellIdx,int particleIdx, Particle** particles){
if(this->particle_ptr==particleIdx){
this->particle_ptr=particles[particleIdx]->getCellPtr();
}
else{
int currParticleIdx=this->particle_ptr;
while (particles[currParticleIdx]->getCellPtr()!=-1){
if(particles[currParticleIdx]->getCellPtr()==particleIdx){
particles[currParticleIdx]->setCellPtr(particles[particleIdx]->getCellPtr());
break;
}
currParticleIdx=particles[currParticleIdx]->getCellPtr();
}
}
if(particles[particleIdx]->getCellIdx()==cellIdx){ //if condition to avoid changing cell index after Particle::updateCellIndex has been called
particles[particleIdx]->setCellIdx(-1);
}
}
vector<int> findParticleNeighbour(int particleIdx, Particle** particles){
vector<int> neighbours;
int flag=0;
int currParticleIdx=this->particle_ptr;
if (currParticleIdx!=particleIdx)
neighbours.push_back(currParticleIdx);
else
flag=1;
currParticleIdx=particles[currParticleIdx]->getCellPtr();
while (currParticleIdx!=-1){
if (currParticleIdx!=particleIdx)
neighbours.push_back(currParticleIdx);
else
flag=1;
currParticleIdx=particles[currParticleIdx]->getCellPtr();
}
if (flag==0)
neighbours.clear();
return neighbours;
}
};
vector<string> split(const string &str, char delimiter) {
@@ -38,66 +121,82 @@ vector<string> split(const string &str, char delimiter) {
return res;
}
// __global__ void particle_kernel(Cell** cells,int N_cells,Particle** particles,int N_particles){
// }
int main() {
int N_cells = 0; //a perfect cube number for symmetric 3d distribution
int N_particles = 0; //number of particles per cell
double Lsys = 0.0;
double T = 0.0;
double Lcell=0.0;
int *cells = nullptr; // Usage of nullptr maybe better than NULL in c++
Cell **cells = nullptr;
Particle **particles = nullptr;
fstream file;
file.open("cell_contents.txt",ios::in); //open a file to perform read operation using file object
file.open("cell_contents.txt", ios::in); //open a file to perform read operation using file object
if (file.is_open()) { //checking whether the file is open
cout << "File open" << endl;
vector<string> vec;
string line;
getline(file, line);
cout << "First line: " << line << endl;
vec = split(line, ' ');
vector<string> vec = split(line, ' ');
// Reading initial data from file
N_cells = stoi(vec[0]);
N_particles = stoi(vec[1]);
Lsys = stod(vec[2]);
T = stod(vec[3]) + 273.0;
Lcell=Lsys/int(cbrt(N_cells));
cout << "N_cells = " << N_cells << endl;
cout << "N_particles per cell= " << N_particles << endl;
cells = new int[N_cells];
for (int i = 0; i < N_cells; ++i) {
cells[i] = -1;
}
cells = new Cell *[N_cells];
// for (int i = 0; i < N_cells; ++i)
// {
// cells[i]->particle_ptr = -1;
// }
particles = new Particle *[N_particles * N_cells];
double origin[3] = {0.0, 0.0, 0.0};
for (int i = 0;
i < N_cells; ++i)// placing molecules at equidistant location in cubical grid of side (Lsys*cbrt(N_cell))
{
int cell_idx=i;
int N_gridAxis=int(cbrt(N_cells));
int N_cellAxis=int(cbrt(N_particles));
origin[0]=(cell_idx%N_gridAxis)*Lsys+Lsys/(N_cellAxis*2);
cell_idx/=N_gridAxis;
origin[1]=(cell_idx%N_gridAxis)*Lsys+Lsys/(N_cellAxis*2);
cell_idx/=N_gridAxis;
origin[2]=(cell_idx%N_gridAxis)*Lsys+Lsys/(N_cellAxis*2);
for (int j = 0; j < N_particles; ++j)
// placing molecules at equidistant location in cubical grid of side (Lsys * length of each cell)
// TODO isn't Lsys already the entire length of simulation?
// TODO: what is interaction range and is Lsys at least twice the size
for (int i = 0; i < N_cells; ++i)
{
int particle_idx=j;
double x=origin[0]+(particle_idx%N_cellAxis)*Lsys/N_cellAxis;
particle_idx/=N_cellAxis;
double y=origin[1]+(particle_idx%N_cellAxis)*Lsys/N_cellAxis;
particle_idx/=N_cellAxis;
double z=origin[2]+(particle_idx%N_cellAxis)*Lsys/N_cellAxis;
particle_idx=i*N_particles+j;
particles[particle_idx] = new Particle();
particles[particle_idx]->p[0]=x;
particles[particle_idx]->p[1]=y;
particles[particle_idx]->p[2]=z;
particles[particle_idx]->cell_ptr = cells[i];
cells[i] = particle_idx;
}
int cell_idx = i;
int N_gridAxis = int(cbrt(N_cells)); //TODO: I think they should be called N_gridLen
int N_cellAxis = int(cbrt(N_particles));
cells[i] = new Cell();
origin[0] = (cell_idx % N_gridAxis) * Lcell + Lcell / (N_cellAxis * 2);
cells[i]->grid_pos[0] = cell_idx % N_gridAxis;
cell_idx /= N_gridAxis;
origin[1] = (cell_idx % N_gridAxis) * Lcell + Lcell / (N_cellAxis * 2);
cells[i]->grid_pos[1] = cell_idx % N_gridAxis;
cell_idx /= N_gridAxis;
origin[2] = (cell_idx % N_gridAxis) * Lcell + Lcell / (N_cellAxis * 2);
cells[i]->grid_pos[2] = cell_idx % N_gridAxis;
cells[i]->neighbours = cells[i]->findCellNeighbour(cells[i]->grid_pos, i, N_gridAxis);
for (int j = 0; j < N_particles; ++j) {
int particle_idx = j;
double x = origin[0] + (particle_idx % N_cellAxis) * Lcell / N_cellAxis;
particle_idx /= N_cellAxis;
double y = origin[1] + (particle_idx % N_cellAxis) * Lcell / N_cellAxis;
particle_idx /= N_cellAxis;
double z = origin[2] + (particle_idx % N_cellAxis) * Lcell / N_cellAxis;
particle_idx = i * N_particles + j;
particles[particle_idx] = new Particle(1, {x, y, z});
cells[i]->addParticle(i,particle_idx,particles[particle_idx]);
// particles[particle_idx]->setCellPtr(cells[i]->particle_ptr);
// particles[particle_idx]->setCellIdx(i);
// cells[i]->particle_ptr = particle_idx;
}
}
// for (int i = 0; i < N_cells; ++i) { //reading cell-particle relations from file and implementing the linked list algorithm
@@ -115,16 +214,65 @@ int main() {
}
cout << "Cells:" << endl;
for (int i = 0; i < N_cells; ++i) {
cout << cells[i] << " ";
cout << "Cell " << i << ":" << cells[i]->particle_ptr << " [";
for (int j = 0; j < cells[i]->neighbours.size(); ++j) {
cout << cells[i]->neighbours[j] << ",";
}
cout << "]" << endl;
}
cout << endl;
cout << "Particles:" << endl;
for (int i = 0; i < N_particles * N_cells; ++i) {
cout << particles[i]->cell_ptr << " ";
cout << "[" << particles[i]->p[0] << "," << particles[i]->p[1] << "," << particles[i]->p[2] << "]" << endl;
cout << "Particle " << i << ": ";
cout << particles[i]->getCellPtr() << " ";
cout << particles[i]->getCellIdx() << endl;
//cout << particles[i]->updateCellIndex(Lsys,N_cells,N_particles) << " ";
// cout << "[" << particles[i]->getP()[0] << "," << particles[i]->getP()[1] << "," << particles[i]->getP()[2]
// << "]" << endl;
}
cout << endl;
//--------EXAMPLE implementation of Particle::updateCellIndex---------
// double x[]={Lsys+0.5,Lsys+0.5,Lsys+0.5};
// int particleIdx=N_cells*N_particles-8;
// particles[particleIdx]->setP(x);
// int oldCellIdx=particles[particleIdx]->updateCellIndex(Lsys,N_cells,N_particles);
// int newCellIdx=particles[particleIdx]->getCellIdx();
// cout << oldCellIdx<<" "<<newCellIdx<<endl;
// cells[oldCellIdx]->removeParticle(oldCellIdx,particleIdx,particles);
// cells[newCellIdx]->addParticle(newCellIdx,particleIdx,particles[particleIdx]);
// cout << "Particle " << particleIdx << ": ";
// cout << particles[particleIdx]->getCellPtr() << " ";
// cout << particles[particleIdx]->getCellIdx() << " ";
// cout << "[" << particles[particleIdx]->getP()[0] << "," << particles[particleIdx]->getP()[1] << "," << particles[particleIdx]->getP()[2]
// << "]" << endl;
// cout << "Cell " << oldCellIdx << ":" << cells[oldCellIdx]->particle_ptr <<endl;
// int currParticleIdx=cells[oldCellIdx]->particle_ptr;
// while (particles[currParticleIdx]->getCellPtr()!=-1){
// cout << "Particle " << currParticleIdx << ": ";
// cout << particles[currParticleIdx]->getCellPtr() << " ";
// currParticleIdx=particles[currParticleIdx]->getCellPtr();
// }
// cout << "Particle " << currParticleIdx << ": ";
// cout << particles[currParticleIdx]->getCellPtr() << endl;
// cout << "Cell " << newCellIdx << ":" << cells[newCellIdx]->particle_ptr <<endl;
// currParticleIdx=cells[newCellIdx]->particle_ptr;
// while (particles[currParticleIdx]->getCellPtr()!=-1){
// cout << "Particle " << currParticleIdx << ": ";
// cout << particles[currParticleIdx]->getCellPtr() << " ";
// currParticleIdx=particles[currParticleIdx]->getCellPtr();
// }
// cout << "Particle " << currParticleIdx << ": ";
// cout << particles[currParticleIdx]->getCellPtr() << endl;
//---------------------------------------------------------------------------------
// --------EXAMPLE implementation of Cell::findParticleNeighbour---------
// vector<int> neighbours=cells[particles[26]->getCellIdx()]->findParticleNeighbour(26,particles);
// cout<<neighbours.size()<<endl;
// for(int i : neighbours)
// cout<<i<<endl;
// ----------------------------------------------------------------------
delete[] cells;
for (int i = 0; i < N_particles * N_cells; ++i) {
delete particles[i];