dgParticleTracker2D.h

#ifndef _DG_PARTICLE_TRACKER_2D_H_
#define _DG_PARTICLE_TRACKER_2D_H_

#include <fstream>
#include <iomanip>
#include <sstream>
#include <iostream>
#include <vector>
#include"math.h"
#include"GModel.h"
#include"SPoint3.h"
#include"MElement.h"
#include"dgDofContainer.h"
#include"dgFunctionEvaluator.h"
#include"dgGroupOfElements.h"
#include"dgNeighbourMap.h"

double randoraw();
void rando(double *randoo, int M);
void Inverse(double Matrix[3][3], double Inverse[3][3]);
void Transpose(double Matrix[3][3]);
void projection(double* v1, double* v2, double* v1_proj_on_v2);
double projectionMag(double* v1, double* v2);
void convertVectorCoords(double target_v[3],double original_v[3], double jac[3][3]);
SPoint2 intersection(SPoint2 p1, SPoint2 p2, SPoint2 p3, SPoint2 p4, int& flag);

class particle {
private:
  double _xcoord, _ycoord;
  int _competency;      //0: pre-competent, 1: competent, 2: post-competent
  int _int_state;       //0: alive, 1: settled, 2: dead, 3: out of domain
  int _source_reef, _destination_reef;
  //Remember element and group Id particle is in:
  int _groupId, _elementId;
  
public:
  // Default constructor
  particle()
  {
    _xcoord = 0.0;
    _ycoord = 0.0;
    _int_state = 0;
    _source_reef = -1;
    _competency = 0;
  }
  
  // Cartesian constructor which SETS mesh element Id's
  particle(double x0, double y0, int iGroup, int iElement, int status)
  {
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = -1;
    _groupId = iGroup;
    _elementId = iElement;
    _competency = 0;
  }
  
  // Cartesian constructor which DOESN'T set mesh element Id's
  particle(double x0, double y0, int status)
  {
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = -1;
    _groupId = -1;
    _elementId = -1;
    _competency = 0;
  }
  
  // Cartesian constructor which sets source reef AND mesh element Id's
  particle(double x0, double y0, int iGroup, int iElement, int srcReef, int status)
  {
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = srcReef;
    _groupId = iGroup;
    _elementId = iElement;
    _competency = 0;
  }
  
  // Cartesian constructor which sets source reef but NOT mesh element Id's (try not to use this, doesn't store element we're in)
  particle(double x0, double y0, int srcReef, int status)
  {
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = srcReef;
    _groupId = -1;
    _elementId = -1;
    _competency = 0;
  }
  
  //Access methods:
  double x()
  {
    return _xcoord;
  }
  
  double y()
  {
    return _ycoord;
  }
  
  int getState()
  {
    return _int_state;
  }

  int getSource() 
  {
    return _source_reef;
  }
  
  void getElementId(int *iGroup, int *iElement)
  {
    *iGroup = _groupId;
    *iElement = _elementId;
  }
  
  int getCompetency()
  {
    return _competency;
  }
  
  void setPosition(double x, double y)
  {
    _xcoord=x;
    _ycoord=y;
  }
  
  void setState(int status)
  {
    _int_state = status;
  }
  
  void setElementId(int iGroup, int iElement)
  {
    _groupId = iGroup;
    _elementId = iElement;
  }
  
  void setCompetency(int neoComp)
  {
    _competency = neoComp;
  }
  
};

class particleArray {
private:
  std::vector<particle> _particle_array;
  int _alive, _settled, _dead, _outOfDomain, _aliveAndCompetent;
  
public:
  //Constructors:
  particleArray()
  {
    std::cout<<"Error: Default particleArray constructor used.\n";
  }
  
  /** Single point constructor which sets element IDs*/
  particleArray(dgGroupCollection* _group, int M, double x0, double y0, int status, int source_reef)
  {
    //Find group & element IDs
    SPoint3 xyz(x0,y0,0.0);
    GModel *gm = _group->getModel();
    MElement *e = gm->getMeshElementByCoord(xyz, 2);
    int iGroup, iElem;
    _group->find(e, iGroup, iElem); 
    
    for (int i=0; i<M; i++)
    {
      //Seed particle
      _particle_array.push_back( particle(x0,y0,iGroup,iElem,source_reef,status) );
    }
    _alive = M;
    _settled = 0;
    _dead = 0;
    _outOfDomain = 0;
    _aliveAndCompetent = 0;
  }
  
