material.cpp

用monte carlo方法计算化学反应的代码Oxygen Ion Ordering

#include "material.h"

Material::Material():
  pos_vec(3,0),//length = 3 filled with 0's
  delta_E(3,0),//delta E for each node
  E(0),J1(2.56)
{;}
//------------------------------------------------------------------------------//
Material::~Material()
{
  ;
}

Material& Material::operator=(const Material &m)//copy constructor
{
  if (this != &m)//check if m=m
    {
      pos_vec = m.pos_vec;
      delta_E = m.delta_E;
      for (int i=0;i<3;i++)
	{
	  p_pos_vec[i] = m.p_pos_vec[i];
	  p_delta_E[i] = m.p_delta_E[i];
	}
      //      cout <<"in copy constructor"<<endl;
    }
  return *this;
}
//-------------------------------------------------------------------------------//
void Material::Init(vector< vector< vector <Material> > >  &lat,
		    double &num_O)
{
  //the position vectors were allocated in the constructor
  //now we point the pointers to the vectors
  int *temp;
  double *temp1;
  int i=0, j=0, k=0;
  endpoint = 888;
  E_endpoint = .333333;

  for( i =0; i< lat.size(); i++)
    for( j =0; j< lat[0].size(); j++)
      for( k =0; k<lat[0][0].size(); k++)
	{
	  lat[i][j][k].set_xyz_val(i,j,k);//this sets our local xzy position so we 
	  //don't need pointers to our neighboors
	  //this will be usefull in cal_E
	  if(k != lat[0][0].size()-1)
	    {
	      temp = lat[i][j][k+1].get_address(1);
	      lat[i][j][k].set_address(temp,0);//set the top pointer
	      temp1 = lat[i][j][k+1].get_address_E(1);//for E
	      lat[i][j][k].set_address(temp1,0);//for E
	    }
	  else
	    {
	      lat[i][j][k].set_address(&endpoint,0);
	      lat[i][j][k].set_address(&E_endpoint,0);//for E
	    }
  	  if(i != lat.size()-1)
	    {
	      temp = lat[i+1][j][k].get_address(2);//set right
	      lat[i][j][k].set_address(temp,1);
	      temp1 = lat[i+1][j][k].get_address_E(2);//set right for E
	      lat[i][j][k].set_address(temp1,1);//for E
	    }
	  else
	    {
	      lat[i][j][k].set_address(&endpoint,1);
	      lat[i][j][k].set_address(&E_endpoint,1);//for E
	    }
      	  if(j != lat[0].size()-1)
	    {
	      temp = lat[i][j+1][k].get_address(4);//set back
	      lat[i][j][k].set_address(temp,2);
	      temp1 = lat[i][j+1][k].get_address_E(4);//set back for E
	      lat[i][j][k].set_address(temp1,2);//for E
	    }
  	  else
	    {
	      lat[i][j][k].set_address(&endpoint,2);
	      lat[i][j][k].set_address(&E_endpoint,2);//for E
	    }
  	}

    //now we can setup the oxygen structure
    //note: we could do this above to increase performance when we init the pointers but
  //this is more readable
  vector <int> fill(6,1);//fill it with oxygens
  //so we need a bunch of if statements
  if (num_O == 2.75)
    {
      //something like this
      fill[0] = -1;//so we are missing one oxygen at the top of the fcc struct
    }
  else if (num_O == 2.5)
    {
      vector <int> pat1(6,1), pat2(6,1);
      pat1[left] = -1; pat1[back] = -1;
      pat2[front] = -1; pat2[right] = -1;

      for(i =1; i< lat.size()-1; i+=2)//should start at 1 not 0
	for(j =1; j< lat[0].size()-1; j++)
	  for(k =1; k< lat[0][0].size()-1; k++)
	    {
	      if(j%2) //if j is odd
		{
		  lat[i][j][k].set_nodes(pat1);
		  lat[i+1][j][k].set_nodes(pat2);
		}
	      else//j is even
		{
		  lat[i+1][j][k].set_nodes(pat1);
		  lat[i][j][k].set_nodes(pat2); 
		}
	    }
    }
  else if (num_O == 2.875)
    {
      fill[0] = -1;
      fill[0] = -1;
    }
  else if(num_O == 2)
    {
      
      
      
    }
  else //num_O ==3
    {
      //do nothing, fill should be full
    
      //etc for all the type of structures that we want to do
      //now we loop through the lattice and pass in the oxygen struct
      for(i =1; i< lat.size()-1; i++)//should start at 1 not 0
	for(j =1; j< lat[0].size()-1; j++)
	  for(k =1; k< lat[0][0].size()-1; k++)
	    lat[i][j][k].set_nodes(fill);
    }
}
//----------------------------------------------------------------------------------------//
void Material::set_nodes(const vector<int> &n)
{
    {//top, bottom left, r, f, back
      *p_pos_vec[0]=n[0]; pos_vec[0]=n[1]; pos_vec[1]=n[2];
      *p_pos_vec[1]=n[3]; pos_vec[2]=n[4]; *p_pos_vec[2]=n[5];
    }
}
//--------------------------------------------------------------------------------------//
void Material::set_node(const int &n, const int &node_num)
{
  switch(node_num)
    {
    case 0: 
      *p_pos_vec[0]=n;
      return;
    case 1:
      pos_vec[0]=n;
      return;
    case 2:
      pos_vec[1]=n;
      return;
    case 3:
      *p_pos_vec[1] = n;
      return;
    case 4:
      pos_vec[2] = n;
      return;
    case 5:
      *p_pos_vec[2] = n;
      return;
    default:
      cout<<"bad node in set_node"<<endl;
      return;
    }
}
//---------------------------------------------------------------------------------------------//
void Material::get_nodes(vector<int> &n)//you can get all nodes or just one
{
  //top, bottom left, r, f, back
  n[0] = *p_pos_vec[0]; n[1] = pos_vec[0]; n[2] = pos_vec[1];
  n[3] = *p_pos_vec[1]; n[4] = pos_vec[2]; n[5] = *p_pos_vec[2];
}
//------------------------------------------------------------------------------//
void Material::get_delta_E(vector<double> &n)//you can get all nodes or just one
{
  //top, bottom left, r, f, back
  n[0] = *p_delta_E[0]; n[1] = delta_E[0]; n[2] = delta_E[1];
  n[3] = *p_delta_E[1]; n[4] = delta_E[2]; n[5] = *p_delta_E[2];
}
//---------------------------------------------------------------------------------//
double Material::get_delta_E(const int &node_num)
{
 switch(node_num)
      {
      case 0: 
	return *p_delta_E[0];
      case 1:
	return delta_E[0];
      case 2:
	return delta_E[1];
      case 3:
	return *p_delta_E[1];
      case 4:
	return delta_E[2];
      case 5:
	return *p_delta_E[2];
      default:
	cout<<"bad node in get delta E"<<endl;
	return 0;
      }
  return 0;
}
//------------------------------------------------------------------------------//
int Material::get_node(int node_num)//you can get all nodes or just one
{
  switch(node_num)
      {
      case 0: 
	return *p_pos_vec[0];
      case 1:
	return pos_vec[0];
      case 2:
	return pos_vec[1];
      case 3:
	return *p_pos_vec[1];
      case 4:
	return pos_vec[2];
      case 5:
	return *p_pos_vec[2];
      default:
	cout<<"bad node in get_node"<<endl;
	return 0;
      }
  return 0;
}

