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

// --- System includes ---
#include <cmath>
#include <cstdlib>
#include <sys/time.h>

using std::ceil;
using std::sqrt;

// --- project includes ---
#include "sphere_fitting/Sphere2DFitter.h"

InsertGenerator2D::InsertGenerator2D()
{}

/*!
  Destructor
*/
InsertGenerator2D::~InsertGenerator2D()
{}

/*!
  Constructor

  \param rmin minimum particle radius
  \param rmax maximum particle radius
  \param ntries max. nr. of tries to insert particle
  \param max_iter maximum iterations within the iterative solvers
  \param prec max. error in iterative solvers
*/
InsertGenerator2D::InsertGenerator2D(double rmin,double rmax,int tries, int max_iter,double prec)
{
  m_rmin=rmin;
  m_rmax=rmax;
  m_max_tries=tries;
  m_max_iter=max_iter;
  m_prec=prec;
  m_old_seeding=false;
}

/*!
  Constructor. Added random seed parameter

  \param rmin minimum particle radius
  \param rmax maximum particle radius
  \param ntries max. nr. of tries to insert particle
  \param max_iter maximum iterations within the iterative solvers
  \param prec max. error in iterative solvers
  \param seed if true, initialize random number generator via time
*/
InsertGenerator2D::InsertGenerator2D(double rmin,double rmax,int tries, int max_iter,double prec, bool seed)
{
  if(seed){
    struct timeval tv;
    gettimeofday(&tv,NULL);    
    int random_seed=tv.tv_usec;
    srand(random_seed);
  }
  m_rmin=rmin;
  m_rmax=rmax;
  m_max_tries=tries;
  m_max_iter=max_iter;
  m_prec=prec;
  m_old_seeding=false;
}

/*!
  toggle between old (pre rev.84, ansiotropic) and new (rev 84 onwards) seeding behaviour

  \param old if true -> old behaviour, if false -> new behaviour
*/
void InsertGenerator2D::setOldSeeding(bool old)
{
  m_old_seeding=old;
}

/*!
  generate packing
  
  \param vol a pointer to the volume in which the packing is generated
  \param ntable a pointer to the neighbour table used 
  \param gid particle group id
  \param tag the particle tag
*/
void InsertGenerator2D::generatePacking (AVolume2D* vol,MNTable2D* ntable,int gid, int tag)
{
  // use this-> to accomodate virtual seed/fill functions
  this->seedParticles(vol,ntable,gid,tag);
  this->fillIn(vol,ntable,gid,tag);
}

/*!
  seed the area with particles

  \param vol a pointer to the packing volume
  \param ntable a pointer to the neighbour table used 
  \param gid particle group id
  \param tag the particle tag
*/
void InsertGenerator2D::seedParticles(AVolume2D* vol,MNTable2D* ntable,int gid, int tag)
{
  // get bounding box
  pair<Vector3,Vector3> bbx=vol->getBoundingBox();
  double dx=(bbx.second.X()-bbx.first.X())-2.0*m_rmin;
  double dy=(bbx.second.Y()-bbx.first.Y())-2.0*m_rmin;
  // get index limits for seeding
  int imax=int(ceil(dx/(m_rmax*2.0)));
  int jmax=int(ceil(dy/(m_rmax*sqrt(3.0))));
  // seed positions
  for(int i=0;i<imax;i++){
    for(int j=0;j<jmax;j++){
      // get position
      double px=bbx.first.X()+m_rmin+(double(i)+0.5*double(j%2))*m_rmax*2.0;
      double py=bbx.first.Y()+m_rmin+double(j)*sqrt(3.0)*m_rmax;

      // get dist to egde
      double dex=(bbx.second.X()-px) < (px-bbx.first.X()) ? bbx.second.X()-px  : px-bbx.first.X();
      double dey=(bbx.second.Y()-py) < (py-bbx.first.Y()) ? bbx.second.Y()-py  : py-bbx.first.Y();
      
      double de=(dex<dey) ? dex : dey;
      
      // check max rad.
      if(de>m_rmin){
	// calc random radius
	double r;
	double jitter;
	if(de<m_rmax) {
	  if(m_old_seeding){
	    r=m_rmin+((de-m_rmin)*((double)(rand())/(double)(RAND_MAX)));
	    jitter=0.0;
	  } else {
	    r=m_rmin+(((de-m_rmin)/2.0)*((double)(rand())/(double)(RAND_MAX)));
	    jitter=de-r;
	  }
	} else {
	  if(m_old_seeding){
	    r=m_rmin+((m_rmax-m_rmin)*((double)(rand())/(double)(RAND_MAX)));
	    jitter=0.0;
	  } else {
	    r=m_rmin+(((m_rmax-m_rmin)/2.0)*((double)(rand())/(double)(RAND_MAX)));
	    jitter=m_rmax-r;
	  }
	}
	r-=m_prec;
	// jitter position 
	double dx=jitter*(2.0*((double)(rand())/(double)(RAND_MAX))-1.0);
	double dy=jitter*(2.0*((double)(rand())/(double)(RAND_MAX))-1.0);
	Sphere S(Vector3(px+dx,py+dy,0.0),r);
	bool fit_in=vol->isIn(S); 
	bool fit_check=ntable->checkInsertable(S,gid);
	if(fit_in && fit_check){
	  S.setTag(tag);
	  ntable->insertChecked(S,gid);	
	} 
      }
    }
  }

}

