dfnTrans
Code for Particle Tracking simulations in 3D DFN
RotateFracture.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 
9 
10 struct posit3d
12 {
13  double cord3[3];
14 };
15 
16 
17 void Convertto2d()
19 {
20  printf(" \n Converting 3d nodes coordinates to 2d xy parallel plane \n");
21  double check = 0;
22  float norm = 0;
23  unsigned int i, j, l;
24  float nve[3] = {0, 0, 0};
25  float angle, anglecos, anglesin;
26 
27  // loop over all fractures in DFN mesh, defining rotational matrices 3D-2D
28  for (j = 0; j < nfract; j++) {
29  angle = fracture[j].theta * pi / 180;
30 
31  for (l = 0; l < 3; l++) {
32  fracture[j].rot2mat[0][l] = 0.0;
33  fracture[j].rot2mat[1][l] = 0.0;
34  fracture[j].rot2mat[2][l] = 0.0;
35  fracture[j].rot3mat[0][l] = 0.0;
36  fracture[j].rot3mat[1][l] = 0.0;
37  fracture[j].rot3mat[2][l] = 0.0;
38  }
39 
40  if (angle != 0.0) {
41  anglecos = cos(angle);
42  anglesin = sin(angle);
43  /* normal vector times e3*/
44  nve[0] = fracture[j].nvect_xy[0];
45  nve[1] = fracture[j].nvect_xy[1];
46  nve[2] = 0.0;
47  norm = sqrt(nve[0] * nve[0] + nve[1] * nve[1] + nve[2] * nve[2]);
48  nve[0] = nve[0] / norm;
49  nve[1] = nve[1] / norm;
50  nve[2] = nve[2] / norm;
51  check = 0.0;
52  // define a rotational matrix form
53  fracture[j].rot2mat[0][0] = nve[0] * nve[0] * (1 - anglecos) + anglecos;
54  fracture[j].rot2mat[0][1] = nve[0] * nve[1] * (1 - anglecos) - nve[2] * anglesin;
55  fracture[j].rot2mat[0][2] = nve[0] * nve[2] * (1 - anglecos) + nve[1] * anglesin;
56  fracture[j].rot2mat[1][0] = nve[0] * nve[1] * (1 - anglecos) + nve[2] * anglesin;
57  fracture[j].rot2mat[1][1] = nve[1] * nve[1] * (1 - anglecos) + anglecos;
58  fracture[j].rot2mat[1][2] = nve[1] * nve[2] * (1 - anglecos) - nve[0] * anglesin;
59  fracture[j].rot2mat[2][0] = nve[0] * nve[2] * (1 - anglecos) - nve[1] * anglesin;
60  fracture[j].rot2mat[2][1] = nve[1] * nve[2] * (1 - anglecos) + nve[0] * anglesin;
61  fracture[j].rot2mat[2][2] = nve[2] * nve[2] * (1 - anglecos) + anglecos;
62  } //end if
63  } //loop j
64 
65  /* loop over all nodes in factrure */
66  for (i = 0; i < nnodes; i++) {
67  if (fracture[node[i].fracture[0] - 1].theta != 0.0) {
68  j = node[i].fracture[0] - 1;
69  node[i].coord_xy[0] = fracture[j].rot2mat[0][0] * node[i].coord[0] + fracture[j].rot2mat[0][1] * node[i].coord[1] + fracture[j].rot2mat[0][2] * node[i].coord[2];
70  node[i].coord_xy[1] = fracture[j].rot2mat[1][0] * node[i].coord[0] + fracture[j].rot2mat[1][1] * node[i].coord[1] + fracture[j].rot2mat[1][2] * node[i].coord[2];
71  node[i].coord_xy[2] = fracture[j].rot2mat[2][0] * node[i].coord[0] + fracture[j].rot2mat[2][1] * node[i].coord[1] + fracture[j].rot2mat[2][2] * node[i].coord[2];
72  } else {
73  /* if angle =0 and fracture is parallel to xy plane, we use the same x and y coordinates */
74  node[i].coord_xy[0] = node[i].coord[0];
75  node[i].coord_xy[1] = node[i].coord[1];
76  node[i].coord_xy[2] = node[i].coord[2];
77  }
78 
80  if (node[i].fracture[1] != 0) {
81  if (fracture[node[i].fracture[1] - 1].theta != 0.0) {
82  j = node[i].fracture[1] - 1;
83  node[i].coord_xy[3] = fracture[j].rot2mat[0][0] * node[i].coord[0] + fracture[j].rot2mat[0][1] * node[i].coord[1] + fracture[j].rot2mat[0][2] * node[i].coord[2];
84  node[i].coord_xy[4] = fracture[j].rot2mat[1][0] * node[i].coord[0] + fracture[j].rot2mat[1][1] * node[i].coord[1] + fracture[j].rot2mat[1][2] * node[i].coord[2];
85  node[i].coord_xy[5] = fracture[j].rot2mat[2][0] * node[i].coord[0] + fracture[j].rot2mat[2][1] * node[i].coord[1] + fracture[j].rot2mat[2][2] * node[i].coord[2];
86  } else {
87  node[i].coord_xy[3] = node[i].coord[0];
88  node[i].coord_xy[4] = node[i].coord[1];
89  node[i].coord_xy[5] = node[i].coord[2];
90  }
91  }
92  } //loop i
93 
94  return;
95 }
96 
98 /*** function calculates rotation matrix to convert 2D coordinates into 3D*****/
99 
100 void Convertto3d()
102 {
103  float norm = 0;
104  unsigned int j;
105  float nve[3] = {0, 0, 0};
106  float angle, anglecos, anglesin;
107 
108  for (j = 0; j < nfract; j++) {
109  angle = fracture[j].theta * pi / 180;
110 
111  if (angle != 0.0) {
112  anglecos = cos(angle);
113  anglesin = sin(angle);
114  /* normal vector times e3*/
115  nve[0] = fracture[j].nvect_z[0];
116  nve[1] = fracture[j].nvect_z[1];
117  nve[2] = 0.0;
118  norm = sqrt(nve[0] * nve[0] + nve[1] * nve[1] + nve[2] * nve[2]);
119  nve[0] = nve[0] / norm;
120  nve[1] = nve[1] / norm;
121  nve[2] = nve[2] / norm;
122  /* define a rotational matrix form */
123  fracture[j].rot3mat[0][0] = nve[0] * nve[0] * (1 - anglecos) + anglecos;
124  fracture[j].rot3mat[0][1] = nve[0] * nve[1] * (1 - anglecos) - nve[2] * anglesin;
125  fracture[j].rot3mat[0][2] = nve[0] * nve[2] * (1 - anglecos) + nve[1] * anglesin;
126  fracture[j].rot3mat[1][0] = nve[0] * nve[1] * (1 - anglecos) + nve[2] * anglesin;
127  fracture[j].rot3mat[1][1] = nve[1] * nve[1] * (1 - anglecos) + anglecos;
128  fracture[j].rot3mat[1][2] = nve[1] * nve[2] * (1 - anglecos) - nve[0] * anglesin;
129  fracture[j].rot3mat[2][0] = nve[0] * nve[2] * (1 - anglecos) - nve[1] * anglesin;
130  fracture[j].rot3mat[2][1] = nve[1] * nve[2] * (1 - anglecos) + nve[0] * anglesin;
131  fracture[j].rot3mat[2][2] = nve[2] * nve[2] * (1 - anglecos) + anglecos;
132  }
133  }
134 
135  return;
136 }
138 void ChangeFracture(int cell_win)
140 {
141  int j;
142  struct posit3d particle3dposit;
143  particle3dposit = CalculatePosition3D(particle3dposit);
144  j = cell[cell_win - 1].fracture - 1;
145 
146  if (fracture[j].theta != 0.0) {
147  particle[np].position[0] = fracture[j].rot2mat[0][0] * particle3dposit.cord3[0] + fracture[j].rot2mat[0][1] * particle3dposit.cord3[1] + fracture[j].rot2mat[0][2] * particle3dposit.cord3[2];
148  particle[np].position[1] = fracture[j].rot2mat[1][0] * particle3dposit.cord3[0] + fracture[j].rot2mat[1][1] * particle3dposit.cord3[1] + fracture[j].rot2mat[1][2] * particle3dposit.cord3[2];
149  } else {
150  particle[np].position[0] = particle3dposit.cord3[0];
151  particle[np].position[1] = particle3dposit.cord3[1];
152  }
153 
154  particle[np].fracture = cell[cell_win - 1].fracture;
155  particle[np].cell = cell_win;
156  return;
157 }
159 
160 
161 struct posit3d CalculatePosition3D()
163 {
164  int j;
165  double thirdcoord = 0.0;
166  struct posit3d particle3dposit;
167  j = particle[np].fracture - 1;
168 
169  if (fracture[j].theta != 0.0) {
170  if (node[fracture[j].firstnode - 1].fracture[0] == particle[np].fracture) {
171  thirdcoord = node[fracture[j].firstnode - 1].coord_xy[2];
172  } else if (node[fracture[j].firstnode - 1].fracture[1] == particle[np].fracture) {
173  thirdcoord = node[fracture[j].firstnode - 1].coord_xy[5];
174  }
175 
176  particle3dposit.cord3[0] = fracture[j].rot3mat[0][0] * particle[np].position[0] + fracture[j].rot3mat[0][1] * particle[np].position[1] + fracture[j].rot3mat[0][2] * thirdcoord;
177  particle3dposit.cord3[1] = fracture[j].rot3mat[1][0] * particle[np].position[0] + fracture[j].rot3mat[1][1] * particle[np].position[1] + fracture[j].rot3mat[1][2] * thirdcoord;
178  particle3dposit.cord3[2] = fracture[j].rot3mat[2][0] * particle[np].position[0] + fracture[j].rot3mat[2][1] * particle[np].position[1] + fracture[j].rot3mat[2][2] * thirdcoord;
179  } else {
180  particle3dposit.cord3[0] = particle[np].position[0];
181  particle3dposit.cord3[1] = particle[np].position[1];
182  particle3dposit.cord3[2] = node[fracture[j].firstnode - 1].coord[2];
183  }
184 
185  return particle3dposit;
186 }
187 
189 
190 void Velocity3D()
192 {
193  char filename[125];
194  sprintf(filename, "%s/Velocity3D", maindir);
195  FILE *w3 = OpenFile (filename, "w");
196  printf("\n Output flow field: Darcy velocities on nodes in 3D \n");
197  fprintf(w3, " node ID number; X-, Y-, Z- node's location; Vx, Vy, Vz of velocity, control cell volume, fracture number \n");
198  double cord3[3] = {0, 0, 0};
199  unsigned int i, j, v;
200 
201  /* loop over all nodes in factrure */
202  for (i = 0; i < nnodes; i++) {
203  j = node[i].fracture[0] - 1;
204 
205  if ((node[i].typeN != 2) && (node[i].typeN != 12)) {
206  if (fracture[j].theta != 0) {
207  // velocity converter
208  cord3[0] = fracture[j].rot3mat[0][0] * node[i].velocity[0][0] + fracture[j].rot3mat[0][1] * node[i].velocity[0][1] + fracture[j].rot3mat[0][2] * 0;
209  cord3[1] = fracture[j].rot3mat[1][0] * node[i].velocity[0][0] + fracture[j].rot3mat[1][1] * node[i].velocity[0][1] + fracture[j].rot3mat[1][2] * 0;
210  cord3[2] = fracture[j].rot3mat[2][0] * node[i].velocity[0][0] + fracture[j].rot3mat[2][1] * node[i].velocity[0][1] + fracture[j].rot3mat[2][2] * 0;
211  } else {
212  cord3[0] = node[i].velocity[0][0];
213  cord3[1] = node[i].velocity[0][1];
214  cord3[2] = 0;
215  }
216 
217  fprintf(w3, " %05d %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e %d\n", i + 1, node[i].coord[0], node[i].coord[1], node[i].coord[2], cord3[0], cord3[1], cord3[2], node[i].pvolume, node[i].fracture[0]);
218  } else {
219  for (v = 0; v < 2; v++) {
220  if (fracture[j].theta != 0) {
221  // velocity converter
222  cord3[0] = fracture[j].rot3mat[0][0] * node[i].velocity[v][0] + fracture[j].rot3mat[0][1] * node[i].velocity[v][1] + fracture[j].rot3mat[0][2] * 0;
223  cord3[1] = fracture[j].rot3mat[1][0] * node[i].velocity[v][0] + fracture[j].rot3mat[1][1] * node[i].velocity[v][1] + fracture[j].rot3mat[1][2] * 0;
224  cord3[2] = fracture[j].rot3mat[2][0] * node[i].velocity[v][0] + fracture[j].rot3mat[2][1] * node[i].velocity[v][1] + fracture[j].rot3mat[2][2] * 0;
225  } else {
226  cord3[0] = node[i].velocity[v][0];
227  cord3[1] = node[i].velocity[v][1];
228  cord3[2] = 0;
229  }
230 
231  fprintf(w3, " %05d %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e %d\n", i + 1, node[i].coord[0], node[i].coord[1], node[i].coord[2], cord3[0], cord3[1], cord3[2], node[i].pvolume, node[i].fracture[0] );
232  }
233  }
234 
235  if ((node[i].fracture[1] != 0) && ((node[i].typeN == 2) || (node[i].typeN == 12))) {
236  j = node[i].fracture[1] - 1;
237 
238  for (v = 2; v < 4; v++) {
239  if (fracture[j].theta != 0) {
240  // velocity converter
241  cord3[0] = fracture[j].rot3mat[0][0] * node[i].velocity[v][0] + fracture[j].rot3mat[0][1] * node[i].velocity[v][1] + fracture[j].rot3mat[0][2] * 0;
242  cord3[1] = fracture[j].rot3mat[1][0] * node[i].velocity[v][0] + fracture[j].rot3mat[1][1] * node[i].velocity[v][1] + fracture[j].rot3mat[1][2] * 0;
243  cord3[2] = fracture[j].rot3mat[2][0] * node[i].velocity[v][0] + fracture[j].rot3mat[2][1] * node[i].velocity[v][1] + fracture[j].rot3mat[2][2] * 0;
244  } else {
245  cord3[0] = node[i].velocity[v][0];
246  cord3[1] = node[i].velocity[v][1];
247  cord3[2] = 0;
248  }
249 
250  fprintf(w3, " %05d %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e %d\n", i + 1, node[i].coord[0], node[i].coord[1], node[i].coord[2], cord3[0], cord3[1], cord3[2], node[i].pvolume, node[i].fracture[1] );
251  }
252  }
253  }
254 
255  fclose(w3);
256  return;
257 }
259 struct posit3d CalculateVelocity3D()
261 {
262  int j = particle[np].fracture - 1;
263  struct posit3d particle3dvelocity;
264 
265  // printf(" fracture %d theat %lf \n", j, fracture[j].theta);
266 
267  if (fracture[j].theta != 0) {
268  particle3dvelocity.cord3[0] = fracture[j].rot3mat[0][0] * particle[np].velocity[0] + fracture[j].rot3mat[0][1] * particle[np].velocity[1] + fracture[j].rot3mat[0][2] * 0;
269  particle3dvelocity.cord3[1] = fracture[j].rot3mat[1][0] * particle[np].velocity[0] + fracture[j].rot3mat[1][1] * particle[np].velocity[1] + fracture[j].rot3mat[1][2] * 0;
270  particle3dvelocity.cord3[2] = fracture[j].rot3mat[2][0] * particle[np].velocity[0] + fracture[j].rot3mat[2][1] * particle[np].velocity[1] + fracture[j].rot3mat[2][2] * 0;
271  } else {
272  particle3dvelocity.cord3[0] = particle[np].velocity[0];
273  particle3dvelocity.cord3[1] = particle[np].velocity[1];
274  particle3dvelocity.cord3[2] = 0;
275  }
276 
277  return particle3dvelocity;
278 }
280 void Coordinations2D ()
283 {
284  printf(" \n Output 2D coordinates of nodes into files \n");
285  unsigned int inode, ifract;
286  double mincx[nfract], maxcx[nfract], mincy[nfract], maxcy[nfract];
287 
288  //first, find min and max of coordinations ofr every fracture
289  for (ifract = 0; ifract < nfract; ifract++) {
290  mincx[ifract] = 100000;
291  maxcx[ifract] = -100000;
292  mincy[ifract] = 100000;
293  maxcy[ifract] = -100000;
294 
295  for (inode = 0; inode < nnodes; inode++) {
296  if (node[inode].fracture[0] == ifract + 1) {
297  if (node[inode].coord_xy[0] < mincx[ifract]) {
298  mincx[ifract] = node[inode].coord_xy[0];
299  }
300 
301  if (node[inode].coord_xy[1] < mincy[ifract]) {
302  mincy[ifract] = node[inode].coord_xy[1];
303  }
304 
305  if (node[inode].coord_xy[0] > maxcx[ifract]) {
306  maxcx[ifract] = node[inode].coord_xy[0];
307  }
308 
309  if (node[inode].coord_xy[1] > maxcy[ifract]) {
310  maxcy[ifract] = node[inode].coord_xy[1];
311  }
312  } else {
313  if (node[inode].fracture[1] == ifract + 1) {
314  if (node[inode].coord_xy[3] < mincx[ifract]) {
315  mincx[ifract] = node[inode].coord_xy[3];
316  }
317 
318  if (node[inode].coord_xy[4] < mincy[ifract]) {
319  mincy[ifract] = node[inode].coord_xy[4];
320  }
321 
322  if (node[inode].coord_xy[3] > maxcx[ifract]) {
323  maxcx[ifract] = node[inode].coord_xy[3];
324  }
325 
326  if (node[inode].coord_xy[4] > maxcy[ifract]) {
327  maxcy[ifract] = node[inode].coord_xy[4];
328  }
329  }
330  }
331  }
332  }
333 
334  char filename[15];
335  double lengthx, lengthy;
336  FILE *fr;
337  mkdir("Coord2D", 0777);
338 
339  for (ifract = 0; ifract < nfract; ifract++) {
340  sprintf(filename, "Coord2D/coord_%d.dat", ifract + 1);
341  fr = OpenFile(filename, "w");
342  lengthx = fabs(maxcx[ifract] - mincx[ifract]);
343  lengthy = fabs(maxcy[ifract] - mincy[ifract]);
344 
345  for (inode = 0; inode < nnodes; inode++) {
346  if (node[inode].fracture[0] == ifract + 1) {
347  fprintf(fr, "%d %5.8e %5.8e %5.8e %5.8e %5.8e\n", inode + 1, node[inode].coord_xy[0], node[inode].coord_xy[1], (node[inode].coord_xy[0] - mincx[ifract]) / lengthx, (node[inode].coord_xy[1] - mincy[ifract]) / lengthy, node[inode].aperture);
348  } else if (node[inode].fracture[1] == ifract + 1) {
349  fprintf(fr, "%d %5.8e %5.8e %5.8e %5.8e %5.8e\n", inode + 1, node[inode].coord_xy[3], node[inode].coord_xy[4], (node[inode].coord_xy[3] - mincx[ifract]) / lengthx, (node[inode].coord_xy[4] - mincy[ifract]) / lengthy, node[inode].aperture);
350  }
351  }
352 
353  fclose(fr);
354  }
355 
356  return;
357 }
nnodes
unsigned int nnodes
Definition: main.c:29
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
maindir
char maindir[125]
Definition: main.c:43
np
unsigned int np
Definition: main.c:34
particle
struct contam * particle
Definition: main.c:41
node
struct vertex * node
Definition: main.c:40
fracture
struct material * fracture
Definition: main.c:39
nfract
unsigned int nfract
Definition: main.c:31
Coordinations2D
void Coordinations2D()
Definition: RotateFracture.c:280
Convertto2d
void Convertto2d()
Definition: RotateFracture.c:17
ChangeFracture
void ChangeFracture(int cell_win)
Definition: RotateFracture.c:138
CalculateVelocity3D
struct posit3d CalculateVelocity3D()
Definition: RotateFracture.c:259
Velocity3D
void Velocity3D()
Definition: RotateFracture.c:190
CalculatePosition3D
struct posit3d CalculatePosition3D()
Definition: RotateFracture.c:161
Convertto3d
void Convertto3d()
Definition: RotateFracture.c:100
contam::fracture
unsigned int fracture
Definition: FuncDef.h:189
contam::velocity
double velocity[2]
Definition: FuncDef.h:180
contam::cell
unsigned int cell
Definition: FuncDef.h:192
contam::position
double position[2]
Definition: FuncDef.h:183
element::fracture
unsigned int fracture
Definition: FuncDef.h:164
material::rot2mat
double rot2mat[3][3]
Definition: FuncDef.h:98
material::nvect_xy
float nvect_xy[2]
Definition: FuncDef.h:80
material::firstnode
unsigned int firstnode
Definition: FuncDef.h:86
material::rot3mat
double rot3mat[3][3]
Definition: FuncDef.h:101
material::nvect_z
float nvect_z[2]
Definition: FuncDef.h:83
material::theta
float theta
Definition: FuncDef.h:77
posit3d::cord3
double cord3[3]
Definition: InitialPartPositions.c:17
vertex::coord
double coord[3]
Definition: FuncDef.h:115
vertex::coord_xy
double coord_xy[6]
Definition: FuncDef.h:118
vertex::velocity
double velocity[4][2]
Definition: FuncDef.h:130
vertex::fracture
unsigned int fracture[2]
Definition: FuncDef.h:112