dgParticleTracker2D.h 10.9 KB
Newer Older
thomas's avatar
thomas committed
1
2
3
4
5
6
7
8
#ifndef _DG_PARTICLE_TRACKER_2D_H_
#define _DG_PARTICLE_TRACKER_2D_H_

#include <fstream>
#include <iomanip>
#include <sstream>
#include <iostream>
#include <vector>
9
10
#include"math.h"
#include"GModel.h"
thomas's avatar
thomas committed
11
#include"SPoint3.h"
12
#include"MElement.h"
thomas's avatar
thomas committed
13
#include"dgDofContainer.h"
14
15
16
#include"dgFunctionEvaluator.h"
#include"dgGroupOfElements.h"
#include"dgNeighbourMap.h"
thomas's avatar
thomas committed
17
18
19
20
21
22
23
24

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]);
25
SPoint2 intersection(SPoint2 p1, SPoint2 p2, SPoint2 p3, SPoint2 p4, int& flag);
thomas's avatar
thomas committed
26
27
28

class particle {
private:
thomas's avatar
thomas committed
29
  double _xcoord, _ycoord;
thomas's avatar
thomas committed
30
  int _competency;      //0: pre-competent, 1: competent, 2: post-competent
31
  int _int_state;       //0: alive, 1: settled, 2: dead, 3: out of domain
thomas's avatar
thomas committed
32
33
34
  int _source_reef, _destination_reef;
  //Remember element and group Id particle is in:
  int _groupId, _elementId;
thomas's avatar
thomas committed
35
36
37
38
39
  
public:
  // Default constructor
  particle()
  {
thomas's avatar
thomas committed
40
41
42
43
    _xcoord = 0.0;
    _ycoord = 0.0;
    _int_state = 0;
    _source_reef = -1;
thomas's avatar
thomas committed
44
    _competency = 0;
thomas's avatar
thomas committed
45
46
  }
  
thomas's avatar
thomas committed
47
48
49
50
51
52
53
54
55
  // 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;
thomas's avatar
thomas committed
56
    _competency = 0;
thomas's avatar
thomas committed
57
58
59
  }
  
  // Cartesian constructor which DOESN'T set mesh element Id's
thomas's avatar
thomas committed
60
61
  particle(double x0, double y0, int status)
  {
thomas's avatar
thomas committed
62
63
64
65
66
67
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = -1;
    _groupId = -1;
    _elementId = -1;
thomas's avatar
thomas committed
68
    _competency = 0;
thomas's avatar
thomas committed
69
70
71
72
73
74
75
76
77
78
79
  }
  
  // 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;
thomas's avatar
thomas committed
80
    _competency = 0;
thomas's avatar
thomas committed
81
82
  }
  
thomas's avatar
thomas committed
83
  // Cartesian constructor which sets source reef but NOT mesh element Id's (try not to use this, doesn't store element we're in)
84
85
  particle(double x0, double y0, int srcReef, int status)
  {
thomas's avatar
thomas committed
86
87
88
89
90
91
    _xcoord = x0;
    _ycoord = y0;
    _int_state = status;
    _source_reef = srcReef;
    _groupId = -1;
    _elementId = -1;
thomas's avatar
thomas committed
92
    _competency = 0;
93
94
  }
  
thomas's avatar
thomas committed
95
96
97
  //Access methods:
  double x()
  {
thomas's avatar
thomas committed
98
    return _xcoord;
thomas's avatar
thomas committed
99
100
101
102
  }
  
  double y()
  {
thomas's avatar
thomas committed
103
    return _ycoord;
thomas's avatar
thomas committed
104
105
106
107
  }
  
  int getState()
  {
thomas's avatar
thomas committed
108
109
110
111
112
113
114
115
116
117
118
119
    return _int_state;
  }

  int getSource() 
  {
    return _source_reef;
  }
  
  void getElementId(int *iGroup, int *iElement)
  {
    *iGroup = _groupId;
    *iElement = _elementId;
thomas's avatar
thomas committed
120
121
  }
  
thomas's avatar
thomas committed
122
123
124
125
126
  int getCompetency()
  {
    return _competency;
  }
  
