calc_concent.cpp

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <cassert>
#include <cmath>
#include <cstdlib>
 
#define INLINE __attribute__((always_inline))
 
template<typename T>
struct vector3{
  T x,y,z;
  
  //constructors
  INLINE vector3<T>(void){
    x = y = z = (T)0.;
  }
  template<typename U>
  INLINE vector3<T>(const vector3<U> &rhs) : x(rhs.x), y(rhs.y), z(rhs.z) {}
  template<typename U>
  INLINE vector3<T>(const U &r){
    x = y = z = (T)r;
  }
  template<typename U>
  INLINE vector3<T>(const U *p){
    x = p[0];
    y = p[1];
    z = p[2];
  }
 
  //operator overloads
  INLINE const vector3<T> operator+(const vector3<T>& obj) const {
    vector3<T> ret;
    ret.x = x + obj.x;
    ret.y = y + obj.y;
    ret.z = z + obj.z;
    return ret;
  }
  INLINE const vector3<T> operator-(const vector3<T>& obj) const {
    vector3<T> ret;
    ret.x = x - obj.x;
    ret.y = y - obj.y;
    ret.z = z - obj.z;
    return ret;
  }
  INLINE const vector3<T>& operator+=(const vector3<T>& obj){
    x += obj.x;
    y += obj.y;
    z += obj.z;
    return *this;
  }
  INLINE const vector3<T>& operator-=(const vector3<T>& obj){
    x -= obj.x;
    y -= obj.y;
    z -= obj.z;
    return *this;
  }
  INLINE T operator*(const vector3<T>& obj) const {
    const T prod = x*obj.x + y*obj.y + z*obj.z;
    return prod;
  }
  INLINE vector3<T> operator*(const T& cf) const {
    const T cf_ = (T)cf;
    vector3<T> ret;
    ret.x = x * cf_;
    ret.y = y * cf_;
    ret.z = z * cf_;
    return ret;
  }
  template<typename U>
  INLINE vector3<T>& operator*=(const U& cf){
    const T cf_ = (T)cf;
    x *= cf_;
    y *= cf_;
    z *= cf_;
    return *this;
  }
  template<typename U>
  INLINE vector3<T>& operator/=(const U& cf){
    const T cf_ = (T)cf;
    x /= cf_;
    y /= cf_;
    z /= cf_;
    return *this;
  }
 
  //inverse square root
  //for IEEE754
  INLINE float m_rsqrt(const float& r) const {
    const float rHalf = 0.5 * r;
    const int   tmp   = 0x5F3759DF - ( *(int*)&r >> 1 );
    float rRes        = *(float*)&tmp;
    rRes *= ( 1.5 - ( rHalf * rRes * rRes ) );
    rRes *= ( 1.5 - ( rHalf * rRes * rRes ) );
    return rRes;
  }
  //for IEEE754
  INLINE double m_rsqrt(const double& r) const {
    const double         rHalf = 0.5 * r;
    const long long int  tmp   = 0x5FE6EB50C7B537AAl - ( *(long long int*)&r >> 1);//initial guess
    double         rRes        = *(double*)&tmp;
    rRes *= ( 1.5 - ( rHalf * rRes * rRes ) );
    rRes *= ( 1.5 - ( rHalf * rRes * rRes ) );
    return rRes;
  }
  
  //ret norm
  INLINE T invnorm2(void) const {
    const T dr2 = x*x + y*y + z*z;
    return m_rsqrt(dr2);
  }
  INLINE T norm2(void) const {
    const T dr2 = x*x + y*y + z*z;
    return sqrt(dr2);
  }
  INLINE T fastnorm2(void) const {
    const T dr2 = x*x + y*y + z*z;
    return m_rsqrt(dr2) * dr2;
  }
 
  //util
  void elemprint(void) const {
    std::cout << "(" << x << "," << y << "," << z << ")" << std::endl;
  }
};
 
typedef vector3<double> double3;
typedef vector3<float > float3;
 
#undef INLINE
 
double grid_leng = 4.0;
int grid_numb[3];
int all_grid_numb;
double L[3];
 
int GenHash(const double3& r){
  int i = (int)(r.x / grid_leng);
  int j = (int)(r.y / grid_leng);
  int k = (int)(r.z / grid_leng);
 
  if(i == grid_numb[0]) i--;
  if(j == grid_numb[1]) j--;
  if(k == grid_numb[2]) k--;
  
  const int hash = i + grid_numb[0] * (j + k * grid_numb[1]);
  assert(hash < all_grid_numb);
  assert(hash >= 0);
  return hash;
}
 
void MinImage(double3& a){
  a.x -= L[0] * (int) (2. * a.x / L[0]);
  a.y -= L[1] * (int) (2. * a.y / L[1]);
  a.z -= L[2] * (int) (2. * a.z / L[2]);
}
 