//-----------------------------------------------------------------------------//
void Material::print_nodes()
{
  vector<int> temp(6);
  this->get_nodes(temp);
  cout <<temp<<endl;
}
//-----------------------------------------------------------------------------//
void Material::set_address(int *p,int pos)
{
    p_pos_vec[pos] = p;
}
void Material::set_address(double *p, int pos)
{
  p_delta_E[pos] = p;
}
//-----------------------------------------------------------------------------//	    
int * Material::get_address(const int &position)
{
  int *temp;
  switch(position)
    {
    case 1://bottom
      temp = &pos_vec[0];
      break;
    case 2://left
      temp = &pos_vec[1];
      break;
    case 4://front
      temp = &pos_vec[2];
      break;
    default:
      cout <<"WRONG POSITION FROM GET_ADDRESS!!!"<<endl;
      break;
    }
  return temp;
}
//-------------------------------------------------------------------------------//
double * Material::get_address_E(const int &position)
{
  double *temp;
  switch(position)
    {
    case 1://bottom
      temp = &delta_E[0];
      break;
    case 2://left
      temp = &delta_E[1];
      break;
    case 4://front
      temp = &delta_E[2];
      break;
    default:
      cout <<"WRONG POSITION FROM GET_ADDRESS!!!"<<endl;
      break;
    }
 return temp;
}