  /**Seed particles using reef map DC*/
  particleArray(dgGroupCollection* groups, dgDofContainer* reefsDC, fullMatrix<double>* initCM, fullMatrix<double>* reefsToSeed, double concentration, int seedOverWhichReefs, int nbOfShallowReefs);
    
  void addParticles(dgGroupCollection* _group, int M, double x0, double y0, int status, int source_reef)
  {
    //Find group & element IDs
    SPoint3 xyz(x0,y0,0.0);
    GModel *gm = _group->getModel();
    MElement *e = gm->getMeshElementByCoord(xyz, 2);
    int iGroup, iElem;
    _group->find(e, iGroup, iElem);
    
    for (int i=0; i<M; i++)
    {      
      //Seed particle
      _particle_array.push_back( particle(x0,y0,iGroup,iElem,source_reef,status) );
    }
    _alive += M;

  }
  
  void setStats(int *stats) {
    _alive = stats[0];
    _settled = stats[1];
    _dead = stats[2];
    _outOfDomain = stats[3];
    _aliveAndCompetent = stats[4];
  }
  
  void getStats(int *stats) {
    stats[0] = _alive;
    stats[1] = _settled;
    stats[2] = _dead;
    stats[3] = _outOfDomain;
    stats[4] = _aliveAndCompetent;
  }  
  
  int printNbParticles() {
    return static_cast<int>(_particle_array.size());
  }
  
  //Access a given particle:
  particle* getParticle(int m)
  {
    return &_particle_array[m];
  }

  double getParticle_x(int m)
  {
    return _particle_array[m].x();
  }
  
  double getParticle_y(int m)
  {
    return _particle_array[m].y();
  } 
  
  /** Print positions of all particles.
   *  Generates 1 file per call with particles colour-coded by source reef.
   *  Only records particles of state 'status'
   *    filetype = dat or pos 
   */
  void printPositions(const std::string filename, const std::string filetype, int status);
  
  /** Print positions of all particles. 
   *  1 file per particle, colour-coded by time-step.
   *  filetype = dat or pos
   */
  void printPositionsFollow(const std::string filename, const std::string filetype, int n, int N, double dt);
  
  /** Print position of one particle. 
   *  Colour-coded by time-step.
   *  filetype = dat or pos
   */
  void printPositionsFollowSingle(const std::string filename, const std::string filetype, int n, int N, double dt, int particleId);
};

class connectivityMatrixPrint {
private:
  
public:
  /** Constructor */
  connectivityMatrixPrint() {
  }
  /** Method to print contents of matrix */
  void CMPrinter(fullMatrix<double>* FM, const std::string filename);
};

class dgParticleTracker2D {
private:
  double _dt;
  int _N;
  int _settleType, _mortType, _swimType, _settleOverWhichReefType, _nbOfShallowReefs;
  double _t_preComp, _t_postComp, _t_settMax, _settRate_max, _compAcquisitionRate, _compLossRate;
  double _mortRate_max, _mortLambda, _mortNu;
  double _swimSpeed, _swimDist, _swimStart;
  dgGroupCollection *_group;
  particleArray *_ParticleArray;
  dgDofContainer* _bathymetry;
  dgDofContainer* _diffusivity;
  dgDofContainer* _reefs;
  dgDofContainer* _distToReefs;
  dgDofContainer* _directionToReefs;
  GModel* _gm;
  dgNeighbourMap* _neighbourMap;
  std::string _LP_OutputDir;
  bool _mortFlag, _settleFlag, _swimFlag, _connectivityExposureTimeFlag;
  
public:
  /** Constructor */
  dgParticleTracker2D(dgGroupCollection* group, particleArray* ParticleArray, dgDofContainer* bathymetry, dgDofContainer* diffusivity,  int N, double dt, const std::string LP_OutputDir) {
    _dt = dt;
    _N = N;
    _group = group;
    _ParticleArray = ParticleArray;
    _bathymetry = bathymetry;
    _diffusivity = diffusivity;
    _LP_OutputDir = LP_OutputDir;
    
    //By default, set mortality, settling, swimming to OFF
    _mortFlag = false;
    _settleFlag = false;
    _swimFlag = false;
    
    //Load model:
    _gm = _group->getModel();
    _neighbourMap = new dgNeighbourMap(*_group);
    
    //Initialise output files
    std::ostringstream fname;
    std::ofstream outStats;
    //Define filename
    fname << _LP_OutputDir << "/particleStats.dat";
    //Open file (replace if already exists)
    std::string temp = fname.str();
    outStats.open( temp.c_str(), std::ios::out | std::ios::trunc );
    outStats.close();
  }
  