int main(int argc,char* argv[]){
  assert(argc == 3);
  using namespace std;
  const char* curdir = argv[1];
  const int btime = atoi(argv[2]);
  stringstream ss;
  ss << curdir << "/particle_data.txt";
  ifstream fin1(ss.str().c_str());
  ss.str("");
  ss << curdir << "/macro_data.txt";
  ifstream fin2(ss.str().c_str());
  ss.str("");
  
  int wN,hbN;
  int all_time,time_step;
  double buf;
  fin2 >> wN >> hbN >> L[0] >> buf >> all_time >> time_step;
  L[1] = L[0];
  L[2] = L[0];
  const int syssize = wN + hbN;
 
  grid_numb[0] = static_cast<int>(L[0] / grid_leng);
  grid_numb[1] = static_cast<int>(L[1] / grid_leng);
  grid_numb[2] = static_cast<int>(L[2] / grid_leng);
  all_grid_numb = grid_numb[0] * grid_numb[1] * grid_numb[2];

  vector<int    > elem_in_grid(all_grid_numb,0);
  vector<double3> pos_in_grid(all_grid_numb,0.0);
 
  grid_leng = L[0] / grid_numb[0];
  
  string str = curdir;
  str += "/concentrate.txt";
  ofstream fout(str.c_str());
  fout << "#time in_numb out_numb in_radius out_radius in_concentrate out_concentrate" << endl;
 
  vector<double3> pos;
  vector<int>     prp;
  vector<bool>    chm;
  
  for(int time=0; time < all_time; time += time_step){
    //load particle data
    for(int i=0; i<syssize; i++){
      double3 buf; int bprp; bool bchm;
      fin1 >> buf.x >> buf.y >> buf.z >> bprp >> bchm;
      if(bprp != 0){
	pos.push_back(buf); prp.push_back(bprp); chm.push_back(bchm);
      }
    }
    
    if(time > btime){
      //get base pos
      for(size_t i=0; i<pos.size(); i++){
	if(prp[i] == 2){
	  const double3 temp = pos[i];
	  const int temp_hash = GenHash(temp);
	  elem_in_grid[temp_hash]++;
	  pos_in_grid[temp_hash] += temp;
	}
      }
      int max_elem=0;
      size_t max_idx=0;
      for(int i=0; i<all_grid_numb; i++){
	if(elem_in_grid[i] > max_elem){
	  max_elem = elem_in_grid[i];
	  max_idx = i;
	}
      }
      double3 base_pos = pos_in_grid[max_idx];
      base_pos /= elem_in_grid[max_idx];
      
      //move origin = base_pos
      for(size_t i=0; i<pos.size(); i++){
	pos[i] -= base_pos;
	MinImage(pos[i]);
      }
      
      //calc center of mass
      double3 eps = 1e-6;
      double3 dif = 1.0;
      double3 bef_cm_pos(100.0);
      do{
	int cm_elem=0;
	double3 cm_pos(0.0);
	for(size_t i=0; i<pos.size(); i++){
	  if(prp[i] == 2){
	    cm_pos += pos[i];
	    cm_elem++;
	  }
	}
	cm_pos /= cm_elem;

	//move origin => cm_pos
	for(size_t i=0; i<pos.size(); i++){
	  pos[i] -= cm_pos;
	  MinImage(pos[i]);
	}
	dif.x = abs(bef_cm_pos.x - cm_pos.x);
	dif.y = abs(bef_cm_pos.y - cm_pos.y);	
	dif.z = abs(bef_cm_pos.z - cm_pos.z);
	
	bef_cm_pos = cm_pos;
	assert(std::isfinite(cm_pos.x));
	assert(std::isfinite(cm_pos.y));
	assert(std::isfinite(cm_pos.z));

      }while(dif.x > eps.x && dif.y > eps.y && dif.z > eps.z);

      //calc min/max radius of oil particle
      //ASSUME: vesicle radius < 0.5 * L[0]
      double max_rad=0.0; double min_rad=100.0;
      for(size_t i=0; i<pos.size(); i++){
	if(prp[i] == 2){
	  const double cm2r = sqrt(pos[i].x * pos[i].x + pos[i].y * pos[i].y + pos[i].z * pos[i].z);
	  assert(std::isfinite(cm2r));
	  if(cm2r < 0.5 * L[0]){
	    if(cm2r > max_rad) max_rad = cm2r;
	    if(cm2r < min_rad) min_rad = cm2r;
	  }
	}
      }
      
      //calc out/in hydrophilic particle number
      int in_numb=0; int out_numb=0;
      for(size_t i=0; i<pos.size(); i++){
	if(prp[i] == 1 && chm[i] == 0){
	  const double cm2r = sqrt(pos[i].x * pos[i].x + pos[i].y * pos[i].y + pos[i].z * pos[i].z);
	  if(cm2r > max_rad) out_numb++;
	  if(cm2r < min_rad) in_numb++;
	}
      }

      //write data to file
      const double out_conc = out_numb / (L[0]*L[1]*L[2] - 4.*M_PI*max_rad*max_rad*max_rad/3.0);
      if(out_conc > 0.0){
	fout << time  << " " <<  in_numb << " " << out_numb << " " << min_rad << " " << max_rad << " " <<  in_numb / (4.*M_PI*min_rad*min_rad*min_rad/3.0) << " " << out_conc << endl;
      }
    }
    //clear vector
    pos.clear();
    prp.clear();
    chm.clear();
    elem_in_grid.clear();
    pos_in_grid.clear();
  }
}