dfnTrans
Code for Particle Tracking simulations in 3D DFN
TrackingPart.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <search.h>
3 #include <stdlib.h>
4 #include <math.h>
5 #include <string.h>
6 #include "FuncDef.h"
7 #include <unistd.h>
8 #include <time.h>
9 #include <sys/stat.h>
10 #include <stdio.h>
11 #include <stdbool.h>
12 
13 struct inpfile {
14  char filename[120];
15  long int flag;
16  double param;
17 };
18 
19 struct lagrangian {
20  double tau;
21  double betta;
22  double initx;
23  double inity;
24  double initz;
25 
26 };
27 
28 struct lagrangian lagvariable;
29 unsigned int FLAG_OUT = 0, all_out = 0;
30 unsigned int np, t, nodeID = 0, avs_o = 0, traj_o = 0, curv_o = 0, no_out = 0, tdrw = 0, mixing_rule = 1;
31 unsigned int marfa = 0, plumec = 0, disp_o = 0, timecounter = 0, frac_o = 0, tfile = 0, tdrw_o = 0, tdrw_limited = 0;
32 double tdrw_porosity = 0.0, tdrw_diffcoeff = 0.0, t_adv0 = 0.0, t_adv = 0.0, timediff = 0.0, tdrw_lambda = 0.0;
33 struct intcoef {
34  double weights[3];
35 };
36 
37 struct posit3d {
38  double cord3[3];
39 };
40 
41 struct tempout
43 {
44  unsigned int times; // count of the time steps
45  double position2d[2]; // 2d position of particles on a fracture plane
46  double position3d[3]; // 3d position of particle in the 3d space
47  double velocity3d[3]; //particles velocity, 3d
48  unsigned int cellp; // number of cell in current particle position
49  unsigned int fracturep; // number of fracture in current particle position
50  double timep; // particle's current travel time
51  double betap; // current beta parameter
52  double length_t; //particle's trajectory length
53  double pressure; // fluid pressure at particle's position
54 };
55 
56 struct tempout *tempdata;
57 
58 static FILE *tmp;
59 static FILE *wpt;
60 static FILE *wpt_att;
61 static FILE *wv;
62 static FILE *wint;
63 static FILE *diff;
64 
66 void ParticleTrack ()
76 {
77  /*** read output options *****/
78  int res;
79  struct inpfile inputfile;
80  char filename[125];
81  unsigned int tort_o = 0;
82  // Output of tortuosity file
83  inputfile = Control_File_Optional("out_tort:", 9);
84 
85  if (inputfile.flag < 0) {
86  tort_o = 0;
87  } else {
88  res = strncmp(inputfile.filename, "yes", 3);
89 
90  if (res == 0) {
91  tort_o = 1;
92  }
93  }
94 
95  // Output of fractures ID
96  inputfile = Control_File_Optional("out_fract:", 10);
97 
98  if (inputfile.flag < 0) {
99  frac_o = 0;
100  } else {
101  res = strncmp(inputfile.filename, "yes", 3);
102 
103  if (res == 0) {
104  frac_o = 1;
105  }
106  }
107 
108  // Determine the subgrid process for particle motion at fracture intersections
109  inputfile = Control_File("streamline_routing:", 19);
110  res = strncmp(inputfile.filename, "yes", 3);
111 
112  if (res == 0) {
113  mixing_rule = 0;
114  printf("\nFracture Intersection Rule: Streamline Routing \n");
115  } else {
116  printf("\nFracture Intersection Rule: Complete Mixing \n");
117  }
118 
119  // Output according to trajectory curvature (not every time step)
120  inputfile = Control_File("out_curv:", 9);
121  res = strncmp(inputfile.filename, "yes", 3);
122 
123  if (res == 0) {
124  curv_o = 1;
125  }
126 
127  // output AVS file for each trajectory
128  inputfile = Control_File("out_avs:", 8);
129  res = strncmp(inputfile.filename, "yes", 3);
130 
131  if (res == 0) {
132  avs_o = 1;
133  }
134 
135  // Output ASCII trajectory file
136  inputfile = Control_File("out_traj:", 9);
137  res = strncmp(inputfile.filename, "yes", 3);
138 
139  if (res == 0) {
140  traj_o = 1;
141  }
142 
143  // Output MARFA input file
144  inputfile = Control_File_Optional("out_marfa:", 10);
145 
146  if (inputfile.flag < 0) {
147  marfa = 0;
148  } else {
149  res = strncmp(inputfile.filename, "yes", 3);
150 
151  if (res == 0) {
152  marfa = 1;
153  traj_o = 1;
154  curv_o = 1;
155  }
156  }
157 
158  // Output PLUMECALC input file
159  inputfile = Control_File_Optional("out_plumecalc:", 14);
160 
161  if (inputfile.flag < 0) {
162  plumec = 0;
163  } else {
164  res = strncmp(inputfile.filename, "yes", 3);
165 
166  if (res == 0) {
167  plumec = 1;
168  traj_o = 1;
169  curv_o = 1;
170  }
171  }
172 
173  // Output Time control planes
174  inputfile = Control_File_Optional("out_disp:", 9);
175 
176  if (inputfile.flag < 0) {
177  disp_o = 0;
178  } else {
179  res = strncmp(inputfile.filename, "yes", 3);
180 
181  if (res == 0 ) {
182  disp_o = 1;
183  }
184  }
185 
186  //if there is an additional output
187  if ((avs_o + traj_o) == 0) {
188  no_out = 1;
189  }
190 
191  // Add diffusion time (TDRW)
192  inputfile = Control_File_Optional("tdrw:", 5);
193 
194  if (inputfile.flag < 0) {
195  tdrw = 0;
196  } else {
197  res = strncmp(inputfile.filename, "yes", 3);
198 
199  if (res == 0) {
200  printf("\nRunning with Time Domain Random Walk\n");
201  tdrw = 1;
202  inputfile = Control_Param("tdrw_porosity:", 14);
203  tdrw_porosity = inputfile.param;
204  inputfile = Control_Param("tdrw_diffcoeff:", 15);
205  tdrw_diffcoeff = inputfile.param;
206  inputfile = Control_File_Optional("tdrw_rate_limited:", 5);
207  res = strncmp(inputfile.filename, "yes", 3);
208 
209  if (res == 0) {
210  tdrw_limited = true;
211  printf("--> Running with Rate Limited TDRW\n");
212  inputfile = Control_Param("tdrw_lambda:", 14);
213  tdrw_lambda = inputfile.param;
214  printf("--> Rate Limited Distance: %0.2f m\n\n", tdrw_lambda);
215  }
216 
217  printf("--> Matrix Porosity: %0.2f\n", tdrw_porosity );
218  printf("--> Matrix Diffusivity: %0.2e m^2/s\n\n", tdrw_diffcoeff);
219  tdrw_diffcoeff = tdrw_diffcoeff * timeunit;
220  }
221  }
222 
223  // open file with traveling time results of all particles
224  inputfile = Control_File("out_time:", 9 );
225  sprintf(filename, "%s/%s", maindir, inputfile.filename);
226  FILE *tp = OpenFile (filename, "w");
227 
228  if (tdrw == 1) {
229  fprintf(tp, "# of time steps, flux weights, total advective travel time, total advective + diffusion time, total diffusion time, beta, total length[m] \n");
230  } else {
231  fprintf(tp, "# of time steps, flux weights, total travel time, x-, y-, z-final pos, beta, total length[m] \n");
232  }
233 
234  /* output of initial and final positions of particle*/
235  FILE *inp;
236  FILE *fnp;
237  inputfile = Control_File_Optional("allparticles_output:", 20);
238 
239  if (inputfile.flag < 0) {
240  all_out = 0;
241  } else {
242  res = strncmp(inputfile.filename, "yes", 3);
243 
244  if (res == 0) {
245  all_out = 1;
246  sprintf(filename, "%s/initpos", maindir);
247  inp = OpenFile (filename, "w");
248  fprintf(inp, " #of particle, # of cell, # of fracture, x-, y-, z- initial position, flux weight ");
249  sprintf(filename, "%s/finpos", maindir);
250  fnp = OpenFile (filename, "w");
251  fprintf(fnp, " #of particle, # of cell, # of fracture, x-, y-, z- final position, flux weight ");
252  }
253  }
254 
255  // open tortuosity file
256  char path[125];
257  FILE *tort;
258 
259  if (tort_o > 0) {
260  sprintf(filename, "%s/torts.dat", maindir);
261  tort = OpenFile (filename, "w");
262  fprintf(tort, "Data for tortuosity calculation: total length of trajectory, x-,y-, z- of initial pos., x-, y-, z- final pos., number of intersections. \n");
263  }
264 
265  // open FractureID file
266  FILE *frac;
267 
268  if (frac_o > 0) {
269  sprintf(filename, "%s/FractureID", maindir);
270  frac = OpenFile (filename, "w");
271  fprintf(frac, "Fractures ID \n");
272  }
273 
274  // output in case of tdrw
275  if (tdrw == 1) {
276  inputfile = Control_File_Optional("tdrw_out:", 9);
277 
278  if (inputfile.flag < 0) {
279  tdrw_o = 0;
280  } else {
281  res = strncmp(inputfile.filename, "yes", 3);
282 
283  if (res == 0) {
284  tdrw_o = 1;
285  }
286  }
287  }
288 
289  // Create path for trajectory outputs
290  if ((no_out != 1) || (tdrw_o == 1)) {
291  inputfile = Control_File("out_path:", 9 );
292  sprintf(path, "%s/%s", maindir, inputfile.filename );
293  mkdir(path, 0777);
294  printf("\n All output trajectory files will be written in %s/ \n", path);
295  printf(" Note: if the directory exists all the files in it will be replaced. \n \n");
296  // temporary otputs go to external file or memory buffer
297  inputfile = Control_File_Optional("out_filetemp:", 13);
298 
299  if (inputfile.flag < 0) {
300  tfile = 0;
301  } else {
302  res = strncmp(inputfile.filename, "yes", 3);
303 
304  if (res == 0) {
305  tfile = 1;
306  }
307  }
308  }
309 
310  //settings for Control Plane/Cylinder Output
311  int out_control = 0, out_plane = 0, out_cylinder = 0;
312  inputfile = Control_File("ControlPlane:", 13);
313  res = String_Compare(inputfile.filename, "yes");
314 
315  if (res == 0) {
316  out_plane = 1;
317  } else {
318  inputfile = Control_File_Optional("ControlCylinder:", 16);
319 
320  if (inputfile.flag < 0) {
321  out_cylinder = 0;
322  } else {
323  res = strncmp(inputfile.filename, "yes", 3);
324 
325  if (res == 0) {
326  out_cylinder = 1;
327  }
328  }
329  }
330 
331  if ((out_cylinder + out_plane) > 0) {
332  out_control = 1;
333  }
334 
335  // reading variables for dispersivity (time-control) calculation
336  unsigned int time_d = 1, kd = 1;
337  double dtime = 0.0;
338  double epsl = 0.0;
339 
340  if (disp_o == 1) {
341  inputfile = Control_Param("out_dtimest:", 12 );
342  time_d = (int)inputfile.param;
343  inputfile = Control_Param("out_dtime:", 10 );
344  dtime = inputfile.param;
345  epsl = 0.05 * dtime;
346  }
347 
348  double sumsquares[time_d][2];//keeps sum of squares, like (x0-x)^2 for all particles
349  FILE * dis;
350  unsigned int nd[time_d];
351  int itime;
352  char filetime[15];
353  double inflowcoord = 0.0; //define automatically
354  double outflowcoord = 0.0;
355  double controllength = 0.0;
356  double wellthick = 0.0;
357  char pathcontrol[125];
358  double deltaCP = 0, current_CP = 0.0;
359  int ic, icl = 0, idist = 0, flowd = 0, welld = 0;
360 
361  if (disp_o == 1) {
362  for (ic = 0; ic < time_d - 1; ic++) {
363  sumsquares[ic][0] = 0.0;
364  sumsquares[ic][1] = 0.0;
365  nd[ic] = 0;
366  sprintf(filename, "%s/ControlTime_t%d", maindir, ic + 1);
367  dis = OpenFile(filename, "w");
368  // fprintf(dis,"positions for longitudinal dispersivity calculation at time %f \n", (ic+1)*dtime);
369  fclose(dis);
370  }
371 
372  inputfile = Control_Param("dflowdir:", 8);
373  flowd = inputfile.param;
374  }
375 
376  // settings for Control Plane / Control Cylinder options
377  if (out_control == 1) {
378  inputfile = Control_File("control_out:", 12 );
379  sprintf(pathcontrol, "%s/%s", maindir, inputfile.filename );
380  mkdir(pathcontrol, 0777);
381  inputfile = Control_Param("delta_Control:", 14);
382  deltaCP = inputfile.param;
383 
384  if (out_plane == 1) {
385  inputfile = Control_Param("flowdir:", 8);
386  flowd = inputfile.param;
387  node[nodezonein[0] - 1].fracture[0];
388  inflowcoord = node[nodezonein[0] - 1].coord[flowd];
389  outflowcoord = node[nodezoneout[0] - 1].coord[flowd];
390  controllength = fabs(outflowcoord) + fabs(inflowcoord);
391  icl = controllength / deltaCP + 1;
392  }
393 
394  if (out_cylinder == 1) {
395  inputfile = Control_Param("welldir:", 8);
396  welld = inputfile.param;
397  inputfile = Control_Param("lengthtowell:", 13);
398  controllength = inputfile.param;
399  inputfile = Control_Param("wellthickness:", 14);
400  wellthick = inputfile.param;
401 
402  if (deltaCP < wellthick) {
403  icl = (controllength - wellthick / 2) / deltaCP + 1;
404  } else {
405  icl = controllength / deltaCP + 1;
406  }
407  }
408  }
409 
410  int cross[icl];
411  /*** define particle's initial positions **/
412  int numbpart, initweight = 0;
413  /*** set up initial positions of particles ***/
414  numbpart = InitPos();
415  printf("\n %ld time steps in the current run \n", timesteps);
416  printf("\n Initial number of particles being placed: %d \n", numbpart);
417  inputfile = Control_File("flux_weight:", 12);
418  res = strncmp(inputfile.filename, "yes", 3);
419 
420  if ((res == 0) && (flag_w == 1)) {
421  initweight = 1;
422  }
423 
424  int i = 1, ins = 0, intersm = 0, curr_n = 1, curr_o = 1;
425  int t_end = 0, fracthit = 0;
426  struct posit3d particle3dposit, particle3dvelocity;
427  int counttimestep = 0, prevcell = 0, prevfract = 0;
428  double xcop, ycop, zcop, currentlength = 0, totallength = 0;
429  double xx = 0, yy = 0, zz = 0, xinit = 0, yinit = 0, zinit = 0;
430 
431  /**** calculate initial flux weight of particles *****/
432  /**** works for first 3 options of particles initial positions ******/
433 
434  if (initweight == 1) {
435  FlowInWeight(numbpart);
436  }
437 
438  unsigned int percent_done = 0;
439 
440  /************ LOOP ON PARTICLES **********/
441  for (np = 0; np < numbpart; np++) {
442  t = 0;
443 
444  if ((np >= (int)(0.01 * numbpart)) && ( percent_done == 0)) {
445  printf("Done %d particles, 1%%. \n", np);
446  percent_done = 1;
447  }
448 
449  if ((np >= (int)(0.05 * numbpart)) && ( percent_done == 1)) {
450  printf("Done %d particles, 5%%. \n", np);
451  percent_done = 2;
452  }
453 
454  if ((np >= (int)(0.25 * numbpart)) && ( percent_done == 2)) {
455  printf("Done %d particles, 25%%. \n", np);
456  percent_done = 3;
457  }
458 
459  if ((np >= (int)(0.5 * numbpart)) && ( percent_done == 3)) {
460  printf("Done %d particles, 50%%. \n", np);
461  percent_done = 4;
462  }
463 
464  if ((np >= (int)(0.75 * numbpart)) && ( percent_done == 4)) {
465  printf("Done %d particles, 75%%. \n", np);
466  percent_done = 5;
467  }
468 
469  if (avs_o == 1) {
470  // AVS output (should be optional)
471  sprintf(filename, "%s/part3D_%d.inp", path, curr_n);
472  wpt = OpenFile(filename, "w");
473  fprintf(wpt, "%10lu %10d %10d %10d %10d\n", timesteps, 0, 0, 0, 0);
474  //open a separate file for attributes, will be attached to the original AVS later
475  sprintf(filename, "%s/part3D_%d.att", path, curr_n);
476  wpt_att = OpenFile(filename, "w");
477  fprintf(wpt_att, "0008 1 1 1 1 1 1 1 1\n");
478  fprintf(wpt_att, "fracture, integer\n");
479  fprintf(wpt_att, "time, real\n");
480  fprintf(wpt_att, "velocity, real\n");
481  fprintf(wpt_att, "vel_x, real\n");
482  fprintf(wpt_att, "vel_y, real\n");
483  fprintf(wpt_att, "vel_z, real\n");
484  fprintf(wpt_att, "aperture, real\n");
485  fprintf(wpt_att, "pressure, real\n");
486  }
487 
488  if (traj_o == 1) {
489  // ascii output of: 3d positions, 3d velocities, cell, fracture, time and beta
490  sprintf(filename, "%s/traject_%d", path, curr_n);
491  wv = OpenFile(filename, "w");
492  fprintf(wv, " Current time step, x-, y-, z- pos., Vx, Vy, Vz at this positions, # of cell, #of fracture, travel time, aperture, beta, intersect. fracture ID, fluid pressure at particle's position \n");
493  // output data on intersections only
494  sprintf(filename, "%s/inters_%d", path, curr_n);
495  wint = OpenFile(filename, "w");
496  fprintf(wint, " Current traj. length, travel time, x-, y-, z- pos., fracture ID, beta, fluid pressure at part.pos. \n");
497  }
498 
499  if ((tdrw == 1) && (tdrw_o == 1)) {
500  sprintf(filename, "%s/tdrw_%d", path, curr_n);
501  diff = OpenFile(filename, "w");
502  fprintf(diff, " Advective travel time on the fracture, Diffusion time on the fracture, Total travel time on the fracture, fracture ID, Accumulative advective travel time, Accumulative total time, Accumulative diffusion time \n");
503  }
504 
505  // define capacity for temp data used for outputs
506  int capacity = (int) timesteps / 10;
507 
508  if (disp_o = ! 1) {
509  time_d = 1;
510  }
511 
512  double part_squares[time_d][3];
513 
514  if (disp_o = ! 1) {
515  for (ic = 0; ic < time_d; ic++) {
516  part_squares[ic][0] = 0.0;
517  part_squares[ic][1] = 0.0;
518  part_squares[ic][2] = 0.0;
519  }
520  }
521 
522  // control plane/cylinder output
523  t_end = 0;
524  intersm = 0;
525  /* define an initial cell */
526  ins = 0;
527 
528  if (particle[np].cell != 0) {
529  ins = 1;
530  } else {
531  ins = InitCell();
532  }
533 
534  if (ins == 0) {
535  printf("Initial cell is not found for particle %d %f %f in fract %d. \n", np + 1, particle[np].position[0], particle[np].position[1], particle[np].fracture);
536  } else {
537  // set up initial values for Lagrangian variables
538  lagvariable.tau = 0.0;
539  lagvariable.betta = 0.0;
540  prevcell = particle[np].cell;
541  prevfract = particle[np].fracture;
542  unsigned int fract_id[nfract], id = 0;
543 
544  for (id = 0; id < nfract; id++) {
545  fract_id[id] = 0;
546  }
547 
548  fract_id[0] = particle[np].fracture;
549  FLAG_OUT = 0;
550  int stuck = 0, stuckcell = 0, cur_ind = 0, cur_node = 0;
551  t = 0;
552  nodeID = 0;
553  counttimestep = 0;
554  xinit = 0.0;
555  yinit = 0.0;
556  zinit = 0.0;
557  fracthit = 0;
558  // output particles initial positions (in 3D)
559  particle3dposit = CalculatePosition3D();
560  xcop = particle3dposit.cord3[0];
561  ycop = particle3dposit.cord3[1];
562  zcop = particle3dposit.cord3[2];
563  xinit = xcop;
564  yinit = ycop;
565  zinit = zcop;
566  lagvariable.initx = xinit;
567  lagvariable.inity = yinit;
568  lagvariable.initz = zinit;
569  kd = 1;
570  double tautau = 0.0;
571  double beta = 0.0;
572 
573  if (all_out > 0) {
574  fprintf(inp, "\n %d %d %d %5.12E %5.12E %5.12E %5.12E", np + 1, particle[np].cell, particle[np].fracture, particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle[np].fl_weight);
575  }
576 
577  totallength = 0.0;
578  currentlength = 0.0;
579  //counts for control plane/time output
580  FILE* tmp2;
581 
582  if (out_control == 1) {
583  for (ic = 0; ic < icl; ic++) {
584  cross[ic] = 0;
585  }
586 
587  idist = 0;
588  sprintf(filename, "%s/part_control_%d", pathcontrol, curr_n);
589  tmp2 = OpenFile(filename, "w");
590 
591  if (tdrw == 1) {
592  fprintf(tmp2, " x-, y-, z- position, Vx, Vy, Vz, trajectory length, fracture ID , aperture, Accumulative advective travel time, Accumulative total time, Accumulative diffusion time \n");
593  } else {
594  fprintf(tmp2, " travel time, x-, y-, z- position, Vx, Vy, Vz, trajectory length, # of current fracture, aperture \n");
595  }
596 
597  if (out_plane == 1) {
598  if (tdrw == 1) {
599  fprintf(tmp2, " %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %5.12E %5.12E %5.12E %5.12E \n", particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], totallength, particle[np].fracture, node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture, particle[np].time, particle[np].t_adv_diff, particle[np].t_diff);
600  } else {
601  fprintf(tmp2, "%5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %5.12E\n", particle[np].time, particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], totallength, particle[np].fracture, node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture);
602  }
603 
604  if (inflowcoord < 0) {
605  current_CP = inflowcoord + deltaCP;
606  } else {
607  current_CP = inflowcoord - deltaCP;
608  }
609  }
610 
611  if (out_cylinder == 1) {
612  current_CP = controllength;
613  }
614  }
615 
616  if (tfile > 0) {
617  sprintf(filename, "%s/tempdata_%d", maindir, np);
618  tmp = OpenFile(filename, "w");
619  }
620 
621  timecounter = 0; //for temp data allocation
622  t_adv0 = 0; //for tdrw calculation; starting time
623  particle[np].t_diff = 0.0;
624  particle[np].t_adv_diff = 0.0;
625 
627 
628  for (t = 0; t < timesteps; t++) {
629  particle3dposit = CalculatePosition3D();
630  xx = particle3dposit.cord3[0] - xcop;
631  yy = particle3dposit.cord3[1] - ycop;
632  zz = particle3dposit.cord3[2] - zcop;
633  currentlength = sqrt(xx * xx + yy * yy + zz * zz);
634  totallength = totallength + currentlength;
635  xcop = particle3dposit.cord3[0];
636  ycop = particle3dposit.cord3[1];
637  zcop = particle3dposit.cord3[2];
638 
639  if (no_out != 1) {
640  particle3dvelocity = CalculateVelocity3D();
641 
642  if (tfile == 1) {
643  fprintf(tmp, "%05d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %05d %5.12E %5.12E %5.12E %5.12E \n", t, particle[np].position[0], particle[np].position[1], particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], particle[np].cell, particle[np].fracture, particle[np].time, beta, totallength, particle[np].pressure);
644  } else {
645  if (timecounter == 0) {
646  tempdata = (struct tempout*) malloc (capacity * sizeof(struct tempout));
647  }
648 
649  if (tempdata == NULL) {
650  printf("Allocation memory problem - tempdata\n");
651  printf("timecounter %d time %d capacity %d\n", timecounter, t, capacity);
652  }
653 
654  particle3dvelocity = CalculateVelocity3D();
655  tempdata[timecounter].times = t;
656  tempdata[timecounter].position2d[0] = particle[np].position[0];
657  tempdata[timecounter].position2d[1] = particle[np].position[1];
658  tempdata[timecounter].position3d[0] = particle3dposit.cord3[0];
659  tempdata[timecounter].position3d[1] = particle3dposit.cord3[1];
660  tempdata[timecounter].position3d[2] = particle3dposit.cord3[2];
661  tempdata[timecounter].velocity3d[0] = particle3dvelocity.cord3[0];
662  tempdata[timecounter].velocity3d[1] = particle3dvelocity.cord3[1];
663  tempdata[timecounter].velocity3d[2] = particle3dvelocity.cord3[2];
664  tempdata[timecounter].cellp = particle[np].cell;
665  tempdata[timecounter].fracturep = particle[np].fracture;
666  tempdata[timecounter].timep = particle[np].time;
667  tempdata[timecounter].betap = beta;
668  tempdata[timecounter].length_t = totallength;
669  tempdata[timecounter].pressure = particle[np].pressure;
670  timecounter++;
671 
672  // if memory should be reallocated
673  if (timecounter == capacity) {
674  capacity = 2 * capacity;
675 
676  if (capacity >= timesteps) {
677  printf("overload\n");
678  FLAG_OUT = 0;
679  t_end = t;
680  break;
681  free(tempdata);
682  }
683 
684  tempdata = (struct tempout*)realloc(tempdata, sizeof(struct tempout) * capacity);
685 
686  if (tempdata == NULL) {
687  printf("REAllocation memory problem - tempdata\n");
688  printf("timecounter %d time %d capacity %d\n", timecounter, t, capacity);
689  }
690  }
691  }
692  }
693 
694  //calculations for dispersivity: square of (xo-x) is calculated for transverse disersivity only
695  // for longitudinal dispersivity we save the actual coordination of the particle (commented out)
696  double ctime = 0.0;
697 
698  if (disp_o == 1) {
699  ctime = kd * dtime;
700 
701  if (((ctime - epsl) <= particle[np].time) && ((ctime + epsl) >= particle[np].time)) {
702  // 3d position of particle
703  part_squares[kd - 1][0] = xcop;
704  part_squares[kd - 1][1] = ycop;
705  part_squares[kd - 1][2] = zcop;
706 
707  if (kd < time_d) {
708  kd++;
709  }
710  }
711  }
712 
713  /***** output data at each control plane ********/
714  double welldist = 0.0; //shortest distance from particle to well
715 
716  if (out_control == 1) {
717  if (out_plane == 1) {
718  if (( (inflowcoord < 0) && (cross[idist] == 0) && (particle3dposit.cord3[flowd] >= current_CP)) || ( (inflowcoord > 0) && (cross[idist] == 0) && (particle3dposit.cord3[flowd] <= current_CP))) {
719  cross[idist] = 1;
720 
721  if (no_out == 1) {
722  particle3dvelocity = CalculateVelocity3D();
723  }
724 
725  if (tdrw == 1) {
726  t_adv = particle[np].time - t_adv0;
727  timediff = TimeDomainRW(t_adv);
728  t_adv0 = particle[np].time;
729  particle[np].t_diff = particle[np].t_diff + timediff;
730  particle[np].t_adv_diff = particle[np].t_adv_diff + t_adv + timediff;
731  fprintf(tmp2, "%5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %5.12E %5.12E %5.12E %5.12E \n", particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], totallength, particle[np].fracture, node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture, particle[np].time, particle[np].t_adv_diff, particle[np].t_diff);
732  } else {
733  fprintf(tmp2, "%5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %5.12E\n", particle[np].time, particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], totallength, particle[np].fracture, node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture);
734  }
735 
736  if (inflowcoord < 0) {
737  current_CP = current_CP + deltaCP;
738  } else {
739  current_CP = current_CP - deltaCP;
740  }
741 
742  idist = idist + 1;
743  }
744  }
745 
746  if ((out_cylinder == 1) && (current_CP >= wellthick / 2)) {
747  if (welld == 0) {
748  welldist = sqrt(ycop * ycop + zcop * zcop);
749  }
750 
751  if (welld == 1) {
752  welldist = sqrt(xcop * xcop + zcop * zcop);
753  }
754 
755  if (welld == 2) {
756  welldist = sqrt(xcop * xcop + ycop * ycop);
757  }
758 
759  if ((cross[idist] == 0) && (welldist <= current_CP)) {
760  cross[idist] = 1;
761 
762  if (no_out == 1) {
763  particle3dvelocity = CalculateVelocity3D();
764  }
765 
766  fprintf(tmp2, "%5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %5.12E\n", particle[np].time, xcop, ycop, zcop, particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], totallength, particle[np].fracture, node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture);
767  current_CP = current_CP - deltaCP;
768  idist = idist + 1;
769  }
770  }
771  }
772 
773  /*** if particle is on intersection cell **/
774  /*** check distance to intersection ***/
775 
776  if ((particle[np].intcell == 1) || (particle[np].intcell == 3)) {
777  intersm = CheckDistance ();
778  prevcell = particle[np].cell;
779  }
780 
781  /*** if particle's new cell was not found move to the next particle ***/
782  if (particle[np].cell == 0) {
783  break;
784  }
785 
786  /*** Get new particle's velocity ***/
787  PredictorStep();
788  /*** Calculate new weights and check if particle is in new cell ***/
789  CheckNewCell();
790 
791  if ((FLAG_OUT == 1)) {
792  break;
793  }
794 
795  if (particle[np].cell != 0) {
796  if ((particle[np].intcell != 1) && (particle[np].intcell != 3))
797  /*** Get new particle's position ***/
798  {
799  CorrectorStep();
800  }
801 
802  lagvariable = CalculateLagrangian(xcop, ycop, zcop, lagvariable.initx, lagvariable.inity, lagvariable.initz);
803  tautau = tautau + lagvariable.tau;
804  beta = beta + lagvariable.betta;
805  particle[np].time = particle[np].time + lagvariable.tau;
806  } else {
807  if ((FLAG_OUT != 1) && (FLAG_OUT != 3)) {
808  // printf("%5.8e % 5.8e %d \n", particle[np].velocity[0], particle[np].velocity[1], particle[np].cell);
809  // printf("cell=0 %d for particle %d at time %d after predictor step \n",particle[np].cell, np+1, t );
810  } else {
811  if (FLAG_OUT == 1) {
812  break;
813  } else if ((particle[np].cell == 0) && (FLAG_OUT == 3)) {
814  Moving2NextCellBound(prevcell);
815  FLAG_OUT = 0;
816  }
817  }
818  }
819 
820  /*** Calculate new weights and check if particle is in new cell ***/
821  if (particle[np].cell != 0) {
822  CheckNewCell();
823  }
824 
825  if ((particle[np].cell == 0) && (FLAG_OUT == 3)) {
826  FLAG_OUT = 0;
827  Moving2NextCellBound(prevcell);
828  }
829 
830  if ((particle[np].cell != prevcell) && (particle[np].cell != 0)) {
831  counttimestep = 0;
832  prevcell = particle[np].cell;
833 
834  if (prevfract != particle[np].fracture) {
835  stuck = 0;
836  stuckcell = 0;
837  cur_ind = 0;
838  cur_node = 0;
839  prevfract = particle[np].fracture;
840  fracthit = fracthit + 1;
841 
842  if (fracthit > nfract) {
843  FLAG_OUT = 0;
844  break;
845  }
846 
847  fract_id[fracthit] = particle[np].fracture;
848  // printf(" %d %d %15.8e \n",np+1, prevfract, particle[np].time);
849  }
850  } else if (particle[np].cell == prevcell) {
851  counttimestep++;
852 
853  /**** The number of time steps that particle spent in one cell is counted.
854  If particles spent more than 1/10 of the total time steps given by user for the entire path, then the particle is considered to be stuck and the time lopp ends. ****/
855  if (counttimestep > (int)timesteps / 10.0) {
856  // printf("stuck %05d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %05d %5.12E %d %d\n",t+1,particle[np].position[0], particle[np].position[1], particle3dposit.cord3[0], particle3dposit.cord3[1],particle3dposit.cord3[2],particle3dvelocity.cord3[0], particle3dvelocity.cord3[1],particle3dvelocity.cord3[2],particle[np].cell, particle[np].fracture, particle[np].time, counttimestep, (int)timesteps/3.0);
857  FLAG_OUT = 0;
858  break;
859  }
860  }
861 
862  t_end = t;
863 
864  /*** if particle's new cell was not found move to the next particle ***/
865  if (particle[np].cell == 0) {
866  FLAG_OUT = 0;
867  break;
868  }
869  } //end of loop on time steps
870 
871  /********** Final outputs ***************/
872 
873  /***** if particle did not go out through flow-out zone ****/
874  if (FLAG_OUT == 0) {
875  if (all_out == 0) {
876  if (tfile == 1) {
877  fclose(tmp);
878  int status;
879  sprintf(filename, "%s/tempdata_%d", maindir, np);
880  status = remove(filename);
881  }
882 
883  if (out_control == 1) {
884  fclose(tmp2);
885  }
886  } else {
887  curr_o++;
888  }
889  }
890 
891  /**** if particle went out through flow out zone ****/
892  if ((FLAG_OUT == 1) || (all_out == 1)) {
893  if (no_out != 1) {
894  if (particle[np].cell != 0) {
895  FinalPosition();
896  }
897 
898  particle3dposit = CalculatePosition3D();
899  particle3dvelocity = CalculateVelocity3D();
900 
901  if (tfile == 1) {
902  fprintf(tmp, "%05d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %05d %5.12E %5.12E %5.12E %5.12E \n", t, particle[np].position[0], particle[np].position[1], particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], particle[np].cell, particle[np].fracture, particle[np].time, beta, totallength, particle[np].pressure);
903  } else {
904  tempdata[timecounter].times = t;
905  tempdata[timecounter].position2d[0] = particle[np].position[0];
906  tempdata[timecounter].position2d[1] = particle[np].position[1];
907  tempdata[timecounter].position3d[0] = particle3dposit.cord3[0];
908  tempdata[timecounter].position3d[1] = particle3dposit.cord3[1];
909  tempdata[timecounter].position3d[2] = particle3dposit.cord3[2];
910  tempdata[timecounter].velocity3d[0] = particle3dvelocity.cord3[0];
911  tempdata[timecounter].velocity3d[1] = particle3dvelocity.cord3[1];
912  tempdata[timecounter].velocity3d[2] = particle3dvelocity.cord3[2];
913  tempdata[timecounter].cellp = particle[np].cell;
914  tempdata[timecounter].fracturep = particle[np].fracture;
915  tempdata[timecounter].timep = particle[np].time;
916  tempdata[timecounter].betap = beta;
917  tempdata[timecounter].length_t = totallength;
918  tempdata[timecounter].pressure = particle[np].pressure;
919  }
920 
921  ParticleOutput(t, 0);
922 
923  if (tfile == 1) {
924  fclose(tmp);
925  int status;
926  sprintf(filename, "%s/tempdata_%d", maindir, np);
927  status = remove(filename);
928  }
929  } else {
930  if (particle[np].cell != 0) {
931  FinalPosition();
932  }
933 
934  particle3dposit = CalculatePosition3D();
935  }
936 
937  if (tdrw == 1) {
938  t_adv = particle[np].time - t_adv0;
939  // printf("%d %lf %lf %lf \n", np+1, t_adv, particle[np].time, t_adv0);
940  timediff = TimeDomainRW(t_adv);
941  t_adv0 = particle[np].time;
942  particle[np].t_diff = particle[np].t_diff + timediff;
943  particle[np].t_adv_diff = particle[np].t_adv_diff + t_adv + timediff;
944 
945  if (tdrw_o == 1) {
946  fprintf(diff, "%5.12E %5.12E %5.12E %d %5.12E %5.12E %5.12E \n", t_adv, timediff, timediff + t_adv, particle[np].fracture, particle[np].time, particle[np].t_adv_diff, particle[np].t_diff);
947  }
948  }
949 
950  //adding data to dispersivity
951 
952  for (ic = 0; ic < kd - 1; ic++) {
953  nd[ic]++;
954 
955  if (disp_o == 1) {
956  sprintf(filename, "%s/dispers_t%d", maindir, ic + 1);
957  dis = OpenFile(filename, "r+");
958  fseek(dis, 0, SEEK_END);
959  fprintf(dis, "%f %f %f \n", part_squares[ic][0], part_squares[ic][1], part_squares[ic][2]);
960  fclose(dis);
961  }
962  }
963 
964  if (particle[np].cell != 0) {
965  lagvariable = CalculateLagrangian(xcop, ycop, zcop, lagvariable.initx, lagvariable.inity, lagvariable.initz);
966  tautau = tautau + lagvariable.tau;
967  beta = beta + lagvariable.betta;
968  }
969 
970  xx = particle3dposit.cord3[0] - xcop;
971  yy = particle3dposit.cord3[1] - ycop;
972  zz = particle3dposit.cord3[2] - zcop;
973  currentlength = sqrt(xx * xx + yy * yy + zz * zz);
974  totallength = totallength + currentlength;
975  sprintf(filename, "%s/initpos", maindir);
976 
977  if (all_out == 1) {
978  if (particle[np].cell != 0) {
979  fprintf(fnp, "\n %d %d %d %5.12E %5.12E %5.12E %5.12E ", np + 1, particle[np].cell, particle[np].fracture, particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle[np].fl_weight);
980  } else {
981  fprintf(fnp, "\n %d %d %d ", np + 1, particle[np].cell, particle[np].fracture);
982  }
983  }
984 
985  /******* output travel time *****/
986  if (tdrw == 1) {
987  fprintf(tp, "%d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E \n", t_end, particle[np].fl_weight, particle[np].time, particle[np].t_adv_diff, particle[np].t_diff, beta, totallength);
988  } else {
989  fprintf(tp, "%d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E \n", t_end, particle[np].fl_weight, particle[np].time, particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], beta, totallength);
990  }
991 
992  if (tort_o > 0) {
993  fprintf(tort, "%5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %d\n", totallength, xinit, yinit, zinit, particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], fracthit);
994  }
995 
996  if (frac_o > 0) {
997  id = 0;
998 
999  do {
1000  fprintf(frac, "%d ", fract_id[id]);
1001  id++;
1002  } while (fract_id[id] != 0);
1003 
1004  fprintf(frac, "\n");
1005  }
1006 
1007  // output of last control plane - outflow plane
1008 
1009  if (out_control == 1) {
1010  if (particle[np].cell == 0) {
1011  particle[np].cell = prevcell;
1012  }
1013 
1014  particle3dvelocity = CalculateVelocity3D();
1015 
1016  if (tdrw == 1) {
1017  t_adv = particle[np].time - t_adv0;
1018  timediff = TimeDomainRW(t_adv);
1019  t_adv0 = particle[np].time;
1020  particle[np].t_diff = particle[np].t_diff + timediff;
1021  particle[np].t_adv_diff = particle[np].t_adv_diff + t_adv + timediff;
1022  fprintf(tmp2, "%5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %5.12E %5.12E %5.12E %5.12E \n", particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], totallength, particle[np].fracture, node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture, particle[np].time, particle[np].t_adv_diff, particle[np].t_diff);
1023  } else {
1024  fprintf(tmp2, "%5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %5.12E\n", particle[np].time, particle3dposit.cord3[0], particle3dposit.cord3[1], particle3dposit.cord3[2], particle3dvelocity.cord3[0], particle3dvelocity.cord3[1], particle3dvelocity.cord3[2], totallength, particle[np].fracture, node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture);
1025  }
1026 
1027  fclose(tmp2);
1028  }
1029 
1030  curr_n++;
1031 
1032  if (avs_o == 1) {
1033  /*** write a connectivity list in inp files ***/
1034  for (i = 0; i < nodeID - 1; i++) {
1035  fprintf(wpt, "%10d %5d line %10d %10d\n", i + 1, 1, i + 1, i + 2);
1036  }
1037 
1038  rewind(wpt);
1039  fprintf(wpt, "%10d %10d %10d %10d %10d\n", nodeID, nodeID - 1, 8, 0, 0);
1040  fclose(wpt);
1041  fclose(wpt_att);
1042  char filename1[125] = {0}, filename2[125] = {0}, filename3[125] = {0}, buffer[500] = {0};
1043  char wspace = ' ';
1044  sprintf(filename1, "%s/part3D_%d.inp", path, curr_n - 1);
1045  sprintf(filename2, "%s/part3D_%d.att", path, curr_n - 1);
1046  sprintf(filename3, "%s/part_%d.inp", path, curr_n - 1);
1047  sprintf(buffer, "cat%c%s%c%s%c>%c%s", wspace, filename1, wspace, filename2, wspace, wspace, filename3);
1048  system (buffer);
1049  sprintf(buffer, "rm%c-f%c%s%c%s ", wspace, wspace, filename1, wspace, filename2);
1050  system (buffer);
1051  }
1052  }
1053  } //end if ins!=0 (the initial cell was found)
1054 
1055  if (traj_o == 1) {
1056  rewind(wv);
1057  fprintf(wv, "%d \n", nodeID);
1058  fclose(wv);
1059  fclose(wint);
1060  }
1061 
1062  if ((tdrw == 1) && (tdrw_o == 1)) {
1063  fclose(diff);
1064  }
1065  } //end of particle loop
1066 
1067  if (all_out == 0) {
1068  printf("\n Number of particles that went out through out-flow boundary: %d \n", curr_n - 1);
1069  } else {
1070  printf("\n Number of particles completed %d, number of particles that went out through out-flow boundary: %d \n", curr_n - 1, curr_n - 1 - curr_o);
1071  }
1072 
1073  fclose(tp);
1074 
1075  if (tort_o > 0) {
1076  fclose(tort);
1077  }
1078 
1079  if (all_out > 0) {
1080  fclose(inp);
1081  }
1082 
1083  if (all_out > 0) {
1084  fclose(fnp);
1085  }
1086 
1087  if (frac_o > 0) {
1088  fclose(frac);
1089  }
1090 
1091  sprintf(filename, "%s/TotalNumberP", maindir);
1092  FILE *tn = OpenFile (filename, "w");
1093 
1094  if (all_out == 0) {
1095  fprintf(tn, " %10d \n", curr_n - 1);
1096  } else {
1097  fprintf(tn, " %10d %10d \n", curr_n - 1, curr_n - curr_o - 1);
1098  }
1099 
1100  fclose(tn);
1101 
1102  if (marfa == 1 || plumec == 1 ) {
1103  printf("\n Working on additional outputs \n");
1104  OutputMarPlumDisp (curr_n - 1, path);
1105  }
1106 
1107  return;
1108 }
1110 
1111 void CheckNewCell()
1113 {
1114  struct intcoef lambda;
1115  int n1 = 0, n2 = 0, n3 = 0, pcell;
1116  double delta_t;
1117  int cb = 0;
1118  n1 = cell[particle[np].cell - 1].node_ind[0];
1119  n2 = cell[particle[np].cell - 1].node_ind[1];
1120  n3 = cell[particle[np].cell - 1].node_ind[2];
1121  delta_t = CalculateCurrentDT();
1122  /**** first, calculate weights in the current cell****/
1123  lambda = CalculateWeights(n1, n2, n3);
1124 
1125  if ((lambda.weights[0] <= 1.) && (lambda.weights[0] >= 0.) && (lambda.weights[1] <= 1.) && (lambda.weights[1] >= 0.) && (lambda.weights[2] <= 1.) && (lambda.weights[2] >= 0.)) {
1126  /**** particle is in the current cell ***/
1127  particle[np].weight[0] = lambda.weights[0];
1128  particle[np].weight[1] = lambda.weights[1];
1129  particle[np].weight[2] = lambda.weights[2];
1130  } else {
1131  /* particle is not in the current cell ***/
1132  pcell = particle[np].cell;
1133  int pfract;
1134  pfract = particle[np].fracture;
1135  particle[np].cell = 0;
1136  SearchNeighborCells(n1, n2, n3);
1137  cb = 0;
1138 
1139  if ((node[n1 - 1].typeN == 210) || (node[n1 - 1].typeN == 212) || (node[n1 - 1].typeN == 200) || (node[n1 - 1].typeN == 202)) {
1140  cb = cb + 1;
1141  }
1142 
1143  if ((node[n2 - 1].typeN == 210) || (node[n2 - 1].typeN == 212) || (node[n2 - 1].typeN == 200) || (node[n2 - 1].typeN == 202)) {
1144  cb = cb + 1;
1145  }
1146 
1147  if ((node[n3 - 1].typeN == 210) || (node[n3 - 1].typeN == 212) || (node[n3 - 1].typeN == 200) || (node[n3 - 1].typeN == 202)) {
1148  cb = cb + 1;
1149  }
1150 
1151  if (cb != 0) {
1152  // printf(" Particle %d IS OUT of flow out zone. \n", np+1);
1153  FLAG_OUT = 1;
1154  }
1155 
1156  if (particle[np].cell == 0) {
1157  cb = 0;
1158 
1159  if (((node[n1 - 1].typeN == 10) || (node[n1 - 1].typeN == 12))) {
1160  cb = cb + 1;
1161  }
1162 
1163  if (((node[n2 - 1].typeN == 10) || (node[n2 - 1].typeN == 12))) {
1164  cb = cb + 1;
1165  }
1166 
1167  if (((node[n3 - 1].typeN == 10) || (node[n3 - 1].typeN == 12))) {
1168  cb = cb + 1;
1169  }
1170 
1171  if ((cb > 0) && (node[n1 - 1].typeN != 210) && (node[n1 - 1].typeN != 212) && (node[n2 - 1].typeN != 210) && (node[n2 - 1].typeN != 212) && (node[n3 - 1].typeN != 210) && (node[n3 - 1].typeN != 212)) {
1172  // printf("Particle %d is out of fracture %f %f %d %d\n", np+1, particle[np].position[0], particle[np].position[1], particle[np].fracture, pcell);
1173  /**** if out of fracture - make a flip in x direction ***/
1174  double cx, cy;
1175  cy = 1.0;
1176  cx = -1.0;
1177  particle[np].position[0] = particle[np].prev_pos[0] + delta_t*particle[np].velocity[0] * cx;
1178  particle[np].position[1] = particle[np].prev_pos[1] + delta_t*particle[np].velocity[1] * cy;
1179  SearchNeighborCells(n1, n2, n3);
1180 
1181  if (particle[np].cell == 0) {
1182  /**** if still out of fracture - make a flip in y direction ***/
1183  cy = -1.0;
1184  cx = 1.0;
1185  particle[np].position[0] = particle[np].prev_pos[0] + delta_t*particle[np].velocity[0] * cx;
1186  particle[np].position[1] = particle[np].prev_pos[1] + delta_t*particle[np].velocity[1] * cy;
1187  SearchNeighborCells(n1, n2, n3);
1188  }
1189 
1190  if (particle[np].cell == 0) {
1191  FLAG_OUT = 3;
1192  }
1193  }
1194 
1195  cb = 0;
1196 
1197  if ((node[n1 - 1].typeN == 210) || (node[n1 - 1].typeN == 212) || (node[n1 - 1].typeN == 200)) {
1198  cb = cb + 1;
1199  }
1200 
1201  if ((node[n2 - 1].typeN == 210) || (node[n2 - 1].typeN == 212) || (node[n2 - 1].typeN == 200)) {
1202  cb = cb + 1;
1203  }
1204 
1205  if ((node[n3 - 1].typeN == 210) || (node[n3 - 1].typeN == 212) || (node[n3 - 1].typeN == 200)) {
1206  cb = cb + 1;
1207  }
1208 
1209  if (cb != 0) {
1210  // printf(" Particle %d IS OUT of flow out zone. \n", np+1);
1211  FLAG_OUT = 1;
1212  } else {
1213  int ii = 0;
1214 
1215  for (ii = 0; ii < node[n1 - 1].numneighb; ii++) {
1216  if ((node[node[n1 - 1].indnodes[ii] - 1].typeN == 212) || (node[node[n1 - 1].indnodes[ii] - 1].typeN == 210)) {
1217  cb++;
1218  FLAG_OUT = 1;
1219  break;
1220  }
1221  }
1222 
1223  if (cb == 0) {
1224  for (ii = 0; ii < node[n2 - 1].numneighb; ii++) {
1225  if ((node[node[n2 - 1].indnodes[ii] - 1].typeN == 212) || (node[node[n2 - 1].indnodes[ii] - 1].typeN == 210)) {
1226  cb++;
1227  FLAG_OUT = 1;
1228  break;
1229  }
1230  }
1231  }
1232 
1233  if (cb == 0) {
1234  for (ii = 0; ii < node[n3 - 1].numneighb; ii++) {
1235  if ((node[node[n3 - 1].indnodes[ii] - 1].typeN == 212) || (node[node[n3 - 1].indnodes[ii] - 1].typeN == 210)) {
1236  cb++;
1237  FLAG_OUT = 1;
1238  break;
1239  }
1240  }
1241  }
1242  }
1243  }
1244  }
1245 
1246  return;
1247 }
1249 
1250 struct intcoef CalculateWeights(int nn1, int nn2, int nn3)
1252 {
1253  double n1x, n1y, n2x, n2y, n3x, n3y, deter;
1254  struct intcoef lambda;
1255  n1x = node[nn1 - 1].coord_xy[Xindex(nn1, np)];
1256  n1y = node[nn1 - 1].coord_xy[Yindex(nn1, np)];
1257  n2x = node[nn2 - 1].coord_xy[Xindex(nn2, np)];
1258  n2y = node[nn2 - 1].coord_xy[Yindex(nn2, np)];
1259  n3x = node[nn3 - 1].coord_xy[Xindex(nn3, np)];
1260  n3y = node[nn3 - 1].coord_xy[Yindex(nn3, np)];
1261  deter = (n2y - n3y) * (n1x - n3x) + (n3x - n2x) * (n1y - n3y);
1262  lambda.weights[0] = ((n2y - n3y) * (particle[np].position[0] - n3x) + (n3x - n2x) * (particle[np].position[1] - n3y)) / deter;
1263  lambda.weights[1] = ((n3y - n1y) * (particle[np].position[0] - n3x) + (n1x - n3x) * (particle[np].position[1] - n3y)) / deter;
1264  lambda.weights[2] = 1 - lambda.weights[0] - lambda.weights[1];
1265  double eps = 10e-5;
1266 
1267  if ((lambda.weights[0] >= -eps) && (lambda.weights[0] <= eps)) {
1268  lambda.weights[0] = 0.0;
1269  }
1270 
1271  if ((lambda.weights[1] >= -eps) && (lambda.weights[1] <= eps)) {
1272  lambda.weights[1] = 0.0;
1273  }
1274 
1275  lambda.weights[2] = 1 - lambda.weights[0] - lambda.weights[1];
1276 
1277  if ((lambda.weights[2] >= -eps) && (lambda.weights[2] <= eps)) {
1278  lambda.weights[2] = 0.0;
1279  }
1280 
1281  return lambda;
1282 }
1284 
1285 void SearchNeighborCells(int nn1, int nn2, int nn3)
1287 {
1288  int k = 0;
1289 
1290  /*search for cells of the vertex with highest intepolation weight - highest probability that particle is in one of it's neighbors*/
1291  if ((particle[np].weight[0] >= particle[np].weight[1]) && (particle[np].weight[0] >= particle[np].weight[2])) {
1292  k = nn1;
1293  }
1294 
1295  if ((particle[np].weight[1] >= particle[np].weight[0]) && (particle[np].weight[1] >= particle[np].weight[2])) {
1296  k = nn2;
1297  }
1298 
1299  if ((particle[np].weight[2] >= particle[np].weight[1]) && (particle[np].weight[2] >= particle[np].weight[0])) {
1300  k = nn3;
1301  }
1302 
1303  if (k != 0) {
1304  NeighborCells (k);
1305  }
1306 
1307  /* if not found, search between neighbors of other two nodes */
1308  if ((particle[np].cell == 0) && (k != nn1)) {
1309  NeighborCells (nn1);
1310  }
1311 
1312  if ((particle[np].cell == 0) && (k != nn2)) {
1313  NeighborCells (nn2);
1314  }
1315 
1316  if ((particle[np].cell == 0) && (k != nn3)) {
1317  NeighborCells (nn3);
1318  }
1319 
1320  return;
1321 }
1323 
1324 int InsideCell (unsigned int numc)
1327 {
1328  struct intcoef lambda;
1329  int nk_1, nk_2, nk_3, inside, nb = 0;
1330  nk_1 = cell[numc - 1].node_ind[0];
1331  nk_2 = cell[numc - 1].node_ind[1];
1332  nk_3 = cell[numc - 1].node_ind[2];
1333  double eps = 1e-5;
1334  lambda = CalculateWeights(nk_1, nk_2, nk_3);
1335  int intc = 0;
1336 
1337  if (particle[np].intcell == 4) {
1338  intc = 1;
1339  }
1340 
1341  if((lambda.weights[0] <= 1. + eps) && (lambda.weights[0] >= -eps) && (lambda.weights[1] <= 1. + eps) && (lambda.weights[1] >= -eps) && (lambda.weights[2] <= 1. + eps) && (lambda.weights[2] >= -eps)) {
1342  inside = 1;
1343  particle[np].weight[0] = fabs(lambda.weights[0]);
1344  particle[np].weight[1] = fabs(lambda.weights[1]);
1345  particle[np].weight[2] = fabs(lambda.weights[2]);
1346  particle[np].cell = numc;
1347 
1348  /* particle.intcell is a flag of particle being in intersection (=1) or boundary(=2) cell */
1349 
1350  if (node[nk_1 - 1].typeN == 10) {
1351  nb++;
1352  }
1353 
1354  if (node[nk_2 - 1].typeN == 10) {
1355  nb++;
1356  }
1357 
1358  if (node[nk_3 - 1].typeN == 10) {
1359  nb++;
1360  }
1361 
1362  if (nb > 1) {
1363  particle[np].intcell = 2;
1364  } else {
1365  particle[np].intcell = 0;
1366  }
1367 
1368  if (((node[nk_1 - 1].typeN == 2) || (node[nk_2 - 1].typeN == 2) || (node[nk_3 - 1].typeN == 2) || (node[nk_1 - 1].typeN == 12) || (node[nk_2 - 1].typeN == 12) || (node[nk_3 - 1].typeN == 12)) && (intc == 0)) {
1369  particle[np].intcell = 1;
1370  }
1371 
1372  if ((nb > 0) && (particle[np].intcell == 1)) {
1373  particle[np].intcell = 3;
1374  }
1375  } else {
1376  inside = 0;
1377  }
1378 
1379  return inside;
1380 }
1382 
1383 void PredictorStep()
1385 {
1386  int n1 = 0, n2 = 0, n3 = 0, v1 = 0, v2 = 0, v3 = 0;
1387  double delta_t;
1388  n1 = cell[particle[np].cell - 1].node_ind[0];
1389  n2 = cell[particle[np].cell - 1].node_ind[1];
1390  n3 = cell[particle[np].cell - 1].node_ind[2];
1391  v1 = cell[particle[np].cell - 1].veloc_ind[0];
1392  v2 = cell[particle[np].cell - 1].veloc_ind[1];
1393  v3 = cell[particle[np].cell - 1].veloc_ind[2];
1394  delta_t = CalculateCurrentDT();
1395  /*** velocity interpolation ***/
1396  particle[np].velocity[0] = particle[np].weight[0] * node[n1 - 1].velocity[v1][0] + particle[np].weight[1] * node[n2 - 1].velocity[v2][0] + particle[np].weight[2] * node[n3 - 1].velocity[v3][0];
1397  particle[np].velocity[1] = particle[np].weight[0] * node[n1 - 1].velocity[v1][1] + particle[np].weight[1] * node[n2 - 1].velocity[v2][1] + particle[np].weight[2] * node[n3 - 1].velocity[v3][1];
1398  particle[np].pressure = particle[np].weight[0] * node[n1 - 1].pressure + particle[np].weight[1] * node[n2 - 1].pressure + particle[np].weight[2] * node[n3 - 1].pressure;
1399  particle[np].prev_pos[0] = particle[np].position[0];
1400  particle[np].prev_pos[1] = particle[np].position[1];
1401  particle[np].position[0] = particle[np].position[0] + delta_t*particle[np].velocity[0];
1402  particle[np].position[1] = particle[np].position[1] + delta_t*particle[np].velocity[1];
1403  return;
1404 }
1405 
1407 
1408 void CorrectorStep()
1410 {
1411  int n1 = 0, n2 = 0, n3 = 0, v1 = 0, v2 = 0, v3 = 0;
1412  double delta_t;
1413  n1 = cell[particle[np].cell - 1].node_ind[0];
1414  n2 = cell[particle[np].cell - 1].node_ind[1];
1415  n3 = cell[particle[np].cell - 1].node_ind[2];
1416  v1 = cell[particle[np].cell - 1].veloc_ind[0];
1417  v2 = cell[particle[np].cell - 1].veloc_ind[1];
1418  v3 = cell[particle[np].cell - 1].veloc_ind[2];
1419  delta_t = CalculateCurrentDT();
1420  particle[np].position[0] = particle[np].prev_pos[0] + delta_t*particle[np].velocity[0];
1421  particle[np].position[1] = particle[np].prev_pos[1] + delta_t*particle[np].velocity[1];
1422  return;
1423 }
1425 
1426 void NeighborCells (int k)
1428 {
1429  int i = 0, j, inscell = 0;
1430  unsigned int nc;
1431 
1432  do {
1433  for (j = 0; j < 4; j++) {
1434  if ((node[k - 1].fracts[i][j] == particle[np].fracture)) {
1435  nc = node[k - 1].cells[i][j];
1436  inscell = InsideCell(nc);
1437  }
1438  }
1439 
1440  i++;
1441  } while ((inscell == 0) && (i < node[k - 1].numneighb));
1442 
1443  return;
1444 }
1446 
1447 void DefineTimeStep()
1449 {
1450  int i;
1451  double dotvel1, epsd = 1e-10;
1452 
1453  for (i = 0; i < nnodes; i++) {
1454  short int j = 0;
1455 
1456  for (j = 0; j < 4; j++) {
1457  dotvel1 = node[i].velocity[j][0] * node[i].velocity[j][0] + node[i].velocity[j][1] * node[i].velocity[j][1];
1458 
1459  if (dotvel1 > epsd) {
1460  if (node[i].typeN == 10) {
1461  node[i].timestep[j] = 0.005 * sqrt(((node[i].pvolume) / node[i].aperture) / dotvel1);
1462  } else {
1463  node[i].timestep[j] = 0.005 * sqrt(((node[i].pvolume) / node[i].aperture) / dotvel1);
1464  }
1465  } else {
1466  node[i].timestep[j] = 0.005 * sqrt(((node[i].pvolume) / node[i].aperture) / epsd);
1467  }
1468  }
1469  }
1470 
1471  return;
1472 }
1474 int CheckDistance()
1476 {
1477  /*** define the intersection segment/line ***/
1478  int dn1, dn2, dn3, int1 = 0, int2 = 0, int3 = 0, i, intm = 0;
1479  int ind_int2, fract_p;
1480  double px, py, dist, delta_t;
1481  double cx1 = 0, cy1 = 0, cx2 = 0, cy2 = 0;
1482  dn1 = cell[particle[np].cell - 1].node_ind[0];
1483  dn2 = cell[particle[np].cell - 1].node_ind[1];
1484  dn3 = cell[particle[np].cell - 1].node_ind[2];
1485  delta_t = CalculateCurrentDT();
1486  int prevcell = 0;
1487  prevcell = particle[np].cell;
1488  px = particle[np].position[0];
1489  py = particle[np].position[1];
1490  /* distance that particle will make during next step */
1491  dist = sqrt(pow((particle[np].velocity[0] * delta_t), 2) + pow((particle[np].velocity[1] * delta_t), 2));
1492 
1493  /* check: one edge of cell belongs to intersection line */
1494  if ((node[dn1 - 1].typeN == 2) || (node[dn1 - 1].typeN == 12)) {
1495  int1 = dn1;
1496  }
1497 
1498  if ((node[dn2 - 1].typeN == 2) || (node[dn2 - 1].typeN == 12)) {
1499  if (int1 != 0) {
1500  int2 = dn2;
1501  } else {
1502  int1 = dn2;
1503  }
1504  }
1505 
1506  if ((node[dn3 - 1].typeN == 2) || (node[dn3 - 1].typeN == 12)) {
1507  if (int1 != 0) {
1508  int2 = dn3;
1509  } else {
1510  int1 = dn3;
1511  }
1512  }
1513 
1514  if ((int1 != 0) && (int2 != 0)) {
1515  int3 = -1;
1516  }
1517 
1518  /* check: only one node of cell belongs to intersection - find the second node in neighboring list */
1519  if (int2 == 0) {
1520  for(i = 0; i < node[int1 - 1].numneighb; i++) {
1521  if (((node[int1 - 1].type[i] == 2) || (node[int1 - 1].type[i] == 12)) && ((node[node[int1 - 1].indnodes[i] - 1].fracture[0] == particle[np].fracture) || (node[node[int1 - 1].indnodes[i] - 1].fracture[1] == particle[np].fracture))) {
1522  if (int2 == 0) {
1523  int2 = node[int1 - 1].indnodes[i];
1524  } else {
1525  int3 = node[int1 - 1].indnodes[i];
1526  }
1527  }
1528 
1529  if ((int3 != 0) && (int3 != -1)) {
1530  break;
1531  }
1532  }
1533  }
1534 
1535  // printf("int1 %d int2 %d int3 %d cell %d frac %d \n", int1, int2, int3, particle[np].cell, particle[np].fracture);
1536 
1537  /*** if two nodes on intersection are found. the intersection edge will be defined *****/
1538  if ((int1 != 0) && (int2 != 0)) {
1539  cx1 = node[int1 - 1].coord_xy[Xindex(int1, np)];
1540  cy1 = node[int1 - 1].coord_xy[Yindex(int1, np)];
1541  cx2 = node[int2 - 1].coord_xy[Xindex(int2, np)];
1542  cy2 = node[int2 - 1].coord_xy[Yindex(int2, np)];
1543 
1544  if ((int1 != 0) && (int2 != 0) && (int3 != 0)) {
1545  /* define height of triangle, where the base is intersection segment */
1546  double base, side1, side2, perim, area, height;
1547  base = sqrt(((cx2 - cx1) * (cx2 - cx1)) + ((cy2 - cy1) * (cy2 - cy1)));
1548  side2 = sqrt(((cx2 - px) * (cx2 - px)) + ((cy2 - py) * (cy2 - py)));
1549  side1 = sqrt(((cx1 - px) * (cx1 - px)) + ((cy1 - py) * (cy1 - py)));
1550  perim = (side1 + side2 + base) / 2.0;
1551  area = sqrt(perim * (perim - base) * (perim - side1) * (perim - side2));
1552  height = area / (base * 0.5);
1553 
1554  /**** if distance between particle and intersection line is small enough, **/
1555  /**** particle makes an additional step and we check did it cross the intersection line or not**/
1556  /**** if particle crossed intersection, find the intersection point and move particle there**/
1557 
1558  if (height <= dist) {
1559  double pr1, pr2, pr3, pr4, px1, py1, px2, py2;
1560  px2 = particle[np].prev_pos[0];
1561  py2 = particle[np].prev_pos[1];
1562  PredictorStep();
1563  px1 = particle[np].position[0];
1564  py1 = particle[np].position[1];
1565  px1 = particle[np].position[0];
1566  py1 = particle[np].position[1];
1567  pr1 = (px1 - cx1) * (py2 - cy1) - (py1 - cy1) * (px2 - cx1);
1568  pr2 = (cx1 - px1) * (cy2 - py1) - (cy1 - py1) * (cx2 - px1);
1569  pr3 = (px1 - cx2) * (py2 - cy2) - (py1 - cy2) * (px2 - cx2);
1570  pr4 = (cx1 - px2) * (cy2 - py2) - (cy1 - py2) * (cx2 - px2);
1571 
1572  if ((pr1 * pr3 < 0) && (pr2 * pr4 < 0)) {
1573  /* particle crossed intersection */
1574  pr1 = cx1 * cy2 - cy1 * cx2;
1575  pr2 = px1 * py2 - py1 * px2;
1576  pr3 = (cx1 - cx2) * (py1 - py2) - (cy1 - cy2) * (px1 - px2);
1577  px = ((px1 - px2) * pr1 - (cx1 - cx2) * pr2) / pr3;
1578  py = ((py1 - py2) * pr1 - (cy1 - cy2) * pr2) / pr3;
1579  particle[np].position[0] = px;
1580  particle[np].position[1] = py;
1581  CheckNewCell();
1582  intm = 1;
1583 
1584  if (particle[np].cell != 0) {
1585  if (no_out != 1) {
1586  if (node[int1 - 1].fracture[0] != particle[np].fracture) {
1587  fract_p = node[int1 - 1].fracture[0];
1588  } else {
1589  fract_p = node[int1 - 1].fracture[1];
1590  }
1591 
1592  ParticleOutput(t, fract_p);
1593  }
1594 
1595  if (tdrw == 1) {
1596  t_adv = particle[np].time - t_adv0;
1597  timediff = TimeDomainRW(t_adv);
1598  t_adv0 = particle[np].time;
1599  particle[np].t_diff = particle[np].t_diff + timediff;
1600  particle[np].t_adv_diff = particle[np].t_adv_diff + t_adv + timediff;
1601 
1602  if (tdrw_o == 1) {
1603  fprintf(diff, "%5.12E %5.12E %5.12E %d %5.12E %5.12E %5.12E \n", t_adv, timediff, timediff + t_adv, particle[np].fracture, particle[np].time, particle[np].t_adv_diff, particle[np].t_diff);
1604  }
1605  }
1606 
1607  AcrossIntersection (prevcell, int1, int2, mixing_rule);
1608  }
1609  }
1610  }//if height
1611  }
1612 
1613  /* particle is in a cell at the ending point of intersection */
1614  if ((int1 != 0) && (int2 != 0) && (int3 == 0)) {
1615  double nposx = 0, nposy = 0, dist_init = 0, dist_fin = 0, inout = 0;
1616  int coutf = 0;
1617  nposx = (cx1 + cx2) / 2;
1618  nposy = (cy1 + cy2) / 2;
1619  dist_init = pow((particle[np].position[0] - nposx), 2) + pow((particle[np].position[1] - nposy), 2);
1620  PredictorStep();
1621  dist_fin = pow((particle[np].position[0] - nposx), 2) + pow((particle[np].position[1] - nposy), 2);
1622 
1623  if (dist_init >= dist_fin) {
1624  /* particle moves toward the intersection */
1625  particle[np].position[0] = nposx;
1626  particle[np].position[1] = nposy;
1627  CheckNewCell();
1628  intm = 1;
1629 
1630  if (particle[np].cell != 0) {
1631  prevcell = particle[np].cell;
1632 
1633  if (no_out != 1) {
1634  if (node[int1 - 1].fracture[0] != particle[np].fracture) {
1635  fract_p = node[int1 - 1].fracture[0];
1636  } else {
1637  fract_p = node[int1 - 1].fracture[1];
1638  }
1639 
1640  ParticleOutput(t, fract_p);
1641  }
1642 
1643  if (tdrw == 1) {
1644  t_adv = particle[np].time - t_adv0;
1645  timediff = TimeDomainRW(t_adv);
1646  t_adv0 = particle[np].time;
1647  particle[np].t_diff = particle[np].t_diff + timediff;
1648  particle[np].t_adv_diff = particle[np].t_adv_diff + t_adv + timediff;
1649 
1650  if (tdrw_o == 1) {
1651  fprintf(diff, "%5.12E %5.12E %5.12E %d %5.12E %5.12E %5.12E \n", t_adv, timediff, timediff + t_adv, particle[np].fracture, particle[np].time, particle[np].t_adv_diff, particle[np].t_diff);
1652  }
1653  }
1654 
1655  AcrossIntersection (prevcell, int1, int2, mixing_rule);
1656  } else {
1657  printf("Particle is lost on end of intersection. \n");
1658  }
1659  } else {
1660  ind_int2 = 0;
1661 
1662  for (i = 0; i < node[int1 - 1].numneighb; i++) {
1663  if (node[int1 - 1].indnodes[i] == int2) {
1664  ind_int2 = i;
1665  break;
1666  }
1667  }
1668 
1669  if (ind_int2 == 0)
1670 
1671  // printf("ind=0 \n");
1672  for (i = 0; i < 4; i++) {
1673  if (node[int1 - 1].fracts[ind_int2][i] == particle[np].fracture) {
1674  int indc = 0;
1675  indc = node[int1 - 1].cells[ind_int2][i];
1676  inout = InOutFlowCell(indc, int1, nposx, nposy);
1677 
1678  if (inout > 0) {
1679  coutf++;
1680  }
1681  }
1682  }
1683 
1684  if (coutf > 1) {
1685  particle[np].position[0] = nposx;
1686  particle[np].position[1] = nposy;
1687  CheckNewCell();
1688  intm = 1;
1689 
1690  if (particle[np].cell != 0) {
1691  prevcell = particle[np].cell;
1692 
1693  if (no_out != 1) {
1694  if (node[int1 - 1].fracture[0] != particle[np].fracture) {
1695  fract_p = node[int1 - 1].fracture[0];
1696  } else {
1697  fract_p = node[int1 - 1].fracture[1];
1698  }
1699 
1700  ParticleOutput(t, fract_p);
1701  }
1702 
1703  if (tdrw == 1) {
1704  t_adv = particle[np].time - t_adv0;
1705  timediff = TimeDomainRW(t_adv);
1706  t_adv0 = particle[np].time;
1707  particle[np].t_diff = particle[np].t_diff + timediff;
1708  particle[np].t_adv_diff = particle[np].t_adv_diff + t_adv + timediff;
1709 
1710  if (tdrw_o == 1) {
1711  fprintf(diff, "%5.12E %5.12E %5.12E %d %5.12E %5.12E %5.12E \n", t_adv, timediff, timediff + t_adv, particle[np].fracture, particle[np].time, particle[np].t_adv_diff, particle[np].t_diff);
1712  }
1713  }
1714 
1715  AcrossIntersection (prevcell, int1, int2, mixing_rule);
1716  }
1717  }
1718  }
1719  }
1720  }
1721 
1722  if ((int1 == 0) || (int2 == 0)) {
1723  /* there is no node of cell belongs to intersection */
1724  // printf("Check if %d cell is on intersection %d %d %d %d %d dn %d %d %d %d. \n", particle[np].cell, node[dn1-1].typeN, node[dn2-1].typeN, node[dn3-1].typeN, int1, int2, dn1, dn2, dn3, particle[np].intcell);
1725  }
1726 
1727  return intm;
1728 }
1730 double InOutFlowCell(int indcell, int int1, double nposx, double nposy)
1732 {
1733  double inoutf = 0;
1734  int n1n, n2n, n3n, v1v, v2v, v3v;
1735  struct intcoef lambda;
1736  double prevpos0 = particle[np].position[0], prevpos1 = particle[np].position[1];
1737  int prevfract = particle[np].fracture, previouscell = particle[np].cell;
1738  double products = 0, product = 0;
1739  particle[np].position[0] = nposx;
1740  particle[np].position[1] = nposy;
1741  n1n = cell[indcell - 1].node_ind[0];
1742  n2n = cell[indcell - 1].node_ind[1];
1743  n3n = cell[indcell - 1].node_ind[2];
1744  v1v = cell[indcell - 1].veloc_ind[0];
1745  v2v = cell[indcell - 1].veloc_ind[1];
1746  v3v = cell[indcell - 1].veloc_ind[2];
1747  int thirdnode = 0;
1748  double tnx = 0, tny = 0;
1749 
1750  if ((node[n1n - 1].typeN != 12) && (node[n1n - 1].typeN != 2)) {
1751  thirdnode = n1n;
1752  }
1753 
1754  if ((node[n2n - 1].typeN != 12) && (node[n2n - 1].typeN != 2)) {
1755  thirdnode = n2n;
1756  }
1757 
1758  if ((node[n3n - 1].typeN != 12) && (node[n3n - 1].typeN != 2)) {
1759  thirdnode = n3n;
1760  }
1761 
1762  tnx = node[thirdnode - 1].coord_xy[0];
1763  tny = node[thirdnode - 1].coord_xy[1];
1764  double vintx = 0, vinty = 0, velocx = 0, velocy = 0;
1765  vintx = node[int1 - 1].coord_xy[XindexC(int1, indcell - 1)];
1766  vinty = node[int1 - 1].coord_xy[YindexC(int1, indcell - 1)];
1767  ChangeFracture(indcell);
1768  lambda = CalculateWeights(n1n, n2n, n3n);
1769  velocx = lambda.weights[0] * node[n1n - 1].velocity[v1v][0] + lambda.weights[1] * node[n2n - 1].velocity[v2v][0] + lambda.weights[2] * node[n3n - 1].velocity[v3v][0];
1770  velocy = lambda.weights[0] * node[n1n - 1].velocity[v1v][1] + lambda.weights[1] * node[n2n - 1].velocity[v2v][1] + lambda.weights[2] * node[n3n - 1].velocity[v3v][1];
1771  /* calculate vector's cross product to define outgoing and incoming flow cells */
1772  product = ((particle[np].position[0] - tnx) * (particle[np].position[1] - vinty)) - ((particle[np].position[0] - vintx) * (particle[np].position[1] - tny));
1773  products = (velocx * (particle[np].position[1] - vinty)) - ((particle[np].position[0] - vintx) * velocy);
1774  inoutf = products * product;
1775  particle[np].position[0] = prevpos0;
1776  particle[np].position[1] = prevpos1;
1777  particle[np].cell = previouscell;
1778  particle[np].fracture = prevfract;
1779  return inoutf;
1780 }
1781 
1783 
1784 void AcrossIntersection (int prevcell, int int1, int int2, int mixing_rule)
1786 {
1787  struct intcoef lambda;
1788  int indj = -1, k, indcell, cell_win = 0, indk = 0;
1789  int n1n, n2n, n3n, v1v, v2v, v3v;
1790  int prevfrac; // the fracture of previous cell;
1791  int finalfrac; // the fracture the particle moved at the intersection
1792  double speedsq[4] = {0.0, 0.0, 0.0, 0.0}, velocx, velocy;
1793  double products[4] = {0.0, 0.0, 0.0, 0.0};
1794  int neighborcellind[4]; //a list of cell indices of all the neighbors
1795  int neighborfracind[4]; // a list of fractures of all the neighbors
1796  //int rule =1; // if 0 streamline, 1 complete mixing
1797  prevfrac = cell[prevcell - 1].fracture;
1798 
1799  if ((int1 != 0) && (int2 != 0)) {
1800  /* defines indj - index of current cell in int1 node list */
1801  k = 0;
1802 
1803  do {
1804  if (node[int1 - 1].indnodes[k] == int2) {
1805  indj = k;
1806  }
1807 
1808  k++;
1809  } while ((indj < 0) && (k < node[int1 - 1].numneighb));
1810 
1811  if (indj < 0) {
1812  printf(" Current cell not found: NODES %d %d pw %f %f %f \n", int1, int2, particle[np].weight[0], particle[np].weight[1], particle[np].weight[2]);
1813  }
1814 
1815  /* the loop on 4 neighboring cells with common edge: int1 - int2 */
1816  for (k = 0; k < 4; k++) {
1817  if (node[int1 - 1].cells[indj][k] != 0) {
1818  indcell = node[int1 - 1].cells[indj][k];
1819  neighborcellind[k] = indcell;
1820  n1n = cell[indcell - 1].node_ind[0];
1821  n2n = cell[indcell - 1].node_ind[1];
1822  n3n = cell[indcell - 1].node_ind[2];
1823  v1v = cell[indcell - 1].veloc_ind[0];
1824  v2v = cell[indcell - 1].veloc_ind[1];
1825  v3v = cell[indcell - 1].veloc_ind[2];
1826  int thirdnode = 0;
1827  double tnx = 0, tny = 0;
1828 
1829  if ((node[n1n - 1].typeN != 12) && (node[n1n - 1].typeN != 2)) {
1830  thirdnode = n1n;
1831  }
1832 
1833  if ((node[n2n - 1].typeN != 12) && (node[n2n - 1].typeN != 2)) {
1834  thirdnode = n2n;
1835  }
1836 
1837  if ((node[n3n - 1].typeN != 12) && (node[n3n - 1].typeN != 2)) {
1838  thirdnode = n3n;
1839  }
1840 
1841  tnx = node[thirdnode - 1].coord_xy[0];
1842  tny = node[thirdnode - 1].coord_xy[1];
1843  neighborfracind[k] = node[thirdnode - 1].fracture[0]; // find the fracture of third node:
1844  double product, vintx = 0, vinty = 0;
1845  vintx = node[int1 - 1].coord_xy[XindexC(int1, indcell - 1)];
1846  vinty = node[int1 - 1].coord_xy[YindexC(int1, indcell - 1)];
1847 
1848  if ((indcell) == prevcell) {
1849  indk = k;
1850  //printf("Found index cell \n");
1851  }
1852 
1853  /**** move to the intersecting fracture and recalculate coordinations ***/
1854  ChangeFracture(indcell);
1855  lambda = CalculateWeights(n1n, n2n, n3n);
1856  velocx = lambda.weights[0] * node[n1n - 1].velocity[v1v][0] + lambda.weights[1] * node[n2n - 1].velocity[v2v][0] + lambda.weights[2] * node[n3n - 1].velocity[v3v][0];
1857  velocy = lambda.weights[0] * node[n1n - 1].velocity[v1v][1] + lambda.weights[1] * node[n2n - 1].velocity[v2v][1] + lambda.weights[2] * node[n3n - 1].velocity[v3v][1];
1858  /* calculate vector's cross product to define outgoing and incoming flow cells */
1859  product = ((particle[np].position[0] - tnx) * (particle[np].position[1] - vinty)) - ((particle[np].position[0] - vintx) * (particle[np].position[1] - tny));
1860  products[k] = (velocx * (particle[np].position[1] - vinty)) - ((particle[np].position[0] - vintx) * velocy);
1861  products[k] = products[k] * product;
1862  speedsq[k] = velocx * velocx + velocy * velocy;
1863  }
1864  }//loop on k
1865 
1866  //printf("Speed of each cell: %lf, %lf, %lf, %lf, %lf \n", sqrt(speedsq[0]), sqrt(speedsq[1]), sqrt(speedsq[2]),(speedsq[3]), sqrt(4));
1867  if (mixing_rule == 1) {
1868  //printf("Complete Mixing \n");
1869  cell_win = CompleteMixingRandomSampling(products, speedsq, indj, int1, indk);
1870  }
1871 
1872  if(mixing_rule == 0) {
1873  //printf("Streamline Routing \n");
1874  cell_win = StreamlineRandomSampling(products, speedsq, indj, int1, indk, neighborcellind, neighborfracind, prevfrac, prevcell);
1875  }
1876 
1877  ChangeFracture(cell_win);
1878  particle[np].intcell = 4;
1879  particle[np].prev_pos[0] = particle[np].position[0];
1880  particle[np].prev_pos[1] = particle[np].position[1];
1881  finalfrac = cell[cell_win - 1].fracture;
1882  }
1883 
1884  return;
1885 }
1887 /**************************************************************************/
1888 /******************** Streamline Routing Intersection Rule ***************************/
1901 int StreamlineRandomSampling(double products[4], double speedsq[4], int indj, int int1, int indk, int neighborcellind[4], int neighborfracind[4], int prevfrac, int prevcell) {
1902  /*********** Streamline Routing Sampling ****************/
1903  int win_cell = 0, k, jj;
1904  int count = 0, outc[4] = {0, 0, 0, 0};
1905  int oppcellind; // the index of the opposite cell
1906  int oppflag = 1; // opp flag set to 1 is we have outflow on the opposite cell
1907  double random_number, totalspeed = 0, minsp = 0.0, incomingspeedsq = 0, speedopp = 0, totalmag = 0;
1908 
1909  /* find outgoing flow cells */
1910  for (k = 0; k < 4; k++) {
1911  if (products[k] < 0) {
1912  outc[count] = k;
1913  count++;
1914  totalspeed = totalspeed + speedsq[k]; //the sum of the speed square of outflowing cells
1915  totalmag = totalmag + sqrt(speedsq[k]); // the sum of speeds of outflowing cells
1916  }
1917  }
1918 
1919  /* if no outgoing flow cells found - move to cell with largest velocity magnitude
1920  (this should not happen, it will mean we have a physical flow sink)*/
1921  if (count == 0) {
1922  // printf("Case 0 \n");
1923  minsp = 0.0;
1924 
1925  for (k = 0; k < 4; k++)
1926  if ((speedsq[k] > minsp) && (k != indk)) {
1927  minsp = speedsq[k];
1928  win_cell = node[int1 - 1].cells[indj][k];
1929  }
1930  }
1931 
1932  /* if only one cell found */
1933  if (count == 1) {
1934  win_cell = node[int1 - 1].cells[indj][outc[0]];
1935  }
1936 
1937  if (count == 2) {
1938  // loop to find opposite cell. If we find an outflowing opposite cell(cell with same fracture as previous cell), then we are in continuous case and oppflag is set to 1;
1939  for (jj = 0; jj < 4; jj++) {
1940  if ((neighborfracind[jj] == prevfrac) && (neighborcellind[jj] == node[int1 - 1].cells[indj][outc[0]] || neighborcellind[jj] == node[int1 - 1].cells[indj][outc[1]])) { //Opposite flag turned on if: the neighboring cell has the same fracture as the previous cell and that cell is one of the two outflowing cells
1941  oppcellind = neighborcellind[jj];
1942  oppflag = 1; //we are in the continous case
1943  }
1944 
1945  //note this next case only will work if we have continuous case.
1946  if ((neighborfracind[jj] == prevfrac) && neighborcellind[jj] != node[int1 - 1].cells[indj][outc[0]] && neighborcellind[jj] != node[int1 - 1].cells[indj][outc[1]]) { // We find the incoming speed by finding the neighbor cell that shares a fracture with previous cell, and is not outflowing. Note we code this way because in special cases the previous cell is not in the list of cells neighboring the intersection. In this case, we say the incoming velocity is the inflowing cell and has the same fracture of the previous cell.
1947  incomingspeedsq = speedsq[jj];
1948  }
1949 
1950  if(neighborcellind[jj] == prevcell) { //This is the easier way to find incoming speed: We find the neighboring cell equal to the previous cell and take the speed of that cell.
1951  incomingspeedsq = speedsq[jj];
1952  }
1953  } // end of jj loop
1954 
1955  if (cell[node[int1 - 1].cells[indj][outc[1]]].fracture == cell[node[int1 - 1].cells[indj][outc[0]]].fracture) { // if the two outflowing cells have the same fracture we are in the discontinuous case
1956  oppflag = 0;
1957  }
1958 
1959  random_number = drand48();
1960 
1961  if (oppflag == 1) { // if oppflag is 1 then one of the outflow cells is opposite (on the same fracture) as the cell we are coming from: The other outflowing cell must be adjacent
1962  // printf("We are in the continous case \n");
1963  // We need to find which outflowing cell is opposite
1964  if (node[int1 - 1].cells[indj][outc[0]] == oppcellind) { // opposite cell has index outc[0] and outc[1] is adjacent
1965  if (incomingspeedsq <= speedsq[outc[1]]) { // The adjacent cell has large velocity and so we are forced to the adjacent branch
1966  win_cell = node[int1 - 1].cells[indj][outc[1]];
1967  } else { // the adjacent velocity is not bigger the incoming
1968  if (random_number <= sqrt(speedsq[outc[1]]) / sqrt(incomingspeedsq)) {
1969  win_cell = node[int1 - 1].cells[indj][outc[1]];
1970  } else {
1971  win_cell = node[int1 - 1].cells[indj][outc[0]];
1972  }
1973  }
1974  } // end first if regarding opposite cell has index outc[0]
1975  else { // other Continous case scenario: opposite has index outc[1] and adjacent is outc[0]
1976  if (incomingspeedsq <= speedsq[outc[0]]) { //The adjacent cell has large velocity than incoming
1977  win_cell = node[int1 - 1].cells[indj][outc[0]];
1978  } else { // the adjacent (outc[0]) velocity is smaller than incoming
1979  if (random_number <= sqrt(speedsq[outc[0]]) / sqrt(incomingspeedsq)) {
1980  win_cell = node[int1 - 1].cells[indj][outc[0]];
1981  } else {
1982  win_cell = node[int1 - 1].cells[indj][outc[1]];
1983  }
1984  }
1985  }//end case where outc[1] is opposite
1986  } // end continous case
1987  else { // start discontinous case now: assume complete mixing for now
1988  //printf("We are in the discontinous case \n");
1989  //printf("Speed incoming = %lf: Other Speed incoming = %lf, Outflowing1 =%lf: Outflowing2 =%lf: \n", sqrt(incomingspeedsq), sqrt(speedsq[outc[0]]) +sqrt(speedsq[outc[1]]) - sqrt(incomingspeedsq),sqrt(speedsq[outc[0]]),sqrt(speedsq[outc[1]]) );
1990  //New Rule
1991  if(sqrt(incomingspeedsq) >= (sqrt(speedsq[outc[0]]) + sqrt(speedsq[outc[1]])) / 2) { //the incoming speed is bigger than oppposite incoming: we are on the larger incoming fracture
1992  if(speedsq[outc[1]] >= speedsq[outc[0]]) { //outgoing cell 1 is bigger than outgoint cell 0
1993  //printf("On high; Prob of switching to high frac: %lf \n", sqrt(speedsq[outc[1]])/sqrt(incomingspeedsq));
1994  if(random_number >= sqrt(speedsq[outc[1]]) / sqrt(incomingspeedsq)) {
1995  win_cell = node[int1 - 1].cells[indj][outc[0]];
1996  } else {
1997  win_cell = node[int1 - 1].cells[indj][outc[1]];
1998  }
1999  } else { //outgoing cell 0 is bigger
2000  //printf("On high; Prob of switching to high frac: %lf \n", sqrt(speedsq[outc[0]])/sqrt(incomingspeedsq));
2001  if(random_number >= sqrt(speedsq[outc[0]]) / sqrt(incomingspeedsq)) {
2002  win_cell = node[int1 - 1].cells[indj][outc[1]];
2003  } else {
2004  win_cell = node[int1 - 1].cells[indj][outc[0]];
2005  }
2006  }
2007  } else { // the incoming speed is smaller than oppoiste incoming
2008  speedopp = sqrt(speedsq[outc[0]]) + sqrt(speedsq[outc[1]]) - sqrt(incomingspeedsq); // the speed of the larger incoming branch
2009 
2010  if(speedsq[outc[1]] <= speedsq[outc[0]]) { //outgoing cell 1 is smaller than outgoing cell 0
2011  //printf("On low; Prob of switching to high frac: %lf \n", (sqrt(speedsq[outc[0]])-speedopp)/sqrt(incomingspeedsq));
2012  if(random_number <= (sqrt(speedsq[outc[0]]) - speedopp) / sqrt(incomingspeedsq) && sqrt(speedsq[outc[0]]) >= speedopp) {
2013  win_cell = node[int1 - 1].cells[indj][outc[0]];
2014  } else {
2015  win_cell = node[int1 - 1].cells[indj][outc[1]];
2016  }
2017  } else { //outgoing cell 1 is bigger
2018  // printf("On low; Prob of switching to high frac: %lf \n", (sqrt(speedsq[outc[1]])-speedopp)/sqrt(incomingspeedsq));
2019  if(random_number <= (sqrt(speedsq[outc[1]]) - speedopp) / sqrt(incomingspeedsq) && sqrt(speedsq[outc[1]]) >= speedopp) {
2020  win_cell = node[int1 - 1].cells[indj][outc[1]];
2021  } else {
2022  win_cell = node[int1 - 1].cells[indj][outc[0]];
2023  }
2024  }
2025  }
2026 
2027  // End New Rule
2028  //Start Old Rule: Complete mixing
2029  /****if (random_number<=(sqrt(speedsq[outc[0]])/totalmag))
2030  win_cell=node[int1-1].cells[indj][outc[0]];
2031  else
2032  win_cell=node[int1-1].cells[indj][outc[1]]; ****/
2033  // End old rule
2034  } // end discontinous case
2035  }
2036 
2037  if (count == 3) {
2038  random_number = drand48();
2039 
2040  if (random_number > ((sqrt(speedsq[outc[0]]) + sqrt(speedsq[outc[1]])) / totalmag))
2041  // if (random_number<0.3)
2042  {
2043  win_cell = node[int1 - 1].cells[indj][outc[2]];
2044  } else {
2045  if (random_number <= (sqrt(speedsq[outc[0]]) / totalmag))
2046  // if (random_number<0.6)
2047  {
2048  win_cell = node[int1 - 1].cells[indj][outc[0]];
2049  } else {
2050  win_cell = node[int1 - 1].cells[indj][outc[1]];
2051  }
2052  }
2053  }
2054 
2055  if (count == 4) {
2056  random_number = drand48();
2057 
2058  if (random_number > ((speedsq[outc[0]] + speedsq[outc[1]] + speedsq[outc[2]]) / totalspeed)) {
2059  win_cell = node[int1 - 1].cells[indj][outc[3]];
2060  } else {
2061  if (random_number > ((speedsq[outc[0]] + speedsq[outc[1]]) / totalspeed)) {
2062  win_cell = node[int1 - 1].cells[indj][outc[2]];
2063  } else {
2064  if (random_number <= (speedsq[outc[0]] / totalspeed)) {
2065  win_cell = node[int1 - 1].cells[indj][outc[0]];
2066  } else {
2067  win_cell = node[int1 - 1].cells[indj][outc[1]];
2068  }
2069  }
2070  }
2071  }
2072 
2073  return win_cell;
2074 }
2075 
2077 /**************************************************************************/
2078 /******************** Complete Mixing Intersection Rule ***************************/
2089 int CompleteMixingRandomSampling(double products[4], double speedsq[4], int indj, int int1, int indk) {
2090  /***********Complete Mixing Sampling ****************/
2091  int win_cell = 0, k;
2092  int count = 0, outc[4] = {0, 0, 0, 0};
2093  double random_number, totalspeed = 0, minsp = 0.0, totalmag = 0;
2094 
2095  /* find outgoing flow cells */
2096  for (k = 0; k < 4; k++) {
2097  if (products[k] < 0) {
2098  outc[count] = k;
2099  count++;
2100  totalspeed = totalspeed + speedsq[k];
2101  totalmag = totalmag + sqrt(speedsq[k]);
2102  }
2103  }
2104 
2105  /* if no outgoing flow cells found - move to cell with largest velocity magnitude
2106  (this should not happen, it will mean we have a physical flow sink)*/
2107  if (count == 0) {
2108  minsp = 0.0;
2109 
2110  for (k = 0; k < 4; k++)
2111  if ((speedsq[k] > minsp) && (k != indk)) {
2112  minsp = speedsq[k];
2113  win_cell = node[int1 - 1].cells[indj][k];
2114  }
2115  }
2116 
2117  /* if only one cell found */
2118  if (count == 1) {
2119  win_cell = node[int1 - 1].cells[indj][outc[0]];
2120  }
2121 
2122  if (count == 2) {
2123  //printf("Case 2 \n");
2124  random_number = drand48();
2125 
2126  if (random_number <= (sqrt(speedsq[outc[0]]) / totalmag))
2127  // if (random_number<0.5)
2128  {
2129  win_cell = node[int1 - 1].cells[indj][outc[0]];
2130  } else {
2131  win_cell = node[int1 - 1].cells[indj][outc[1]];
2132  }
2133  }
2134 
2135  if (count == 3) {
2136  random_number = drand48();
2137 
2138  if (random_number > ((sqrt(speedsq[outc[0]]) + sqrt(speedsq[outc[1]])) / totalmag)) {
2139  // if (random_number<0.3)
2140  win_cell = node[int1 - 1].cells[indj][outc[2]];
2141  //printf("winning cell = %d \n", win_cell);
2142  } else {
2143  if (random_number <= (sqrt(speedsq[outc[0]]) / totalmag))
2144  // if (random_number<0.6)
2145  {
2146  win_cell = node[int1 - 1].cells[indj][outc[0]];
2147  } else {
2148  win_cell = node[int1 - 1].cells[indj][outc[1]];
2149  }
2150  }
2151  }
2152 
2153  if (count == 4) {
2154  random_number = drand48();
2155 
2156  if (random_number > ((speedsq[outc[0]] + speedsq[outc[1]] + speedsq[outc[2]]) / totalspeed)) {
2157  win_cell = node[int1 - 1].cells[indj][outc[3]];
2158  } else {
2159  if (random_number > ((speedsq[outc[0]] + speedsq[outc[1]]) / totalspeed)) {
2160  win_cell = node[int1 - 1].cells[indj][outc[2]];
2161  } else {
2162  if (random_number <= (speedsq[outc[0]] / totalspeed)) {
2163  win_cell = node[int1 - 1].cells[indj][outc[0]];
2164  } else {
2165  win_cell = node[int1 - 1].cells[indj][outc[1]];
2166  }
2167  }
2168  }
2169  }
2170 
2171  return win_cell;
2172 }
2173 
2175 void Moving2Center (int nnp, int cellnumber)
2177 {
2178  double centx = 0, centy = 0, n1x = 0, n2x = 0, n3x = 0, n1y = 0, n2y = 0, n3y = 0;
2179  int n1 = 0, n2 = 0, n3 = 0;
2180  n1 = cell[cellnumber - 1].node_ind[0];
2181  n2 = cell[cellnumber - 1].node_ind[1];
2182  n3 = cell[cellnumber - 1].node_ind[2];
2183  n1x = node[n1 - 1].coord_xy[Xindex(n1, nnp)];
2184  n1y = node[n1 - 1].coord_xy[Yindex(n1, nnp)];
2185  n2x = node[n2 - 1].coord_xy[Xindex(n2, nnp)];
2186  n2y = node[n2 - 1].coord_xy[Yindex(n2, nnp)];
2187  n3x = node[n3 - 1].coord_xy[Xindex(n3, nnp)];
2188  n3y = node[n3 - 1].coord_xy[Yindex(n3, nnp)];
2189  centx = n1x + n2x + n3x;
2190  centy = n1y + n2y + n3y;
2191  particle[nnp].position[0] = centx / 3;
2192  particle[nnp].position[1] = centy / 3;
2193  int in;
2194  in = InsideCell (particle[nnp].cell);
2195  return;
2196 }
2198 int Moving2NextCell (int stuck, int k)
2200 {
2201  int current_index = 0;
2202  int j, nc = 0, i = 0;
2203 
2204  if (stuck < node[k - 1].numneighb) {
2205  i = stuck;
2206 
2207  do {
2208  for (j = 0; j < 4; j++) {
2209  if (node[k - 1].fracts[i][j] == particle[np].fracture) {
2210  if (node[k - 1].cells[i][j] != particle[np].cell) {
2211  nc = node[k - 1].cells[i][j];
2212  current_index = i + 1;
2213  // printf("current_index %d i %d j %d \n", current_index, i,j);
2214  break;
2215  }
2216  }
2217  }
2218 
2219  if (nc != 0) {
2220  // printf("particle %d moved from cell %d to cell %d \n", np+1, particle[np].cell, nc);
2221  particle[np].cell = nc;
2222  Moving2Center (np, nc);
2223  break;
2224  }
2225 
2226  i++;
2227  } while ((nc == 0) && (i < node[k - 1].numneighb));
2228  } else {
2229  // printf(" 'stuck' is too big \n");
2230  }
2231 
2232  if (nc == 0) {
2233  // printf("moving cell was not found part %d fract %d cell %d k %d\n", np+1, particle[np].fracture, particle[np].cell, k);
2234  }
2235 
2236  return current_index;
2237 }
2239 double CalculateCurrentDT()
2241 {
2242  int dn1 = 0, dn2 = 0, dn3 = 0, dv1 = 0, dv2 = 0, dv3 = 0;
2243  double current_delta_t;
2244  dn1 = cell[particle[np].cell - 1].node_ind[0];
2245  dn2 = cell[particle[np].cell - 1].node_ind[1];
2246  dn3 = cell[particle[np].cell - 1].node_ind[2];
2247  dv1 = cell[particle[np].cell - 1].veloc_ind[0];
2248  dv2 = cell[particle[np].cell - 1].veloc_ind[1];
2249  dv3 = cell[particle[np].cell - 1].veloc_ind[2];
2250  current_delta_t = node[dn1 - 1].timestep[dv1] * particle[np].weight[0] + node[dn2 - 1].timestep[dv2] * particle[np].weight[1] + node[dn3 - 1].timestep[dv3] * particle[np].weight[2];
2251  return current_delta_t;
2252 }
2254 int Xindex(int nodenum, int nnp)
2256 {
2257  int xind = 0;
2258 
2259  if (node[nodenum - 1].fracture[0] == particle[nnp].fracture) {
2260  xind = 0;
2261  }
2262 
2263  if (node[nodenum - 1].fracture[1] == particle[nnp].fracture) {
2264  xind = 3;
2265  }
2266 
2267  return xind;
2268 }
2270 int Yindex(int nodenum, int nnp)
2272 {
2273  int yind = 0;
2274 
2275  if (node[nodenum - 1].fracture[0] == particle[nnp].fracture) {
2276  yind = 1;
2277  }
2278 
2279  if (node[nodenum - 1].fracture[1] == particle[nnp].fracture) {
2280  yind = 4;
2281  }
2282 
2283  return yind;
2284 }
2286 void Moving2NextCellBound(int prevcell)
2288 {
2289  int n1 = 0, n2 = 0, n3 = 0;
2290  n1 = cell[prevcell - 1].node_ind[0];
2291  n2 = cell[prevcell - 1].node_ind[1];
2292  n3 = cell[prevcell - 1].node_ind[2];
2293  int k, j, nc = 0, i = 0, nc0 = 0, nc10 = 0, nc12 = 0, ncb = 0;
2294  k = 0;
2295 
2296  if (node[n1 - 1].typeN == 0) {
2297  k = n1;
2298  } else {
2299  if (node[n2 - 1].typeN == 0) {
2300  k = n2;
2301  } else {
2302  k = n3;
2303  }
2304  }
2305 
2306  do {
2307  for (j = 0; j < 4; j++) {
2308  if (node[k - 1].fracts[i][j] == particle[np].fracture) {
2309  if (node[k - 1].cells[i][j] != prevcell) {
2310  nc = node[k - 1].cells[i][j];
2311 
2312  if (node[cell[nc - 1].node_ind[0] - 1].typeN + node[cell[nc - 1].node_ind[1] - 1].typeN + node[cell[nc - 1].node_ind[2] - 1].typeN == 0) {
2313  nc0 = nc;
2314  break;
2315  }
2316 
2317  if (node[cell[nc - 1].node_ind[0] - 1].typeN + node[cell[nc - 1].node_ind[1] - 1].typeN + node[cell[nc - 1].node_ind[2] - 1].typeN == 10) {
2318  nc10 = nc;
2319  }
2320 
2321  if (node[cell[nc - 1].node_ind[0] - 1].typeN + node[cell[nc - 1].node_ind[1] - 1].typeN + node[cell[nc - 1].node_ind[2] - 1].typeN == 12) {
2322  nc12 = nc;
2323  }
2324 
2325  if (node[cell[nc - 1].node_ind[0] - 1].typeN + node[cell[nc - 1].node_ind[1] - 1].typeN + node[cell[nc - 1].node_ind[2] - 1].typeN > 15) {
2326  ncb = nc;
2327  }
2328  }
2329  }
2330  }
2331 
2332  if (nc0 != 0) {
2333  // printf("particle %d moved from boundary cell %d to cell %d \n", np+1, prevcell, nc0);
2334  particle[np].cell = nc0;
2335  Moving2Center (np, nc0);
2336  break;
2337  }
2338 
2339  i++;
2340  } while (i < node[k - 1].numneighb);
2341 
2342  if ((nc0 == 0) && (nc10 != 0)) {
2343  // printf("particle %d moved from boundary cell %d to cell %d \n", np+1, prevcell, nc10);
2344  particle[np].cell = nc10;
2345  Moving2Center (np, nc10);
2346  } else {
2347  if ((nc0 == 0) && (nc10 == 0) && (nc12 != 0)) {
2348  // printf("particle %d moved from boundary cell %d to cell %d \n", np+1, prevcell, nc12);
2349  particle[np].cell = nc12;
2350  Moving2Center (np, nc12);
2351  } else {
2352  if ((nc0 == 0) && (nc10 == 0) && (nc12 == 0) && (ncb != 0)) {
2353  // printf("particle %d moved from boundary cell %d to cell %d \n", np+1, prevcell, ncb);
2354  particle[np].cell = ncb;
2355  Moving2Center (np, ncb);
2356  } else {
2357  if (nc0 + nc10 + nc12 + ncb == 0) {
2358  // printf("moving cell from bound was not found part %d fract %d cell %d k %d\n", np+1, particle[np].fracture, prevcell, k);
2359  }
2360  }
2361  }
2362  }
2363 
2364  return;
2365 }
2367 void ParticleOutput (int currentt, int fract_p)
2369 {
2370  FILE *tmpp;
2371  double posit[3] = {0.0, 0.0, 0.0}, veloc[3] = {0.0, 0.0, 0.0};
2372  double startx, starty, endx, endy, midx, midy, time, velocity_t = 0.0, obeta = 0.0, length_t = 0.0;
2373  int i, tstart = -1, pcell = 0, pfrac = 0, tend = 0, tmid;
2374  int time_l, kdiv = 2;
2375  double eps = 0.05, pressure = 0.0;
2376  struct posit3d particle3dp, particle3dv;
2377  char filename[125];
2378 
2379  if (tfile == 1) {
2380  fclose(tmp);
2381  sprintf(filename, "%s/tempdata_%d", maindir, np);
2382  tmpp = OpenFile(filename, "r");
2383  fscanf(tmpp, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %lf %lf %lf %lf\n ", &tstart, &startx, &starty, &posit[0], &posit[1], &posit[2], &veloc[0], &veloc[1], &veloc[2], &pcell, &pfrac, &time, &obeta, &length_t, &pressure);
2384  } else {
2385  startx = tempdata[0].position2d[0];
2386  starty = tempdata[0].position2d[1];
2387  posit[0] = tempdata[0].position3d[0];
2388  posit[1] = tempdata[0].position3d[1];
2389  posit[2] = tempdata[0].position3d[2];
2390  veloc[0] = tempdata[0].velocity3d[0];
2391  veloc[1] = tempdata[0].velocity3d[1];
2392  veloc[2] = tempdata[0].velocity3d[2];
2393  pcell = tempdata[0].cellp;
2394  pfrac = tempdata[0].fracturep;
2395  time = tempdata[0].timep;
2396  obeta = tempdata[0].betap;
2397  length_t = tempdata[0].length_t;
2398  pressure = tempdata[0].pressure;
2399  tstart = tempdata[0].times;
2400  }
2401 
2402  if (traj_o == 1) {
2403  fprintf(wint, "%5.12E %5.12E %5.12E %5.12E %5.12E %d %5.12E\n", length_t, time, posit[0], posit[1], posit[2], pfrac, obeta);
2404  }
2405 
2406  if (tstart >= 0) {
2407  if (curv_o == 1) {
2408  /* output according to trajectory's curvature */
2409  endx = particle[np].position[0];
2410  endy = particle[np].position[1];
2411  tend = currentt;
2412 
2413  if (tstart != tend) {
2414  if (traj_o == 1) {
2415  fprintf(wv, "%05d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %05d %5.12E %5.12E %5.12E %d %5.12E\n", tstart, posit[0], posit[1], posit[2], veloc[0], veloc[1], veloc[2], pcell, pfrac, time, node[cell[pcell - 1].node_ind[0] - 1].aperture, obeta, 0, pressure);
2416  }
2417 
2418  nodeID++;
2419 
2420  if (avs_o == 1) {
2421  fprintf(wpt, "%05d %5.12E %5.12E %5.12E \n", nodeID, posit[0], posit[1], posit[2]);
2422  velocity_t = sqrt(pow(veloc[0], 2) + pow(veloc[1], 2) + pow(veloc[2], 2));
2423  fprintf(wpt_att, "%010d %06d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E\n", nodeID, pfrac, time, velocity_t, veloc[0], veloc[1], veloc[2], node[cell[pcell - 1].node_ind[0] - 1].aperture, pressure);
2424  }
2425 
2426  int tstep, flag = 0, isch = 0;
2427  double angle_m;
2428 
2429  do {
2430  time_l = tend - tstart;
2431  kdiv = kdiv * 2;
2432  tstep = (int) (time_l / 2.0);
2433 
2434  if (tfile == 1) {
2435  rewind(tmp);
2436  }
2437 
2438  for (isch = 0; isch < time_l; isch++) {
2439  if (tfile == 1) {
2440  fscanf(tmpp, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %lf %lf %lf %lf\n ", &tmid, &midx, &midy, &posit[0], &posit[1], &posit[2], &veloc[0], &veloc[1], &veloc[2], &pcell, &pfrac, &time, &obeta, &length_t, &pressure);
2441  } else {
2442  tmid = tempdata[isch].times;
2443  midx = tempdata[isch].position2d[0];
2444  midy = tempdata[isch].position2d[1];
2445  posit[0] = tempdata[isch].position3d[0];
2446  posit[1] = tempdata[isch].position3d[1];
2447  posit[2] = tempdata[isch].position3d[2];
2448  veloc[0] = tempdata[isch].velocity3d[0];
2449  veloc[1] = tempdata[isch].velocity3d[1];
2450  veloc[2] = tempdata[isch].velocity3d[2];
2451  pcell = tempdata[isch].cellp;
2452  pfrac = tempdata[isch].fracturep;
2453  time = tempdata[isch].timep;
2454  obeta = tempdata[isch].betap;
2455  pressure = tempdata[isch].pressure;
2456  length_t = tempdata[isch].length_t;
2457  }
2458 
2459  if (tmid == tstart + tstep) {
2460  break;
2461  }
2462  }
2463 
2464  angle_m = DefineAngle(startx - midx, starty - midy, endx - midx, endy - midy);
2465 
2466  if (((angle_m > -eps) && (angle_m < eps)) || ((angle_m > pi - eps) && (angle_m < pi + eps))) {
2467  flag = 1;
2468  break;
2469  } else {
2470  tend = tmid;
2471  endx = midx;
2472  endy = midy;
2473  }
2474 
2475  if (tend == tstart) {
2476  flag = 1;
2477  break;
2478  }
2479  } while (flag == 0);
2480 
2481  if (kdiv < 10) {
2482  kdiv = kdiv * 4;
2483  }
2484 
2485  time_l = t - tstart;
2486 
2487  if (kdiv != 0) {
2488  tstep = (int)(time_l / kdiv);
2489  }
2490 
2491  if (kdiv > time_l) {
2492  tstep = 2;
2493  kdiv = (int)time_l / 2.0;
2494  }
2495 
2496  if (kdiv == 0) {
2497  kdiv = 2;
2498  tstep = 1;
2499  }
2500 
2501  if (tfile == 1) {
2502  rewind(tmp);
2503  }
2504 
2505  for(i = 0; i < kdiv - 1; i++) {
2506  for (isch = 0; isch < time_l; isch++) {
2507  if (tfile == 1) {
2508  fscanf(tmpp, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %lf %lf %lf %lf\n ", &tmid, &midx, &midy, &posit[0], &posit[1], &posit[2], &veloc[0], &veloc[1], &veloc[2], &pcell, &pfrac, &time, &obeta, &length_t, &pressure);
2509  } else {
2510  tmid = tempdata[isch].times;
2511  midx = tempdata[isch].position2d[0];
2512  midy = tempdata[isch].position2d[1];
2513  posit[0] = tempdata[isch].position3d[0];
2514  posit[1] = tempdata[isch].position3d[1];
2515  posit[2] = tempdata[isch].position3d[2];
2516  veloc[0] = tempdata[isch].velocity3d[0];
2517  veloc[1] = tempdata[isch].velocity3d[1];
2518  veloc[2] = tempdata[isch].velocity3d[2];
2519  pcell = tempdata[isch].cellp;
2520  pfrac = tempdata[isch].fracturep;
2521  time = tempdata[isch].timep;
2522  obeta = tempdata[isch].betap;
2523  pressure = tempdata[isch].pressure;
2524  length_t = tempdata[isch].length_t;
2525  }
2526 
2527  if (tmid == tstart + tstep * (i + 1)) {
2528  break;
2529  }
2530  }
2531 
2532  if (traj_o == 1) {
2533  fprintf(wv, "%05d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %05d %5.12E %5.12E %5.12E %d %5.12E\n", tmid, posit[0], posit[1], posit[2], veloc[0], veloc[1], veloc[2], pcell, pfrac, time, node[cell[pcell - 1].node_ind[0] - 1].aperture, obeta, 0, pressure);
2534  }
2535 
2536  nodeID++;
2537 
2538  if (avs_o == 1) {
2539  fprintf(wpt, "%05d %5.12E %5.12E %5.12E \n", nodeID, posit[0], posit[1], posit[2]);
2540  velocity_t = sqrt(pow(veloc[0], 2) + pow(veloc[1], 2) + pow(veloc[2], 2));
2541  fprintf(wpt_att, "%010d %06d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E\n", nodeID, pfrac, time, velocity_t, veloc[0], veloc[1], veloc[2], node[cell[pcell - 1].node_ind[0] - 1].aperture, pressure);
2542  }
2543  }
2544  }
2545  } else {
2546  /* in case of every time step output */
2547  if (curv_o != 1) {
2548  if (tfile == 1) {
2549  fscanf(tmpp, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %lf %lf %lf %lf\n ", &tstart, &startx, &starty, &posit[0], &posit[1], &posit[2], &veloc[0], &veloc[1], &veloc[2], &pcell, &pfrac, &time, &obeta, &length_t, &pressure);
2550  rewind(tmp);
2551  time_l = currentt - tstart;
2552  } else {
2553  tstart = tempdata[0].times;
2554  time_l = currentt - tstart;
2555  }
2556 
2557  for (i = 0; i < time_l; i++) {
2558  if (tfile == 1) {
2559  fscanf(tmpp, "%d %lf %lf %lf %lf %lf %lf %lf %lf %d %d %lf %lf %lf %lf\n ", &tstart, &startx, &starty, &posit[0], &posit[1], &posit[2], &veloc[0], &veloc[1], &veloc[2], &pcell, &pfrac, &time, &obeta, &length_t, &pressure);
2560  } else {
2561  tstart = tempdata[i].times;
2562  startx = tempdata[i].position2d[0];
2563  starty = tempdata[i].position2d[1];
2564  posit[0] = tempdata[i].position3d[0];
2565  posit[1] = tempdata[i].position3d[1];
2566  posit[2] = tempdata[i].position3d[2];
2567  veloc[0] = tempdata[i].velocity3d[0];
2568  veloc[1] = tempdata[i].velocity3d[1];
2569  veloc[2] = tempdata[i].velocity3d[2];
2570  pcell = tempdata[i].cellp;
2571  pfrac = tempdata[i].fracturep;
2572  time = tempdata[i].timep;
2573  obeta = tempdata[i].betap;
2574  pressure = tempdata[i].pressure;
2575  length_t = tempdata[i].length_t;
2576  }
2577 
2578  nodeID++;
2579 
2580  if (avs_o == 1) {
2581  fprintf(wpt, "%05d %5.12E %5.12E %5.12E \n", nodeID, posit[0], posit[1], posit[2]);
2582  velocity_t = sqrt(pow(veloc[0], 2) + pow(veloc[1], 2) + pow(veloc[2], 2));
2583  fprintf(wpt_att, "%010d %06d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E\n", nodeID, pfrac, time, velocity_t, veloc[0], veloc[1], veloc[2], node[cell[pcell - 1].node_ind[0] - 1].aperture, pressure);
2584  }
2585 
2586  if (traj_o == 1) {
2587  fprintf(wv, "%05d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %05d %5.12E %5.12E %5.12E %d %5.12E\n", tstart, posit[0], posit[1], posit[2], veloc[0], veloc[1], veloc[2], pcell, pfrac, time, node[cell[pcell - 1].node_ind[0] - 1].aperture, obeta, 0, pressure);
2588  }
2589  }
2590  }
2591  }
2592 
2593  particle3dp = CalculatePosition3D();
2594  particle3dv = CalculateVelocity3D();
2595 
2596  if (traj_o == 1) {
2597  fprintf(wv, "%05d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %05d %05d %5.12E %5.12E %5.12E %d %5.12E\n", t, particle3dp.cord3[0], particle3dp.cord3[1], particle3dp.cord3[2], particle3dv.cord3[0], particle3dv.cord3[1], particle3dv.cord3[2], particle[np].cell, particle[np].fracture, particle[np].time, node[cell[pcell - 1].node_ind[0] - 1].aperture, obeta, fract_p, particle[np].pressure);
2598  }
2599 
2600  nodeID++;
2601 
2602  if (avs_o == 1) {
2603  fprintf(wpt, "%05d %5.12E %5.12E %5.12E \n", nodeID, particle3dp.cord3[0], particle3dp.cord3[1], particle3dp.cord3[2]);
2604  velocity_t = sqrt(pow(particle3dv.cord3[0], 2) + pow(particle3dv.cord3[1], 2) + pow(particle3dv.cord3[2], 2));
2605  fprintf(wpt_att, "%010d %06d %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E %5.12E\n", nodeID, particle[np].fracture, particle[np].time, velocity_t, particle3dv.cord3[0], particle3dv.cord3[1], particle3dv.cord3[2], node[cell[pcell - 1].node_ind[0] - 1].aperture, particle[np].pressure);
2606  }
2607  }
2608 
2609  if (tfile == 1) {
2610  fclose(tmpp);
2611  int status;
2612  sprintf(filename, "%s/tempdata_%d", maindir, np);
2613  status = remove(filename);
2614 
2615  if (no_out != 1) {
2616  if( status != 0 ) {
2617  printf("Unable to delete the file %s\n", filename);
2618  perror("Error");
2619  }
2620 
2621  {
2622  sprintf(filename, "%s/tempdata_%d", maindir, np);
2623  tmp = OpenFile(filename, "w");
2624  }
2625  }
2626  } else {
2627  timecounter = 0;
2628  free(tempdata);
2629  }
2630 
2631  return;
2632 }
2634 void FinalPosition()
2636 {
2637  int n1, n2, n3;
2638  int n1out = 0, n2out = 0;
2639  n1 = cell[particle[np].cell - 1].node_ind[0];
2640  n2 = cell[particle[np].cell - 1].node_ind[1];
2641  n3 = cell[particle[np].cell - 1].node_ind[2];
2642  double cx1, cx2, cy1, cy2, px1, px2, py1, py2, ap, bp, cp, as, bs, cs, deter, xint, yint;
2643 
2644  if ((node[n1 - 1].typeN == 210) || (node[n1 - 1].typeN == 212) || (node[n1 - 1].typeN == 200) || (node[n1 - 1].typeN == 202)) {
2645  n1out = n1;
2646  }
2647 
2648  if ((node[n2 - 1].typeN == 210) || (node[n2 - 1].typeN == 212) || (node[n2 - 1].typeN == 200) || (node[n2 - 1].typeN == 202)) {
2649  if (n1out == 0) {
2650  n2out = n2;
2651  } else {
2652  n1out = n2;
2653  }
2654  }
2655 
2656  if ((node[n3 - 1].typeN == 210) || (node[n3 - 1].typeN == 212) || (node[n3 - 1].typeN == 200) || (node[n3 - 1].typeN == 202)) {
2657  if (n2out == 0) {
2658  n2out = n3;
2659  }
2660  }
2661 
2662  if ((n1out != 0) && (n2out != 0)) {
2663  // two vertices of particles cell are on out flow boundary
2664  double cx1, cx2, cy1, cy2, px1, px2, py1, py2, ap, bp, cp, as, bs, cs, deter, xint, yint;
2665  cx1 = node[n1out - 1].coord_xy[Xindex(n1out, np)];
2666  cy1 = node[n1out - 1].coord_xy[Yindex(n1out, np)];
2667  cx2 = node[n2out - 1].coord_xy[Xindex(n2out, np)];
2668  cy2 = node[n2out - 1].coord_xy[Yindex(n2out, np)];
2669  as = cy2 - cy1;
2670  bs = cx1 - cx2;
2671  cs = as * cx1 + bs * cy1;
2672  px1 = particle[np].position[0];
2673  py1 = particle[np].position[1];
2674  px2 = particle[np].prev_pos[0];
2675  py2 = particle[np].prev_pos[1];
2676  ap = py2 - py1;
2677  bp = px1 - px2;
2678  cp = ap * px1 + bp * py1;
2679  deter = ap * bs - as * bp;
2680  xint = (bs * cp - bp * cs) / deter;
2681  yint = (ap * cs - as * cp) / deter;
2682  double distance = 0, tfinal;
2683  distance = sqrt((xint - px1) * (xint - px1) + (yint - py1) * (yint - py1));
2684  tfinal = distance / sqrt(particle[np].velocity[0] * particle[np].velocity[0] + particle[np].velocity[1] * particle[np].velocity[1]);
2685  particle[np].prev_pos[0] = particle[np].position[0];
2686  particle[np].prev_pos[1] = particle[np].position[1];
2687  particle[np].position[0] = xint;
2688  particle[np].position[1] = yint;
2689  particle[np].time = tfinal + particle[np].time;
2690  // printf("x %12.5e y %12.5e z %12.5e %d time %5.12e\n", particle3dposit.cord3[0], particle3dposit.cord3[1],particle3dposit.cord3[2], particle[np].cell, particle[np].time);
2691  } else {
2692  int ncent = 0, nnext = 0;
2693  ncent = n1out + n2out;
2694 
2695  if (ncent != 0) {
2696  int ii = 0;
2697 
2698  for (ii = 0; ii < node[ncent - 1].numneighb; ii++) {
2699  if ((node[node[ncent - 1].indnodes[ii] - 1].typeN == 212) || (node[node[ncent - 1].indnodes[ii] - 1].typeN == 210) || (node[node[ncent - 1].indnodes[ii] - 1].typeN == 202)) {
2700  nnext = node[ncent - 1].indnodes[ii];
2701  int newcell = 0;
2702  newcell = node[ncent - 1].cells[ii][0];
2703  cx1 = node[ncent - 1].coord_xy[Xindex(ncent, np)];
2704  cy1 = node[ncent - 1].coord_xy[Yindex(ncent, np)];
2705  cx2 = node[nnext - 1].coord_xy[Xindex(nnext, np)];
2706  cy2 = node[nnext - 1].coord_xy[Yindex(nnext, np)];
2707  as = cy2 - cy1;
2708  bs = cx1 - cx2;
2709  cs = as * cx1 + bs * cy1;
2710  px1 = particle[np].position[0];
2711  py1 = particle[np].position[1];
2712  px2 = particle[np].prev_pos[0];
2713  py2 = particle[np].prev_pos[1];
2714  ap = py2 - py1;
2715  bp = px1 - px2;
2716  cp = ap * px1 + bp * py1;
2717  deter = ap * bs - as * bp;
2718  xint = (bs * cp - bp * cs) / deter;
2719  yint = (ap * cs - as * cp) / deter;
2720  double distance = 0, tfinal;
2721  distance = sqrt((xint - px1) * (xint - px1) + (yint - py1) * (yint - py1));
2722  tfinal = distance / sqrt(particle[np].velocity[0] * particle[np].velocity[0] + particle[np].velocity[1] * particle[np].velocity[1]);
2723  particle[np].prev_pos[0] = particle[np].position[0];
2724  particle[np].prev_pos[1] = particle[np].position[1];
2725  particle[np].position[0] = xint;
2726  particle[np].position[1] = yint;
2727  particle[np].time = tfinal + particle[np].time;
2728  particle[np].cell = newcell;
2729  // printf("INSIDE! x %12.5e y %12.5e z %12.5e %d time %5.12e\n", particle3dposit.cord3[0], particle3dposit.cord3[1],particle3dposit.cord3[2], newcell, particle[np].time);
2730  }
2731 
2732  if (nnext != 0) {
2733  break;
2734  }
2735  }
2736  }
2737  }
2738 
2739  return;
2740 }
2742 struct lagrangian CalculateLagrangian(double xcurrent, double ycurrent, double zcurrent, double xprev, double yprev, double zprev)
2744 {
2745  struct lagrangian lagvariable;
2746  struct posit3d particle3dv;
2747  double currentdistance = 0.0, deltatau = 0.0, deltabeta = 0.0, velsquare = 0.0;
2748  particle3dv = CalculateVelocity3D();
2749  currentdistance = pow((xcurrent - xprev), 2) + pow((ycurrent - yprev), 2) + pow((zcurrent - zprev), 2);
2750  velsquare = (pow(particle3dv.cord3[0], 2) + pow(particle3dv.cord3[1], 2) + pow(particle3dv.cord3[2], 2));
2751  deltatau = currentdistance / velsquare;
2752 
2753  if (currentdistance > 0.0) {
2754  deltabeta = sqrt(currentdistance) / (sqrt(velsquare) * (node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture * 0.5));
2755  lagvariable.tau = sqrt(deltatau);
2756  } else {
2757  deltabeta = 0.0;
2758  lagvariable.tau = 0.0;
2759  }
2760 
2761  lagvariable.betta = deltabeta;
2762  lagvariable.initx = xcurrent;
2763  lagvariable.inity = ycurrent;
2764  lagvariable.initz = zcurrent;
2765  return lagvariable;
2766 }
2768 
2769 double TimeDomainRW (double time_advect)
2782 {
2783  double randomnumber = 0.0;
2784  double term_a = 0;
2785  double b = 0;
2786  double inverse_erfc = 0.0;
2787  double z;
2788  double timediff = 0.0;
2789  unsigned int cnt = 0.0;
2790 
2791  if (particle[np].cell != 0) {
2792  if ((node[cell[particle[np].cell - 1].node_ind[0] - 1].typeN != 2) && (node[cell[particle[np].cell - 1].node_ind[0] - 1].typeN != 12)) {
2793  b = node[cell[particle[np].cell - 1].node_ind[0] - 1].aperture;
2794  } else {
2795  if ((node[cell[particle[np].cell - 1].node_ind[1] - 1].typeN != 2) && (node[cell[particle[np].cell - 1].node_ind[1] - 1].typeN != 12)) {
2796  b = node[cell[particle[np].cell - 1].node_ind[1] - 1].aperture;
2797  } else {
2798  b = node[cell[particle[np].cell - 1].node_ind[2] - 1].aperture;
2799  }
2800  }
2801  } else {
2802  b = node[fracture[particle[np].fracture - 1].firstnode - 1].aperture;
2803  }
2804 
2805  term_a = (tdrw_porosity * sqrt(tdrw_diffcoeff)) / b;
2806  randomnumber = drand48();
2807  z = 1.0 - randomnumber;
2808  // Power expansion of inverse ERFC
2809  inverse_erfc = 0.5 * sqrt(pi) * (z + (pi / 12) * pow(z, 3) + ((7 * pow(pi, 2)) / 480) * pow(z, 5) + ((127 * pow(pi, 3)) / 40320) * pow(z, 7) + ((4369 * pow(pi, 4)) / 5806080) * pow(z, 9) + ((34807 * pow(pi, 5)) / 182476800) * pow(z, 11));
2810  timediff = pow(((term_a * time_advect) / inverse_erfc), 2);
2811 
2812  /* If using rate-limited TDRW, check if the particle diffuses too far into the matrix.
2813  If it does, then limit the time to maximum value.
2814  */
2815  if (tdrw_limited == 1) {
2816  double xdiff = sqrt(2 * tdrw_diffcoeff * timediff);
2817 
2818  if (xdiff > tdrw_lambda) {
2819  timediff = 0.5 * (pow(tdrw_lambda, 2) / tdrw_diffcoeff);
2820  }
2821  }
2822 
2823  return timediff;
2824 }
nnodes
unsigned int nnodes
Definition: main.c:29
String_Compare
int String_Compare(char string1[], char string2[])
Definition: ReadGridInit.c:1089
YindexC
int YindexC(int nodenum, int ii)
Definition: VelocityReconstruction.c:791
timesteps
unsigned long int timesteps
Definition: main.c:48
pi
#define pi
Definition: FuncDef.h:2
OpenFile
FILE * OpenFile(char filen[120], char fileopt[2])
Definition: ReadGridInit.c:664
cell
struct element * cell
Definition: main.c:42
Control_File
struct inpfile Control_File(char fileobject[], int ctr)
Definition: ReadGridInit.c:998
DefineAngle
double DefineAngle(double u1, double u2, double v1, double v2)
Definition: VelocityReconstruction.c:347
maindir
char maindir[125]
Definition: main.c:43
FlowInWeight
void FlowInWeight(int numbpart)
Definition: InitialPartPositions.c:747
InitCell
int InitCell()
Definition: InitialPartPositions.c:542
Control_File_Optional
struct inpfile Control_File_Optional(char fileobject[], int ctr)
Definition: ReadGridInit.c:1045
XindexC
int XindexC(int nodenum, int ii)
Definition: VelocityReconstruction.c:776
timeunit
double timeunit
Definition: main.c:51
OutputMarPlumDisp
void OutputMarPlumDisp(int currentnum, char path[125])
Definition: output.c:19
ChangeFracture
void ChangeFracture(int cell_win)
Definition: RotateFracture.c:138
particle
struct contam * particle
Definition: main.c:41
CalculateVelocity3D
struct posit3d CalculateVelocity3D()
Definition: RotateFracture.c:259
nodezoneout
unsigned int * nodezoneout
Definition: main.c:37
Control_Param
struct inpfile Control_Param(char fileobject[], int ctr)
Definition: ReadGridInit.c:1157
node
struct vertex * node
Definition: main.c:40
CalculatePosition3D
struct posit3d CalculatePosition3D()
Definition: RotateFracture.c:161
fracture
struct material * fracture
Definition: main.c:39
nfract
unsigned int nfract
Definition: main.c:31
flag_w
unsigned int flag_w
Definition: main.c:38
nodezonein
unsigned int * nodezonein
Definition: main.c:36
InitPos
int InitPos()
Definition: InitialPartPositions.c:21
tdrw_diffcoeff
double tdrw_diffcoeff
Definition: TrackingPart.c:32
tfile
unsigned int tfile
Definition: TrackingPart.c:31
frac_o
unsigned int frac_o
Definition: TrackingPart.c:31
CalculateCurrentDT
double CalculateCurrentDT()
Definition: TrackingPart.c:2239
tempdata
struct tempout * tempdata
Definition: TrackingPart.c:56
FLAG_OUT
unsigned int FLAG_OUT
Definition: TrackingPart.c:29
StreamlineRandomSampling
int StreamlineRandomSampling(double products[4], double speedsq[4], int indj, int int1, int indk, int neighborcellind[4], int neighborfracind[4], int prevfrac, int prevcell)
Definition: TrackingPart.c:1901
CalculateWeights
struct intcoef CalculateWeights(int nn1, int nn2, int nn3)
Definition: TrackingPart.c:1250
plumec
unsigned int plumec
Definition: TrackingPart.c:31
CorrectorStep
void CorrectorStep()
Definition: TrackingPart.c:1408
SearchNeighborCells
void SearchNeighborCells(int nn1, int nn2, int nn3)
Definition: TrackingPart.c:1285
nodeID
unsigned int nodeID
Definition: TrackingPart.c:30
Moving2NextCellBound
void Moving2NextCellBound(int prevcell)
Definition: TrackingPart.c:2286
no_out
unsigned int no_out
Definition: TrackingPart.c:30
tdrw_limited
unsigned int tdrw_limited
Definition: TrackingPart.c:31
timecounter
unsigned int timecounter
Definition: TrackingPart.c:31
tdrw
unsigned int tdrw
Definition: TrackingPart.c:30
t_adv0
double t_adv0
Definition: TrackingPart.c:32
t
unsigned int t
Definition: TrackingPart.c:30
traj_o
unsigned int traj_o
Definition: TrackingPart.c:30
curv_o
unsigned int curv_o
Definition: TrackingPart.c:30
tdrw_o
unsigned int tdrw_o
Definition: TrackingPart.c:31
InsideCell
int InsideCell(unsigned int numc)
Definition: TrackingPart.c:1324
AcrossIntersection
void AcrossIntersection(int prevcell, int int1, int int2, int mixing_rule)
Definition: TrackingPart.c:1784
Xindex
int Xindex(int nodenum, int nnp)
Definition: TrackingPart.c:2254
tdrw_lambda
double tdrw_lambda
Definition: TrackingPart.c:32
np
unsigned int np
Definition: TrackingPart.c:30
Moving2Center
void Moving2Center(int nnp, int cellnumber)
Definition: TrackingPart.c:2175
all_out
unsigned int all_out
Definition: TrackingPart.c:29
InOutFlowCell
double InOutFlowCell(int indcell, int int1, double nposx, double nposy)
Definition: TrackingPart.c:1730
ParticleTrack
void ParticleTrack()
Definition: TrackingPart.c:66
marfa
unsigned int marfa
Definition: TrackingPart.c:31
CompleteMixingRandomSampling
int CompleteMixingRandomSampling(double products[4], double speedsq[4], int indj, int int1, int indk)
Definition: TrackingPart.c:2089
CheckDistance
int CheckDistance()
Definition: TrackingPart.c:1474
Moving2NextCell
int Moving2NextCell(int stuck, int k)
Definition: TrackingPart.c:2198
CalculateLagrangian
struct lagrangian CalculateLagrangian(double xcurrent, double ycurrent, double zcurrent, double xprev, double yprev, double zprev)
Definition: TrackingPart.c:2742
avs_o
unsigned int avs_o
Definition: TrackingPart.c:30
lagvariable
struct lagrangian lagvariable
Definition: TrackingPart.c:28
CheckNewCell
void CheckNewCell()
Definition: TrackingPart.c:1111
mixing_rule
unsigned int mixing_rule
Definition: TrackingPart.c:30
TimeDomainRW
double TimeDomainRW(double time_advect)
Definition: TrackingPart.c:2769
tdrw_porosity
double tdrw_porosity
Definition: TrackingPart.c:32
NeighborCells
void NeighborCells(int k)
Definition: TrackingPart.c:1426
ParticleOutput
void ParticleOutput(int currentt, int fract_p)
Definition: TrackingPart.c:2367
disp_o
unsigned int disp_o
Definition: TrackingPart.c:31
DefineTimeStep
void DefineTimeStep()
Definition: TrackingPart.c:1447
t_adv
double t_adv
Definition: TrackingPart.c:32
Yindex
int Yindex(int nodenum, int nnp)
Definition: TrackingPart.c:2270
PredictorStep
void PredictorStep()
Definition: TrackingPart.c:1383
timediff
double timediff
Definition: TrackingPart.c:32
FinalPosition
void FinalPosition()
Definition: TrackingPart.c:2634
contam::fracture
unsigned int fracture
Definition: FuncDef.h:189
contam::time
double time
Definition: FuncDef.h:201
contam::intcell
unsigned int intcell
Definition: FuncDef.h:195
contam::pressure
double pressure
Definition: FuncDef.h:207
contam::weight
double weight[3]
Definition: FuncDef.h:198
contam::velocity
double velocity[2]
Definition: FuncDef.h:180
contam::cell
unsigned int cell
Definition: FuncDef.h:192
contam::t_adv_diff
double t_adv_diff
Definition: FuncDef.h:210
contam::prev_pos
double prev_pos[2]
Definition: FuncDef.h:186
contam::t_diff
double t_diff
Definition: FuncDef.h:213
contam::fl_weight
double fl_weight
Definition: FuncDef.h:204
contam::position
double position[2]
Definition: FuncDef.h:183
element::node_ind
unsigned int node_ind[3]
Definition: FuncDef.h:167
element::veloc_ind
unsigned int veloc_ind[3]
Definition: FuncDef.h:170
element::fracture
unsigned int fracture
Definition: FuncDef.h:164
inpfile::filename
char filename[120]
Definition: InitialPartPositions.c:11
intcoef::weights
double weights[3]
Definition: TrackingPart.c:34
lagrangian::initz
double initz
Definition: TrackingPart.c:24
lagrangian::tau
double tau
Definition: TrackingPart.c:20
lagrangian::inity
double inity
Definition: TrackingPart.c:23
lagrangian::initx
double initx
Definition: TrackingPart.c:22
lagrangian::betta
double betta
Definition: TrackingPart.c:21
material::firstnode
unsigned int firstnode
Definition: FuncDef.h:86
posit3d::cord3
double cord3[3]
Definition: InitialPartPositions.c:17
tempout::position3d
double position3d[3]
Definition: TrackingPart.c:46
tempout::length_t
double length_t
Definition: TrackingPart.c:52
tempout::position2d
double position2d[2]
Definition: TrackingPart.c:45
tempout::times
unsigned int times
Definition: TrackingPart.c:44
tempout::cellp
unsigned int cellp
Definition: TrackingPart.c:48
tempout::betap
double betap
Definition: TrackingPart.c:51
tempout::pressure
double pressure
Definition: TrackingPart.c:53
tempout::fracturep
unsigned int fracturep
Definition: TrackingPart.c:49
tempout::velocity3d
double velocity3d[3]
Definition: TrackingPart.c:47
tempout::timep
double timep
Definition: TrackingPart.c:50
vertex::coord
double coord[3]
Definition: FuncDef.h:115
vertex::coord_xy
double coord_xy[6]
Definition: FuncDef.h:118
vertex::pressure
double pressure
Definition: FuncDef.h:124
vertex::cells
unsigned int ** cells
Definition: FuncDef.h:139
vertex::timestep
double timestep[4]
Definition: FuncDef.h:133
vertex::aperture
double aperture
Definition: FuncDef.h:157
vertex::indnodes
unsigned int * indnodes
Definition: FuncDef.h:136
vertex::velocity
double velocity[4][2]
Definition: FuncDef.h:130
vertex::fracture
unsigned int fracture[2]
Definition: FuncDef.h:112
vertex::numneighb
unsigned int numneighb
Definition: FuncDef.h:127