thomas's avatar
thomas committed
127
128
  void setPosition(double x, double y)
  {
thomas's avatar
thomas committed
129
130
    _xcoord=x;
    _ycoord=y;
thomas's avatar
thomas committed
131
132
133
134
  }
  
  void setState(int status)
  {
thomas's avatar
thomas committed
135
136
137
138
139
140
141
    _int_state = status;
  }
  
  void setElementId(int iGroup, int iElement)
  {
    _groupId = iGroup;
    _elementId = iElement;
thomas's avatar
thomas committed
142
  }
thomas's avatar
thomas committed
143
144
145
146
147
148
  
  void setCompetency(int neoComp)
  {
    _competency = neoComp;
  }
  
thomas's avatar
thomas committed
149
150
151
152
};

class particleArray {
private:
thomas's avatar
thomas committed
153
  std::vector<particle> _particle_array;
thomas's avatar
thomas committed
154
  int _alive, _settled, _dead, _outOfDomain, _aliveAndCompetent;
thomas's avatar
thomas committed
155
156
157
158
159
160
161
162
  
public:
  //Constructors:
  particleArray()
  {
    std::cout<<"Error: Default particleArray constructor used.\n";
  }
  
163
  /** Single point constructor which sets element IDs*/
thomas's avatar
thomas committed
164
  particleArray(dgGroupCollection* _group, int M, double x0, double y0, int status, int source_reef)
thomas's avatar
thomas committed
165
  {
thomas's avatar
thomas committed
166
167
168
169
170
171
172
    //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); 
    
thomas's avatar
thomas committed
173
174
175
    for (int i=0; i<M; i++)
    {
      //Seed particle
thomas's avatar
thomas committed
176
      _particle_array.push_back( particle(x0,y0,iGroup,iElem,source_reef,status) );
thomas's avatar
thomas committed
177
178
179
180
    }
    _alive = M;
    _settled = 0;
    _dead = 0;
thomas's avatar
thomas committed
181
    _outOfDomain = 0;
thomas's avatar
thomas committed
182
    _aliveAndCompetent = 0;
thomas's avatar
thomas committed
183
184
  }
  
185
  /**Seed particles using reef map DC*/
thomas's avatar
thomas committed
186
  particleArray(dgGroupCollection* groups, dgDofContainer* reefsDC, fullMatrix<double>* initCM, fullMatrix<double>* reefsToSeed, double concentration, int seedOverWhichReefs, int nbOfShallowReefs);
thomas's avatar
thomas committed
187
    
thomas's avatar
thomas committed
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
  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;

thomas's avatar
thomas committed
204
205
  }
  
thomas's avatar
thomas committed
206
207
208
209
210
211
212
213
  void setStats(int *stats) {
    _alive = stats[0];
    _settled = stats[1];
    _dead = stats[2];
    _outOfDomain = stats[3];
    _aliveAndCompetent = stats[4];
  }
  
214
  void getStats(int *stats) {
thomas's avatar
thomas committed
215
216
217
    stats[0] = _alive;
    stats[1] = _settled;
    stats[2] = _dead;
thomas's avatar
thomas committed
218
    stats[3] = _outOfDomain;
thomas's avatar
thomas committed
219
    stats[4] = _aliveAndCompetent;
thomas's avatar
thomas committed
220
221
222
  }  
  
  int printNbParticles() {
thomas's avatar
thomas committed
223
    return static_cast<int>(_particle_array.size());
224
225
  }
  
thomas's avatar
thomas committed
226
227
228
  //Access a given particle:
  particle* getParticle(int m)
  {
thomas's avatar
thomas committed
229
    return &_particle_array[m];
thomas's avatar
thomas committed
230
231
232
233
  }

  double getParticle_x(int m)
  {
thomas's avatar
thomas committed
234
    return _particle_array[m].x();
thomas's avatar
thomas committed
235
236
237
238
  }
  
  double getParticle_y(int m)
  {
thomas's avatar
thomas committed
239
    return _particle_array[m].y();
thomas's avatar
thomas committed
240
241
  } 
  
thomas's avatar
thomas committed
242
243
244
245
246
  /** 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 
   */
247
  void printPositions(const std::string filename, const std::string filetype, int status);
thomas's avatar
thomas committed
248
  
thomas's avatar
thomas committed
249
250
251
252
  /** Print positions of all particles. 
   *  1 file per particle, colour-coded by time-step.
   *  filetype = dat or pos
   */