//-------------------------------------------------------------------------------//
void Material::set_xyz_val(const int &i,const int &j,const int &k)
{
  x = i;
  y = j;
  z = k;
}
//-------------------------------------------------------------------------------//
double Material::calc_E(const int &node, vector< vector< vector< Material> > > &lat)
{
  //we are only taking the first neighboors right now
  //we always have 8 neighbors unless we are an endpoint
  double sum=0;
  double nu = 20;
  vector<int> neighboors(8,0);
  this->get_neighboors_of_node(node,neighboors,lat); //get our neighboors
  int my_value = this->get_node(node);  //get my value

  for (int i=0; i<neighboors.size();i++)
    sum += my_value * neighboors[i];
  sum *= J1;
  sum = sum - (my_value*nu);

  //  cout <<"mynode = "<<node<<" sum = "<<sum<<" neighboors = " <<neighboors<<" size = " <<lat.size()<< " "<<lat[0].size()<<" " <<lat[0][0].size()<<endl;
  //  sum =  (-2)*sum;
  this->set_delta_E(node,(2*(nu*my_value - sum)));//delta_E = finalE - initialE = -2*initialE
  return sum;
}
//---------------------------------------------------------------//
double Material::get_super_node_E()
{
  return E;
}
//---------------------------------------------------------------//
void Material::set_delta_E(const int &node, const double &value)
{
  switch(node)
    {
    case 0: 
      *p_delta_E[0]=value;
      return;
    case 1:
      delta_E[0]=value;
      return;
    case 2:
      delta_E[1]=value;
      return;
    case 3:
      *p_delta_E[1] = value;
      return;
    case 4:
      delta_E[2] = value;
      return;
    case 5:
      *p_delta_E[2] = value;
      return;
    default:
      cout<<"bad node in set_node"<<endl;
      return;
    }
}
//------------------------------------------------------------------------------//
void Material::get_neighboors_of_node(const int &node, vector<int> &v,
				      vector< vector< vector< Material> > > &lat)
{
  vector<int> temp;

  if(node == 0)//we want the neighbors of the top node
    {
      v[0] = lat[x][y][z+1].get_node(2);//get the topleft node
      v[1] = lat[x][y][z+1].get_node(3);//get the topright node
      v[2] = lat[x][y][z].get_node(2);//get the bottomleft node
      v[3] = lat[x][y][z].get_node(3);//get the bottomright node

      v[4] = lat[x][y][z+1].get_node(4);//get the fronttop node
      v[5] = lat[x][y][z].get_node(4);//get the frontbottom node
      v[6] = lat[x][y][z+1].get_node(5);//get the backtop node
      v[7] = lat[x][y][z].get_node(5);//get the backbottom node
    }
  if(node == 1)//we want the neighboors of the bottom node
    {
      v[0] = lat[x][y][z].get_node(2);//get the topleft node
      v[1] = lat[x][y][z].get_node(3);//get the topright node
      v[2] = lat[x][y][z-1].get_node(2);//get the bottomleft node
      v[3] = lat[x][y][z-1].get_node(3);//get the bottomright node

      v[4] = lat[x][y][z].get_node(4);//get the fronttop node
      v[5] = lat[x][y][z-1].get_node(4);//get the frontbottom node
      v[6] = lat[x][y][z].get_node(5);//get the backtop node
      v[7] = lat[x][y][z-1].get_node(5);//get the backbottom node
    }
  if(node == 2)//we want the neighboors of the left node
    {
      v[0] = lat[x-1][y][z].get_node(0);//get the topleft node
      v[1] = lat[x][y][z].get_node(0);//get the topright node
      v[2] = lat[x-1][y][z].get_node(1);//get the bottomleft node
      v[3] = lat[x][y][z].get_node(1);//get the bottomright node

      v[4] = lat[x-1][y][z].get_node(4);//get the leftfront node
      v[5] = lat[x-1][y][z].get_node(5);//get the leftback node
      v[6] = lat[x][y][z].get_node(4);//get the rightfront node
      v[7] = lat[x][y][z].get_node(5);//get the rightback node
    }
  if(node == 3)//we want the neighboors of the right node
    {
      v[0] = lat[x][y][z].get_node(0);//get the topleft node
      v[1] = lat[x+1][y][z].get_node(0);//get the topright node
      v[2] = lat[x][y][z].get_node(1);//get the bottomleft node
      v[3] = lat[x+1][y][z].get_node(1);//get the bottomright node

      v[4] = lat[x][y][z].get_node(4);//get the leftfront node
      v[5] = lat[x][y][z].get_node(5);//get the leftback node
      v[6] = lat[x+1][y][z].get_node(4);//get the rightfront node
      v[7] = lat[x+1][y][z].get_node(5);//get the rightback node
    }
  if(node == 4)//we want to get the neighboors of the front node
    {
      v[0] = lat[x][y-1][z].get_node(0);//get the fronttop node
      v[1] = lat[x][y-1][z].get_node(1);//get the frontbottom node
      v[2] = lat[x][y][z].get_node(0);//get the backtop node
      v[3] = lat[x][y][z].get_node(1);//get the backbottom node
      
      v[4] = lat[x][y-1][z].get_node(2);//get the frontleft node
      v[5] = lat[x][y-1][z].get_node(3);//get the frontright node
      v[6] = lat[x][y][z].get_node(2);//get the backleft node
      v[7] = lat[x][y][z].get_node(3);//get the backright node
    }
  if(node ==5)//we want to get the neighboors of the back node
    {
      v[0] = lat[x][y][z].get_node(0);//get the fronttop node
      v[1] = lat[x][y][z].get_node(1);//get the frontbottom node
      v[2] = lat[x][y+1][z].get_node(0);//get the backtop node
      v[3] = lat[x][y+1][z].get_node(1);//get the backbottom node
      
      v[4] = lat[x][y][z].get_node(2);//get the frontleft node
      v[5] = lat[x][y][z].get_node(3);//get the frontright node
      v[6] = lat[x][y+1][z].get_node(2);//get the backleft node
      v[7] = lat[x][y+1][z].get_node(3);//get the backright node
    }
}