/*!
  Fill in particles

  \param vol a pointer to the packing volume
  \param ntable a pointer to the neighbour table used 
  \param gid particle group id
  \param tag the particle tag
*/
void InsertGenerator2D::fillIn(AVolume2D* vol,MNTable2D* ntable,int gid, int tag)
{
  Sphere nsph;

  int total_tries=0;
  int count_insert=0;

  int nvol=vol->getNumberSubVolumes();
  for(int ivol=0;ivol<nvol;ivol++){
    int countfail=0; // number of failed attenpts since last successfull insertion
    while(countfail<m_max_tries){
      bool findfit=false;
      Vector3 P=vol->getAPoint(ivol); // get random point within volume
      multimap<double,const Sphere*> close_particles=ntable->getSpheresClosestTo(P,3); // get 3 nearest spheres
      map<double,const Line2D*> close_lines=vol->getClosestPlanes(P,2); // get 2 nearest planes
      map<double,const AGeometricObject*> close_objects=vol->getClosestObjects(P,2); // get 2 nearest "objects"


      map<double,const AGeometricObject*> geomap;
      geomap.insert(close_particles.begin(),close_particles.end());
      geomap.insert(close_lines.begin(),close_lines.end());
      geomap.insert(close_objects.begin(),close_objects.end());
      
      if(geomap.size()>=3){
	map<double,const AGeometricObject*>::iterator iter=geomap.begin();
	const AGeometricObject* GO1=iter->second;iter++;
	const AGeometricObject* GO2=iter->second;iter++;
	const AGeometricObject* GO3=iter->second;
	nsph=FitSphere2D(GO1,GO2,GO3,P,m_max_iter,m_prec);
	findfit=true;
      }

      if(findfit){
	// check if within radius range
	bool is_radius=(m_rmin<nsph.Radius()) && (m_rmax>nsph.Radius());
	
	// check inside volume, intersections ...
	bool fit=vol->isIn(nsph) && ntable->checkInsertable(nsph,gid);
	
	if(fit && is_radius){ // acceptable particle 
	  nsph.setTag(tag);
	  ntable->insertChecked(nsph,gid); // insert
	  count_insert++;
	  if((count_insert%100)==0) std::cout << "inserted: " << count_insert << std::endl;
	  total_tries+=countfail;
	  if(countfail>m_max_tries/10) std::cout << countfail << " tries" << std::endl;
	  countfail=0; // reset failure counter
	} else countfail++; 
      } else countfail++;  
    }
  }
  std::cout << "total tries: " << total_tries << std::endl;
}

ostream& operator << (ostream& ost,const InsertGenerator2D& T) 
{
   return ost;
}