~esys-p-dev/esys-particle/gengeo

/////////////////////////////////////////////////////////////
//                                                         //
// Copyright (c) 2007-2014 by The University of Queensland //
// Centre for Geoscience Computing                         //
// http://earth.uq.edu.au/centre-geoscience-computing      //
//                                                         //
// Primary Business: Brisbane, Queensland, Australia       //
// Licensed under the Open Software License version 3.0    //
// http://www.opensource.org/licenses/osl-3.0.php          //
//                                                         //
/////////////////////////////////////////////////////////////

#include "CircMNTableXY2D.h"

// --- System includes ---
#include <cmath>

// --- IO includes ---
#include <iostream>

using std::endl;

using std::floor;

CircMNTableXY2D::CircMNTableXY2D()
{}

/*!
  Construct CircMNTableXY2D. 

  \param MinPt minimum point (z component ignored)
  \param MaxPt maximum point (z component ignored)
  \param cd cell dimension
  \param ngroups initial number of particle groups
*/
CircMNTableXY2D::CircMNTableXY2D(const Vector3& MinPt,const Vector3& MaxPt,double cd,unsigned int ngroups):
  CircMNTable2D(MinPt, MaxPt, cd, ngroups)
{
  set_y_circ();

  // check if grid spacing fits size in circular direction:
  double ny=(MaxPt-MinPt).Y()/m_celldim;
  // error message if not
  if(ny!=floor(ny)){
    std::cerr << "WARNING! grid spacing " << m_celldim << " doesn't fit periodic y-dimension " << (MaxPt-MinPt).Y() << std::endl;
  }
  double shift_y=(m_ny-2)*m_celldim;
  m_shift_vec_y=Vector3(0.0,shift_y,0.0);
}

/*!
  Destruct CircMNTable2D. Just calls MNTable2D destructor.
*/
CircMNTableXY2D::~CircMNTableXY2D()
{}

/*!
  set circularity of y-dimension to 1
*/
void CircMNTableXY2D::set_y_circ()
{
  m_y_periodic=1;
} 

/*!
  get the cell index for a given position

  \param Pos the position
  \return the cell index if Pos is inside the table, -1 otherwise
*/
int CircMNTableXY2D::getIndex(const Vector3& Pos) const
{
  int ret;

  int ix=int(floor((Pos.x()-m_x0)/m_celldim));
  int iy=int(floor((Pos.y()-m_y0)/m_celldim));

  // check if pos is in table excl. padding
  if((ix>=0) && (ix<=m_nx-1) && (iy>=0) && (iy<=m_ny-1)){
    ret=idx(ix,iy);
  } else {
    ret=-1;
  }
  
  return ret;
}

/*!
  Insert sphere. Insert clone into other side of Table.

  \param S the Sphere
  \param gid the group id
*/
bool CircMNTableXY2D::insert(const Sphere& S,unsigned int gid)
{
  bool res;

  int id=getIndex(S.Center());
  int idx=getXIndex(S.Center());
  int idy=getYIndex(S.Center());


  if((id!=-1) && (idx!=0) && (idx!=m_nx-1) && (idy!=0) && (idy!=m_ny-1) && (gid<m_ngroups)){ // valid index
    // insert sphere
    m_data[id].insert(S,gid);
    res=true;
    int xidx=getXIndex(S.Center());
    int yidx=getYIndex(S.Center());
    // insert x-clone
    if (xidx==1){
      Sphere SClone=S;
      SClone.shift(m_shift_vec_x);
      int clone_id=getFullIndex(SClone.Center());
      m_data[clone_id].insert(SClone,gid);
    } else if  (xidx==m_nx-2){
      Sphere SClone=S;
      SClone.shift(-1.0*m_shift_vec_x);
      int clone_id=getFullIndex(SClone.Center());
      m_data[clone_id].insert(SClone,gid);
    } 
    // insert y-clone
    if (yidx==1){
      Sphere SClone=S;
      SClone.shift(m_shift_vec_y);
      int clone_id=getFullIndex(SClone.Center());
      m_data[clone_id].insert(SClone,gid);
    } else if  (yidx==m_ny-2){
      Sphere SClone=S;
      SClone.shift(-1.0*m_shift_vec_y);
      int clone_id=getFullIndex(SClone.Center());
      m_data[clone_id].insert(SClone,gid);
    } 
  } else {
    res=false;
  }

  return res;
}