253
  void printPositionsFollow(const std::string filename, const std::string filetype, int n, int N, double dt);
thomas's avatar
thomas committed
254
255
256
257
258
  
  /** Print position of one particle. 
   *  Colour-coded by time-step.
   *  filetype = dat or pos
   */
259
  void printPositionsFollowSingle(const std::string filename, const std::string filetype, int n, int N, double dt, int particleId);
thomas's avatar
thomas committed
260
261
};

thomas's avatar
thomas committed
262
263
class connectivityMatrixPrint {
private:
264
  
thomas's avatar
thomas committed
265
266
267
268
269
public:
  /** Constructor */
  connectivityMatrixPrint() {
  }
  /** Method to print contents of matrix */
270
  void CMPrinter(fullMatrix<double>* FM, const std::string filename);
thomas's avatar
thomas committed
271
272
273
274
275
276
};

class dgParticleTracker2D {
private:
  double _dt;
  int _N;
thomas's avatar
thomas committed
277
278
279
  int _settleType, _mortType, _swimType, _settleOverWhichReefType, _nbOfShallowReefs;
  double _t_preComp, _t_postComp, _t_settMax, _settRate_max, _compAcquisitionRate, _compLossRate;
  double _mortRate_max, _mortLambda, _mortNu;
280
  double _swimSpeed, _swimDist, _swimStart;
thomas's avatar
thomas committed
281
282
283
284
285
  dgGroupCollection *_group;
  particleArray *_ParticleArray;
  dgDofContainer* _bathymetry;
  dgDofContainer* _diffusivity;
  dgDofContainer* _reefs;
286
287
  dgDofContainer* _distToReefs;
  dgDofContainer* _directionToReefs;
thomas's avatar
thomas committed
288
289
290
  GModel* _gm;
  dgNeighbourMap* _neighbourMap;
  std::string _LP_OutputDir;
291
  bool _mortFlag, _settleFlag, _swimFlag, _connectivityExposureTimeFlag;
thomas's avatar
thomas committed
292
293
  
public:
294
  /** Constructor */
295
  dgParticleTracker2D(dgGroupCollection* group, particleArray* ParticleArray, dgDofContainer* bathymetry, dgDofContainer* diffusivity,  int N, double dt, const std::string LP_OutputDir) {
thomas's avatar
thomas committed
296
297
298
299
300
301
302
303
    _dt = dt;
    _N = N;
    _group = group;
    _ParticleArray = ParticleArray;
    _bathymetry = bathymetry;
    _diffusivity = diffusivity;
    _LP_OutputDir = LP_OutputDir;
    
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
    //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;
331
332
    _mortType = static_cast<int>(mortParams->get(0,0)+0.1);
    _mortRate_max = mortParams->get(0,1);
thomas's avatar
thomas committed
333
334
    _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)
335
336
337
338
339
  }
  
  /** Activate particle settling and load parameters */
  void addSettling(fullMatrix<double>* settleParams, dgDofContainer* reefs, int nbOfShallowReefs) {
    _settleFlag = true;
340
341
342
343
344
    _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);
thomas's avatar
thomas committed
345
346
347
348
349
350
    _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;
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
    _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;
thomas's avatar
thomas committed
367
368
    //TEMP: disactivate swimming if swimParams(0,0)=0.0
    if (swimParams < 0.1) _swimFlag=false;
369
370
371
372
373
374
    //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);
375
376
377
    _reefs = reefs;
    _distToReefs = distToReefs;
    _directionToReefs = directionToReefs;
thomas's avatar
thomas committed
378
379
  }
  
380
381
  /** Method to move a particleArray on by a time-step */
  void particleMove(dgDofContainer* solution, fullMatrix<double>* connectivityMatrix, int n);
thomas's avatar
thomas committed
382
383
  /** Method to check if particle changes competency during time dt */
  int assessCompetency(int particleComp, double t);
384
  /** Method to check if particle settles or not during time dt */
thomas's avatar
thomas committed
385
  int isSettled(double t, int m);
386
387
  /** Method to check if particle dies or not during time dt */
  int isDead(double t);
thomas's avatar
thomas committed
388
};
thomas's avatar
thomas committed
389
390

#endif