  ~dgParticleTracker2D() {
    delete _neighbourMap;
  }
  
  /** Activate particle mortality and load parameters */
  void addMortality(fullMatrix<double>* mortParams) {
    _mortFlag = true;
    _mortType = static_cast<int>(mortParams->get(0,0)+0.1);
    _mortRate_max = mortParams->get(0,1);
    _mortLambda = mortParams->get(0,2);   //(Connolly & Baird 2010)lambda = 1/b(Pinder et al. 1978)
    _mortNu = mortParams->get(0,3);       //(Connolly & Baird 2010)nu = c(Pinder et al. 1978)
  }
  
  /** Activate particle settling and load parameters */
  void addSettling(fullMatrix<double>* settleParams, dgDofContainer* reefs, int nbOfShallowReefs) {
    _settleFlag = true;
    _settleType = static_cast<int>(settleParams->get(0,0)+0.1);
    _t_preComp = settleParams->get(0,1);
    _t_postComp = settleParams->get(0,2);
    _t_settMax = settleParams->get(0,3);
    _settRate_max = settleParams->get(0,4);
    _compAcquisitionRate = settleParams->get(0,5);    //(Connolly & Baird 2010): a
    _compLossRate = settleParams->get(0,6);           //(Connolly & Baird 2010): b
    _settleOverWhichReefType = settleParams->get(0,7);
    //Convert competency acquisition & loss rates to (dt-1)
    _compAcquisitionRate = (_dt/(24.0*3600.0)) * _compAcquisitionRate;
    _compLossRate = (_dt/(24.0*3600.0)) * _compLossRate;
    _reefs = reefs;
    _nbOfShallowReefs = nbOfShallowReefs;
  }
  
  /** Calculate connectivity matrix using exposure time method */
  void addConnectivityExposureTime(double t_preComp, double t_postComp) {
    _connectivityExposureTimeFlag = true;
    _settleFlag = false;
    _t_preComp = t_preComp;
    _t_postComp = t_postComp;
    printf("Connectivity by EXPOSURE TIME is enabled. Settling is therefore being disabled.\n");
  }
  
  /** Activate particle swimming and load parameters */
  void addSwimming(fullMatrix<double>* swimParams, dgDofContainer* reefs, dgDofContainer* distToReefs, dgDofContainer* directionToReefs) {
    _swimFlag = true;
    //TEMP: disactivate swimming if swimParams(0,0)=0.0
    if (swimParams < 0.1) _swimFlag=false;
    //Load swimming parameters:
    _swimType = static_cast<int>(swimParams->get(0,0)+0.1);
    _swimSpeed = swimParams->get(0,1);
    _swimDist = swimParams->get(0,2);
    _swimStart = swimParams->get(0,3);
    printf("swimType: %d, swimSpeed: %.2f(m/s), swimDist: %.2f(km), swimStart: %.2f(hrs)\n", _swimType, _swimSpeed, _swimDist, _swimStart);
    _reefs = reefs;
    _distToReefs = distToReefs;
    _directionToReefs = directionToReefs;
  }
  
  /** Method to move a particleArray on by a time-step */
  void particleMove(dgDofContainer* solution, fullMatrix<double>* connectivityMatrix, int n);
  /** Method to check if particle changes competency during time dt */
  int assessCompetency(int particleComp, double t);
  /** Method to check if particle settles or not during time dt */
  int isSettled(double t, int m);
  /** Method to check if particle dies or not during time dt */
  int isDead(double t);
};

#endif