/*!
  insert sphere if it doesn't collide with other spheres

  \param S the Sphere
  \param gid the group id
*/
bool CircMNTableXY2D::checkInsertable(const Sphere& S,unsigned int gid)
{
 bool res;

  int id=getIndex(S.Center());
  int idx=getXIndex(S.Center());
  int idy=getYIndex(S.Center());

  if((id!=-1) && (idx!=0) && (idx!=m_nx-1) && (idy!=0) && (idy!=m_ny-1) && (gid<m_ngroups)){
    multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(S.Center(),S.Radius()-s_small_value,gid);
    if(close_spheres.size()==0){ 
      res=true;
    } else res=false;
  } else {
    res=false;
  }

  return res;
}


/*!
  Insert sphere if it doesn't collide with other spheres. Insert clone into other side of Table.

  \param S the Sphere
  \param gid the group id
*/
bool CircMNTableXY2D::insertChecked(const Sphere& S,unsigned int gid,double tol)
{
  bool res;
  int id=getIndex(S.Center());
  int idx=getXIndex(S.Center());
  int idy=getYIndex(S.Center());
  
//   std::cerr << "CircMNTable2D::insertChecked(" << S << ") id=" << id << " idx= " << idx << std::endl;

  tol+=s_small_value;
  if((id!=-1) && (idx!=0) && (idx!=m_nx-1) && (idy!=0) && (idy!=m_ny-1) && (gid<m_ngroups)){
    // insert original
    multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(S.Center(),S.Radius()-tol,gid);
    if(close_spheres.size()==0){
      m_data[id].insert(S,gid);
      res=true;
    } else res=false;
    // insert x-clone
    if (idx==1){
      Sphere SClone=S;
      SClone.shift(m_shift_vec_x);
      multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(SClone.Center(),SClone.Radius()-tol,gid);
      if(close_spheres.size()==0){
	int clone_id=getFullIndex(SClone.Center());
	m_data[clone_id].insert(SClone,gid);
      }
    } else if  (idx==m_nx-2){
      Sphere SClone=S;
      SClone.shift(-1.0*m_shift_vec_x);
      multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(SClone.Center(),SClone.Radius()-tol,gid);
      if(close_spheres.size()==0){
	int clone_id=getFullIndex(SClone.Center());
	m_data[clone_id].insert(SClone,gid);
      }
    }
    // insert y-clone
    if (idy==1){
      Sphere SClone=S;
      SClone.shift(m_shift_vec_y);
      multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(SClone.Center(),SClone.Radius()-tol,gid);
      if(close_spheres.size()==0){
	int clone_id=getFullIndex(SClone.Center());
	m_data[clone_id].insert(SClone,gid);
      }
    } else if  (idy==m_ny-2){
      Sphere SClone=S;
      SClone.shift(-1.0*m_shift_vec_y);
      multimap<double,const Sphere*> close_spheres=getSpheresFromGroupNear(SClone.Center(),SClone.Radius()-tol,gid);
      if(close_spheres.size()==0){
	int clone_id=getFullIndex(SClone.Center());
	m_data[clone_id].insert(SClone,gid);
      }
    }
  } else {
    res=false;
  }

  return res;
}

