~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 "CircMNTable2D.h"

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

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

using std::endl;

using std::floor;

CircMNTable2D::CircMNTable2D()
{}

/*!
  Construct CircMNTable2D. 

  \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
*/
CircMNTable2D::CircMNTable2D(const Vector3& MinPt,const Vector3& MaxPt,double cd,unsigned int ngroups):
  MNTable2D(MinPt, MaxPt, cd, ngroups)
{
  set_x_circ();
  // check if grid spacing fits size in circular direction:
  double nx=(MaxPt-MinPt).X()/m_celldim;
  // error message if not
  if(nx!=floor(nx)){
    std::cout << "WARNING! grid spacing " << m_celldim << " doesn't fit periodic x-dimension " << (MaxPt-MinPt).X() << std::endl;
  }
  double shift_x=(m_nx-2)*m_celldim;
  m_shift_vec_x=Vector3(shift_x,0.0,0.0);
}

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

/*!
  set circularity of x-dimension to 1
*/
void CircMNTable2D::set_x_circ()
{
  m_x_periodic=1;
} 

int CircMNTable2D::getXIndex(const Vector3& Pos) const
{
  return int(floor((Pos.x()-m_x0)/m_celldim));
}

int CircMNTable2D::getYIndex(const Vector3& Pos) const
{
  return int(floor((Pos.y()-m_y0)/m_celldim));
}


int CircMNTable2D::getFullIndex(const Vector3& Pos) const
{
  int ix=int(floor((Pos.x()-m_x0)/m_celldim));
  int iy=int(floor((Pos.y()-m_y0)/m_celldim));

  return idx(ix,iy);
}

/*!
  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 CircMNTable2D::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 CircMNTable2D::insert(const Sphere& S,unsigned int gid)
{
  bool res;

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

  if((id!=-1) && (idx!=0) && (idx!=m_nx-1) && (gid<m_ngroups)){ // valid index
    // insert sphere
    m_data[id].insert(S,gid);
    res=true;
    int xidx=getXIndex(S.Center());
    // insert 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);
    } 
  } 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 CircMNTable2D::insertChecked(const Sphere& S,unsigned int gid,double tol)
{
  bool res;
  int id=getIndex(S.Center());
  int idx=getXIndex(S.Center());
  
//   std::cout << "CircMNTable2D::insertChecked(" << S << ") id=" << id << " idx= " << idx << std::endl;

  tol+=s_small_value;
  if((id!=-1) && (idx!=0) && (idx!=m_nx-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;
    int xidx=getXIndex(S.Center());
    // insert clone
    if (xidx==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  (xidx==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);
      }
    } 
  } 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 CircMNTable2D::checkInsertable(const Sphere& S,unsigned int gid)
{
  bool res;

  int id=getIndex(S.Center());
  int idx=getXIndex(S.Center());
  
 if((id!=-1) && (idx!=0) && (idx!=m_nx-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;
}


/*!
  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 CircMNTable2D::generateBonds(int gid,double tol,int btag)
{
  std::cout << "CircMNTable2D::generateBonds( " << gid << " , " << tol << " , " << btag << " )" << std::endl;
  // loop over all cells 
  for(int i=0;i<m_nx-1;i++){
    for(int j=1;j<m_ny-1;j++){
      int id=idx(i,j);
      // loop over "upper" neighbors of each cell
      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((id2==id) && (i!=0)){ // intra-cell, not for boundary 
	    //std::cout << id << " - " << id << std::endl;
	    bonds=m_data[id].getBonds(gid,tol);
	  } else if(id2>id){ // inter-cell
	    //std::cout << id << " - " << id2 << std::endl;
	    bonds=m_data[id].getBonds(gid,tol,m_data[id2]);
	  }
	  for(vector<pair<int,int> >::iterator iter=bonds.begin();
	      iter!=bonds.end();
	      iter++){
	    //std::cout << iter->first << " | " << iter->second << "   "; 
	    m_bonds[btag].insert(*iter);
	  }
	  //std::cout << std::endl;
	}
      }
      //std::cout << std::endl;
    }
  }
}

/*!
  Generate bonds between particles of a group with identical particle tag. Takes cloned particles into account.

  \param gid the group ID
  \param tol max. difference between bond length and equilibrium dist.
  \param btag bond tag
  \param ptag the particle tag
  \param mask the mask for the particle tag
*/
void CircMNTable2D::generateBondsWithMask(int gid,double tol,int btag, int ptag, int mask)
{
  std::cout << "CircMNTable2D::generateBondsWithMask( " << gid << " , " << tol << " , " << btag << " " 
	    << ptag <<" " << mask<<" )" << std::endl;
  // loop over all cells 
  for(int i=0;i<m_nx-1;i++){
    for(int j=1;j<m_ny-1;j++){
      int id=idx(i,j);
      // loop over "upper" neighbors of each cell
      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((id2==id) && (i!=0)){ // intra-cell, not for boundary 
	    //std::cout << id << " - " << id << std::endl;
	    bonds=m_data[id].getBonds(gid,tol,ptag,mask);
	  } else if(id2>id){ // inter-cell
	    //std::cout << id << " - " << id2 << std::endl;
	    bonds=m_data[id].getBonds(gid,tol,m_data[id2],ptag,mask);
	  }
	  for(vector<pair<int,int> >::iterator iter=bonds.begin();
	      iter!=bonds.end();
	      iter++){
	    //std::cout << iter->first << " | " << iter->second << "   "; 
	    m_bonds[btag].insert(*iter);
	  }
	  //std::cout << std::endl;
	}
      }
      //std::cout << std::endl;
    }
  }
}

/*!
  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 CircMNTable2D& T)
{
  if(CircMNTable2D::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=1;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 << " 0 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;
}