/*!
  Generate bonds between particles of a group. Takes cloned particles into account.

  \param gid the group ID
  \param tol max. difference between bond length and equilibrium dist.
  \param btag bond tag
*/
void CircMNTableXY2D::generateBonds(int gid,double tol,int btag)
{
  std::cerr << "CircMNTableXY2D::generateBonds( " << gid << " , " << tol << " , " << btag << " )" << std::endl;
  // loop over all cells 
  for(int i=0;i<m_nx-1;i++){
    for(int j=0;j<m_ny-1;j++){
      int id=idx(i,j);
      for(int ii=-1;ii<=1;ii++){
	for(int jj=-1;jj<=1;jj++){
	  int id2=idx(i+ii,j+jj);
	  vector<pair<int,int> > bonds;
	  if(((ii+jj)==0)){ // intra-cell, not for boundary 
	    bonds=m_data[id].getBonds(gid,tol);
	  } else if(id2>id){ // inter-cell
   	    bonds=m_data[id].getBonds(gid,tol,m_data[id2]);
	  } 
	  for(vector<pair<int,int> >::iterator iter=bonds.begin();
	      iter!=bonds.end();
	      iter++){	    
	    m_bonds[btag].insert(*iter);
	  }
	}
      }
    }
  }
}

/*!
  Output the content of a MNTable2D to an ostream. The output format depends on the value of MNTable2D::s_output_style (see  MNTable2D::SetOutputStyle). If it is set to 0, output suitable for debugging is generated, if it is set to 1 output in the esys .geo format is generated. If MNTable2D::s_output_style is set to 2, the output format is VTK-XML.

  \param ost the output stream
  \param T the table
*/
ostream& operator << (ostream& ost,const CircMNTableXY2D& T)
{
  if(CircMNTableXY2D::s_output_style==0){ // debug style
    // don't print padding
    MNTCell::SetOutputStyle(0);
    for(int i=0;i<T.m_nx;i++){
      for(int j=0;j<T.m_ny-1;j++){
	ost << "=== Cell " << i << " , " << j << " === " << endl;
	ost << T.m_data[T.idx(i,j)];
      }
    }
  } else if (CircMNTable2D::s_output_style==1){ // geo style
    int nparts=0;
    for(int i=1;i<T.m_nx-1;i++){
      for(int j=1;j<T.m_ny-1;j++){
	nparts+=T.m_data[T.idx(i,j)].NParts();
      }
    }
    // header
    ost << "LSMGeometry 1.2" << endl;
    ost << "BoundingBox " << T.m_x0+T.m_celldim << " " <<  T.m_y0+T.m_celldim << " 0.0 " << T.m_x0+double(T.m_nx-1)*T.m_celldim << " " <<  T.m_y0+double(T.m_ny-1)*T.m_celldim << " 0.0 " << endl;
    ost << "PeriodicBoundaries "  << T.m_x_periodic << " " << T.m_y_periodic << "  0" << endl;
    ost << "Dimension 2D" << endl;
    // particles
    ost << "BeginParticles" << endl;
    ost << "Simple" << endl;
    ost << nparts << endl;
    MNTCell::SetOutputStyle(1);
     for(int i=1;i<T.m_nx-1;i++){
      for(int j=1;j<T.m_ny-1;j++){
	ost <<  T.m_data[T.idx(i,j)];
      }
    }
    ost << "EndParticles" << endl;
    // bonds
    ost << "BeginConnect" << endl;
    // sum bonds
    int nbonds=0;
    for(map<int,set<pair<int,int> > >::const_iterator iter=T.m_bonds.begin();
	iter!=T.m_bonds.end();
	iter++){
      nbonds+=iter->second.size();
    }
    ost << nbonds << endl;
    for(map<int,set<pair<int,int> > >::const_iterator iter=T.m_bonds.begin();
	iter!=T.m_bonds.end();
	iter++){
      for(set<pair<int,int> >::const_iterator v_iter=iter->second.begin();
	  v_iter!=iter->second.end();
	  v_iter++){
	ost << v_iter->first << " " << v_iter->second << " " << iter->first << endl;
      }
    }
    ost << "EndConnect" << endl;
  } else if (CircMNTable2D::s_output_style==2){ // vtk-xml
    T.WriteAsVtkXml(ost);
  }

  return ost;
}