dfnTrans
Code for Particle Tracking simulations in 3D DFN
VelocityReconstruction.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 struct lb {
10  double length_b;
11  double angle;
12  double norm[2];
13 };
14 
15 struct matr {
16  double matrinvG[2][40];
17  double matinv[2][2];
18 };
19 
20 struct inpfile {
21  char filename[120];
22  long int flag;
23  double param;
24 };
25 
26 
27 
28 
29 void DarcyVelocity()
37 {
38  printf("\n Darcy's velocities reconstruction \n");
39  double normxarea11[max_neighb - 1][2];
40  unsigned int i, j, l, k1, k2;
41  unsigned long int fracture1 = 0, fracture2 = 0;
42  unsigned int fract_j1[max_neighb];
43  unsigned int fract_j2[max_neighb];
44  double length = 1.0;
45  struct lb lbound;
46  unsigned short int flag1 = 0, flag2 = 0;
47 
48  for (i = 0; i < nnodes; i++) {
49  for (j = 0; j < 4; j++) {
50  node[i].velocity[j][0] = 0.;
51  node[i].velocity[j][1] = 0.;
52  }
53 
54  if ((node[i].typeN == 0) || (node[i].typeN == 310) || (node[i].typeN == 210) || (node[i].typeN == 300) || (node[i].typeN == 200))
55  /* velocity reconstruction for interior nodes */
56  /* 310 is type of exterior nodes in flow-in zone - as interior node */
57  /* 210 is type of exterior nodes in flow-out zone - as interior node */
58  {
59  /* calculating norm to edge times edge's area */
60  for (j = 0; j < node[i].numneighb; j++) {
61  fract_j1[j] = j;
62 
63  if (node[i].fracture[0] == node[node[i].indnodes[fract_j1[j]] - 1].fracture[0]) {
64  if (pflotran == 1) {
65  length = sqrt(pow(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[0], 2) + pow(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[1], 2));
66  normxarea11[j][0] = -1.0 * (node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[0]) * (node[i].area[j] / (length));
67  normxarea11[j][1] = -1.0 * (node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[1]) * (node[i].area[j] / (length));
68  }
69 
70  if (fehm == 1) {
71  length = sqrt(pow(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[0], 2) + pow(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[1], 2));
72  normxarea11[j][0] = -1.0 * (node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[0]) * (node[i].area[j]);
73  normxarea11[j][1] = -1.0 * (node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[1]) * (node[i].area[j]);
74  }
75  }
76 
77  if (node[i].fracture[0] == node[node[i].indnodes[fract_j1[j]] - 1].fracture[1]) {
78  if (pflotran == 1) {
79  length = sqrt(pow(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[3], 2) + pow(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[4], 2));
80  normxarea11[j][0] = -1.0 * (node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[3]) * (node[i].area[j] / (length));
81  normxarea11[j][1] = -1.0 * (node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[4]) * (node[i].area[j] / (length));
82  }
83 
84  if (fehm == 1) {
85  length = sqrt(pow(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[3], 2) + pow(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[4], 2));
86  normxarea11[j][0] = -1.0 * (node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[3]) * (node[i].area[j]);
87  normxarea11[j][1] = -1.0 * (node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[4]) * (node[i].area[j]);
88  }
89  }
90  }
91 
92  VelocityInteriorNode (normxarea11, i, node[i].numneighb, fract_j1, 0);
93  }
94 
95  /* velocity reconstruction for exterior nodes with Newman b.c.***********/
96  if ((node[i].typeN == 10)) {
97  unsigned short int edge01 = 200, edge02 = 200, s = 0, kk;
98 
99  for (j = 0; j < node[i].numneighb; j++) {
100  fract_j1[j] = j;
101  /* define boundary edges */
102  s = 0;
103 
104  for (kk = 0; kk < 4; kk++) {
105  if (node[i].cells[j][kk] != 0) {
106  s = s + 1;
107  }
108  }
109 
110  if (s == 1) {
111  if (edge01 == 200) {
112  edge01 = j;
113  } else {
114  edge02 = j;
115  }
116  }
117  }
118 
119  /* calculating angle between two boundary edges, norm and length */
120  if ((edge01 != 200) && (edge02 != 200)) {
121  lbound = DefineBoundaryAngle (i, edge01, edge02, node[i].fracture[0], 0);
122  } else {
123  printf(" Two boundary edges for node %d not found ! \n", i + 1);
124  }
125 
126  /* calculating norm to edge times edge's area, G */
127  for (j = 0; j < node[i].numneighb; j++) {
128  if (node[i].fracture[0] == node[node[i].indnodes[j] - 1].fracture[0]) {
129  if (pflotran == 1) {
130  length = sqrt(pow(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[0], 2) + pow(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[1], 2));
131  normxarea11[j][0] = -1.*(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[0]) * (node[i].area[j] / (length));
132  normxarea11[j][1] = -1.*(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[1]) * (node[i].area[j] / (length));
133  }
134 
135  if (fehm == 1) {
136  length = sqrt(pow(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[0], 2) + pow(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[1], 2));
137  normxarea11[j][0] = -1.*(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[0]) * (node[i].area[j]);
138  normxarea11[j][1] = -1.*(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[1]) * (node[i].area[j]);
139  }
140  }
141 
142  if (node[i].fracture[0] == node[node[i].indnodes[j] - 1].fracture[1]) {
143  if (pflotran == 1) {
144  length = sqrt(pow(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[3], 2) + pow(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[4], 2));
145  normxarea11[j][0] = -1.*(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[3]) * (node[i].area[j] / (length));
146  normxarea11[j][1] = -1.*(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[4]) * (node[i].area[j] / (length));
147  }
148 
149  if (fehm == 1) {
150  length = sqrt(pow(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[3], 2) + pow(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[4], 2));
151  normxarea11[j][0] = -1.*(node[i].coord_xy[0] - node[node[i].indnodes[j] - 1].coord_xy[3]) * (node[i].area[j]);
152  normxarea11[j][1] = -1.*(node[i].coord_xy[1] - node[node[i].indnodes[j] - 1].coord_xy[4]) * (node[i].area[j]);
153  }
154  }
155  }
156 
157  VelocityExteriorNode (normxarea11, i, node[i].numneighb, fract_j1, lbound, 0 );
158  }
159 
160  /* velocity reconstruction for nodes on intersection: ***********/
161  /* divide the polygons on 2 parts for each intersecting fracture******/
162  /* then divide each part on two again. ******/
163 
164  if ((node[i].typeN == 12) || (node[i].typeN == 2) || (node[i].typeN == 302) || (node[i].typeN == 202) || (node[i].typeN == 312) || (node[i].typeN == 212)) {
165  /* separating polygons into two: two fractures*/
166  fracture1 = node[i].fracture[0];
167  fracture2 = node[i].fracture[1];
168  fract_j1[0] = 0;
169  fract_j2[0] = 0;
170  k1 = 0;
171  k2 = 0;
172 
173  for (j = 0; j < node[i].numneighb; j++) {
174  l = 0;
175  flag1 = 0;
176  flag2 = 0;
177 
178  while (node[i].fracts[j][l] != 0) {
179  if (node[i].fracts[j][l] == fracture1) {
180  if (flag1 == 0) {
181  fract_j1[k1] = j;
182  k1++;
183  flag1 = 1;
184  }
185  }
186 
187  if (node[i].fracts[j][l] == fracture2)
188  if (flag2 == 0) {
189  fract_j2[k2] = j;
190  k2++;
191  flag2 = 1;
192  }
193 
194  l++;
195  } // loop on l
196  } //loop j
197 
198  HalfPolygonVelocity(i, k1, fracture1, 0, fract_j1);
199  HalfPolygonVelocity(i, k2, fracture2, 2, fract_j2);
200  }
201  } //loop i
202 
203  BoundaryCells();
204  printf(" Velocities on nodes are calculated \n" );
205  // int res;
206  // struct inpfile inputfile;
207  // inputfile=Control_File("out_2dflow:",11);
208  // res=strncmp(inputfile.filename,"yes",3);
209  // if (res==0)
210  // OutputVelocities();
211  return;
212 }
214 void OutputVelocities() {
218  char filename[125];
219  sprintf(filename, "%s/plot_vel", maindir);
220  FILE *wq = fopen (filename, "w");
221  sprintf(filename, "%s/fract_nodes", maindir);
222  FILE *wn = fopen (filename, "w");
223  int f, sumf = 0, i;
224  double speed = 0.;
225  fprintf(wn, " %d \n", nfract);
226 
227  for (f = 1; f < nfract + 1; f++) {
228  sumf = 0;
229 
230  for(i = 0; i < nnodes; i++) {
231  if (node[i].fracture[0] == f) {
232  speed = node[i].velocity[0][0] * node[i].velocity[0][0] + node[i].velocity[0][1] * node[i].velocity[0][1];
233  fprintf(wq, " %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e\n", node[i].coord_xy[0], node[i].coord_xy[1], node[i].velocity[0][0], node[i].velocity[0][1], node[i].pressure, speed);
234  sumf = sumf + 1;
235 
236  if (node[i].velocity[1][0] != 0.) {
237  speed = node[i].velocity[1][0] * node[i].velocity[1][0] + node[i].velocity[1][1] * node[i].velocity[1][1];
238  fprintf(wq, " %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e\n", node[i].coord_xy[0], node[i].coord_xy[1], node[i].velocity[1][0], node[i].velocity[1][1], node[i].pressure, speed);
239  sumf = sumf + 1;
240  }
241  } else {
242  if (node[i].fracture[1] == f) {
243  speed = node[i].velocity[2][0] * node[i].velocity[2][0] + node[i].velocity[2][1] * node[i].velocity[2][1];
244  fprintf(wq, " %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e\n", node[i].coord_xy[3], node[i].coord_xy[4], node[i].velocity[2][0], node[i].velocity[2][1], node[i].pressure, speed);
245  sumf = sumf + 1;
246 
247  if (node[i].velocity[3][0] != 0.) {
248  speed = node[i].velocity[3][0] * node[i].velocity[3][0] + node[i].velocity[3][1] * node[i].velocity[3][1];
249  fprintf(wq, " %5.8e %5.8e %5.8e %5.8e %5.8e %5.8e\n", node[i].coord_xy[3], node[i].coord_xy[4], node[i].velocity[3][0], node[i].velocity[3][1], node[i].pressure, speed);
250  sumf = sumf + 1;
251  }
252  }
253  }
254  }//loop i
255 
256  fprintf(wn, "%d \n", sumf);
257  }// loop f
258 
259  fclose(wq);
260  fclose(wn);
261  return;
262 }
264 
265 struct matr MatrixProducts (double normxarea[][2], int number) {
266  struct matr matrices;
267  double dp[2][2] = {{0, 0}, {0, 0}};
268  unsigned int j, m, n;
269 
270  for (j = 0; j < number; j++) {
271  /* dot product (Gt.G) */
272  for (m = 0; m < 2; m++) {
273  for (n = 0; n < 2; n++) {
274  dp[n][m] = dp[n][m] + normxarea[j][n] * normxarea[j][m];
275  }
276  }
277  }
278 
279  /* calculating Inverse */
280  double determinant = dp[0][0] * dp[1][1] - dp[0][1] * dp[1][0];
281  double eps = 1e-75;
282 
283  if ((determinant < eps) && (determinant > -eps)) {
284  // printf("Det of matrix is close to zero %15.8e \n",determinant);
285  determinant = eps;
286  }
287 
288  double invdet = 1.0 / determinant;
289  matrices.matinv[0][0] = dp[1][1] * invdet;
290  matrices.matinv[1][0] = -dp[1][0] * invdet;
291  matrices.matinv[0][1] = -dp[0][1] * invdet;
292  matrices.matinv[1][1] = dp[0][0] * invdet;
293 
294  /* inverse dot product Gt */
295  for (m = 0; m < 2; m++) {
296  for (j = 0; j < number; j++) {
297  matrices.matrinvG[m][j] = matrices.matinv[m][0] * normxarea[j][0] + matrices.matinv[m][1] * normxarea[j][1];
298  }
299  }
300 
301  return matrices;
302 }
304 struct lb DefineBoundaryAngle(int i, unsigned int edge_1, unsigned int edge_2, int f1, int coorf) {
306  struct lb lbound;
307  double normu, normv;
308  double u[2] = {0, 0}, v[2] = {0, 0};
309 
310  if (node[node[i].indnodes[edge_1] - 1].fracture[0] == f1) {
311  u[0] = node[node[i].indnodes[edge_1] - 1].coord_xy[0] - node[i].coord_xy[coorf];
312  u[1] = node[node[i].indnodes[edge_1] - 1].coord_xy[1] - node[i].coord_xy[coorf + 1];
313  }
314 
315  if (node[node[i].indnodes[edge_1] - 1].fracture[1] == f1) {
316  u[0] = node[node[i].indnodes[edge_1] - 1].coord_xy[3] - node[i].coord_xy[coorf];
317  u[1] = node[node[i].indnodes[edge_1] - 1].coord_xy[4] - node[i].coord_xy[coorf + 1];
318  }
319 
320  if (node[node[i].indnodes[edge_2] - 1].fracture[1] == f1) {
321  v[0] = node[node[i].indnodes[edge_2] - 1].coord_xy[3] - node[i].coord_xy[coorf];
322  v[1] = node[node[i].indnodes[edge_2] - 1].coord_xy[4] - node[i].coord_xy[coorf + 1];
323  }
324 
325  if (node[node[i].indnodes[edge_2] - 1].fracture[0] == f1) {
326  v[0] = node[node[i].indnodes[edge_2] - 1].coord_xy[0] - node[i].coord_xy[coorf];
327  v[1] = node[node[i].indnodes[edge_2] - 1].coord_xy[1] - node[i].coord_xy[coorf + 1];
328  }
329 
330  lbound.angle = DefineAngle(u[0], u[1], v[0], v[1]);
331  normu = sqrt(u[1] * u[1] + u[0] * u[0]);
332  normv = sqrt(v[1] * v[1] + v[0] * v[0]);
333 
334  if (normu > normv) {
335  lbound.norm[0] = -u[1];
336  lbound.norm[1] = u[0];
337  } else {
338  lbound.norm[0] = -v[1];
339  lbound.norm[1] = v[0];
340  }
341 
342  lbound.length_b = normu + normv;
343  return lbound;
344 }
346 
347 double DefineAngle(double u1, double u2, double v1, double v2) {
348  double angle, normu, normv;
349  normu = sqrt(u1 * u1 + u2 * u2);
350  normv = sqrt(v1 * v1 + v2 * v2);
351  angle = acosf(((u1 * v1 + u2 * v2) / (normu * normv)));
352  return angle;
353 }
355 
356 void VelocityInteriorNode (double normx_area[][2], int i, int number, unsigned int indj[max_neighb], int vi )
357 
358 {
359  short int m, k, j;
360  double qhat[2] = {0, 0};
361  struct matr matrices;
362  double massbalance = 0;
363  matrices.matinv[0][0] = 0.0;
364  matrices.matinv[0][1] = 0.0;
365  matrices.matinv[1][0] = 0.0;
366  matrices.matinv[1][1] = 0.0;
367 
368  for (k = 0; k < 40; k++) {
369  matrices.matrinvG[0][k] = 0.0;
370  matrices.matrinvG[1][k] = 0.0;
371  }
372 
373  matrices = MatrixProducts (normx_area, number);
374 
375  for (m = 0; m < 2; m++) {
376  qhat[m] = 0;
377 
378  for (j = 0; j < number; j++) {
379  qhat[m] = qhat[m] + matrices.matrinvG[m][j] * node[i].flux[indj[j]];
380  massbalance = massbalance + node[i].flux[indj[j]];
381  }
382 
383  // velocity is Darcy's velocity qhat / density * porosity and converted to required time units
384  node[i].velocity[vi][m] = (qhat[m] / (density * porosity)) * timeunit;
385  }
386 
387  return;
388 }
390 
391 void VelocityExteriorNode (double norm_xarea[][2], int i, int number, unsigned int indj[max_neighb], struct lb lbound, int vi) {
392  short int m, k, j;
393  double qhat[2] = {0, 0};
394  struct matr matrices;
395  double Bmatr[2] = {0, 0};
396  double q_b[2] = {0, 0}, b1[2] = {0, 0}, inv, bqhat, btw[2] = {0, 0};
397  double massbalance = 0;
398  matrices.matinv[0][0] = 0.0;
399  matrices.matinv[0][1] = 0.0;
400  matrices.matinv[1][0] = 0.0;
401  matrices.matinv[1][1] = 0.0;
402 
403  for (k = 0; k < 40; k++) {
404  matrices.matrinvG[0][k] = 0.0;
405  matrices.matrinvG[1][k] = 0.0;
406  }
407 
408  matrices = MatrixProducts (norm_xarea, number);
409 
410  for (m = 0; m < 2; m++) {
411  qhat[m] = 0;
412 
413  for (j = 0; j < number; j++) {
414  qhat[m] = qhat[m] + matrices.matrinvG[m][j] * node[i].flux[indj[j]];
415  massbalance = massbalance + node[i].flux[indj[j]];
416  }
417 
418  Bmatr[m] = lbound.norm[m] * lbound.length_b;
419  }
420 
421  // B*qhat
422  bqhat = Bmatr[0] * qhat[0] + Bmatr[1] * qhat[1];
423  // (GtG)-1 Bt
424  b1[0] = Bmatr[0] * matrices.matinv[0][0] + Bmatr[1] * matrices.matinv[1][0];
425  b1[1] = Bmatr[0] * matrices.matinv[0][1] + Bmatr[1] * matrices.matinv[1][1];
426  // inv [ Bt dot (GtG)-1 B ]
427  inv = 1. / (b1[0] * Bmatr[0] + b1[1] * Bmatr[1]);
428  // btw= B*inv
429  btw[0] = Bmatr[0] * inv;
430  btw[1] = Bmatr[1] * inv;
431  // GtG*btw
432  b1[0] = matrices.matinv[0][0] * btw[0] + matrices.matinv[0][1] * btw[1];
433  b1[1] = matrices.matinv[1][0] * btw[0] + matrices.matinv[1][1] * btw[1];
434  q_b[0] = qhat[0] - b1[0] * bqhat;
435  q_b[1] = qhat[1] - b1[1] * bqhat;
436 
437  for (m = 0; m < 2; m++) {
438  // velocity is Darcy's velocity qhat / density * porosity and converted to required time units
439  node[i].velocity[vi][m] = (q_b[m] / (density * porosity)) * timeunit;
440  }
441 
442  return;
443 }
445 
447 void BoundaryCells () {
448  long int i, j, n1, n2, n3, celln = 0;
449 
450  for(i = 0; i < ncells; i++) {
451  /* define k1 and k2 - two boundary nodes in a boundary cell */
452  long int k1 = 0, k2 = 0, v1 = 0, v2 = 0, k3 = 0, v3 = 0, node3 = 0, vel3 = 0, kk3 = 0, vv3 = 0;
453  n1 = cell[i].node_ind[0];
454 
455  if ((node[n1 - 1].typeN == 10) || (node[n1 - 1].typeN == 12) || (node[n1 - 1].typeN == 310) || (node[n1 - 1].typeN == 312)) {
456  k1 = n1;
457  v1 = cell[i].veloc_ind[0];
458  }
459 
460  n2 = cell[i].node_ind[1];
461  n3 = cell[i].node_ind[2];
462 
463  if ((node[n2 - 1].typeN == 10) || (node[n2 - 1].typeN == 12) || (node[n2 - 1].typeN == 310) || (node[n2 - 1].typeN == 312)) {
464  if (k1 == 0) {
465  k1 = n2;
466  v1 = cell[i].veloc_ind[1];
467  node3 = n1;
468  vel3 = cell[i].veloc_ind[0];
469  } else {
470  k2 = n2;
471  v2 = cell[i].veloc_ind[1];
472  node3 = n3;
473  vel3 = cell[i].veloc_ind[2];
474  }
475  }
476 
477  if ((node[n3 - 1].typeN == 10) || (node[n3 - 1].typeN == 12) || (node[n3 - 1].typeN == 310) || (node[n3 - 1].typeN == 312)) {
478  if (k1 == 0) {
479  k1 = n3;
480  v1 = cell[i].veloc_ind[2];
481  } else {
482  if (k2 == 0) {
483  k2 = n3;
484  v2 = cell[i].veloc_ind[2];
485  } else {
486  kk3 = n3;
487  vv3 = cell[i].veloc_ind[2];
488  }
489 
490  if (k1 == n1) {
491  node3 = n2;
492  vel3 = cell[i].veloc_ind[1];
493  }
494  }
495  }
496 
497  if ((k1 * k2 != 0)) {
498  /****** first, reconstruct velocities on boundary corners***************/
499  /* first node k1 */
500  long int jk3 = 0, jk2 = 0, jk1 = 0, notcorner = 0, notcorner1 = 0, k;
501 
502  for (j = 0; j < node[k1 - 1].numneighb; j++) {
503  if (((node[k1 - 1].type[j] == 10) || (node[k1 - 1].type[j] == 12) || (node[node[k1 - 1].indnodes[j] - 1].typeN == 310)) && (node[k1 - 1].indnodes[j] != k2)) {
504  for (k = 0; k < 4; k++) {
505  if (node[k1 - 1].fracts[j][k] == cell[i].fracture) {
506  if (node[k1 - 1].cells[j][1] == 0) {
507  k3 = node[k1 - 1].indnodes[j];
508  jk3 = j;
509  celln = node[k1 - 1].cells[j][k];
510  }
511  }
512  }
513  }
514 
515  if (node[k1 - 1].indnodes[j] == k2) {
516  jk2 = j;
517  }
518  }// loop on j
519 
520  if ((k3 != 0) && (node[k1 - 1].cells[jk2][1] == 0) && (node[k1 - 1].cells[jk3][1] == 0))
521  /* consider only boundary cells on corner */
522  {
523  if (cell[celln - 1].node_ind[0] == k3) {
524  v3 = cell[celln - 1].veloc_ind[0];
525  }
526 
527  if (cell[celln - 1].node_ind[1] == k3) {
528  v3 = cell[celln - 1].veloc_ind[1];
529  }
530 
531  if (cell[celln - 1].node_ind[2] == k3) {
532  v3 = cell[celln - 1].veloc_ind[2];
533  }
534 
535  notcorner1 = CornerVelocity(i, k1, k2, k3, v1, v2, v3);
536  }
537 
538  /* check k2, second node on boundary */
539  jk3 = 0;
540  jk1 = 0;
541 
542  for (j = 0; j < node[k2 - 1].numneighb; j++) {
543  if (((node[k2 - 1].type[j] == 10) || (node[k2 - 1].type[j] == 12) || (node[node[k2 - 1].indnodes[j] - 1].typeN == 310)) && (node[k2 - 1].indnodes[j] != k1)) {
544  for (k = 0; k < 4; k++) {
545  if (node[k2 - 1].fracts[j][k] == cell[i].fracture) {
546  if (node[k2 - 1].cells[j][1] == 0) {
547  k3 = node[k2 - 1].indnodes[j];
548  jk3 = j;
549  celln = node[k2 - 1].cells[j][0];
550  }
551  }
552  }
553  }
554 
555  if (node[k2 - 1].indnodes[j] == k1) {
556  jk1 = j;
557  }
558  }// loop on j
559 
560  if ((k3 != 0) && (node[k2 - 1].cells[jk1][1] == 0) && (node[k2 - 1].cells[jk3][1] == 0))
561  // consider only boundary cells on corner
562  {
563  if (cell[celln - 1].node_ind[0] == k3) {
564  v3 = cell[celln - 1].veloc_ind[0];
565  }
566 
567  if (cell[celln - 1].node_ind[1] == k3) {
568  v3 = cell[celln - 1].veloc_ind[1];
569  }
570 
571  if (cell[celln - 1].node_ind[2] == k3) {
572  v3 = cell[celln - 1].veloc_ind[2];
573  }
574 
575  notcorner = CornerVelocity(i, k2, k1, k3, v2, v1, v3);
576  }
577 
578  /**** reconstruct velocities on boundary of fracture, end point of intersection ***/
579  /*** when two boundary velocities are antiparallel and pointing to each other,
580  one of them being turned onto flow direction ****************************/
581  if ((notcorner + notcorner1) > 0) {
582  double vk1[2] = {0, 0}, vk2[2] = {0, 0}, vk3[2] = {0, 0}, velx, vely, angleuv;
583  double angle, eps = 0.0001, node3x = 0, node3y = 0;
584  vk1[0] = node[k1 - 1].velocity[v1][0];
585  vk1[1] = node[k1 - 1].velocity[v1][1];
586  vk2[0] = node[k2 - 1].velocity[v2][0];
587  vk2[1] = node[k2 - 1].velocity[v2][1];
588  vk3[0] = node[k3 - 1].velocity[v3][0];
589  vk3[1] = node[k3 - 1].velocity[v3][1];
590  double angle_1;
591  angle_1 = DefineAngle(vk1[0], vk1[1], vk2[0], vk2[1]);
592 
593  /* check those cells on boundaies but not on intersection */
594  if ((vk1[0]*vk2[0] < 0) && (vk1[1]*vk2[1] < 0) && (angle_1 < (pi + eps)) && (angle_1 > (pi - eps)) && ((node[k1 - 1].typeN + node[k2 - 1].typeN) == 20)) {
595  // printf("corner1, %d %d %lf %lf\n", node[k2-1].fracture[0], k2, node[k2-1].coord_xy[0], node[k2-1].coord_xy[0]);
596  double minpressure = 1000;
597  int ni = 0;
598 
599  for (ni = 0; ni < node[k1 - 1].numneighb; ni++) {
600  if ((node[node[k1 - 1].indnodes[ni] - 1].pressure < minpressure) && (node[k1 - 1].indnodes[ni] != k2)) {
601  minpressure = node[node[k1 - 1].indnodes[ni] - 1].pressure;
602  node3 = node[k1 - 1].indnodes[ni];
603  node3x = node[node3 - 1].coord_xy[XindexC(node3, i)] - node[k1 - 1].coord_xy[0];
604  node3y = node[node3 - 1].coord_xy[YindexC(node3, i)] - node[k1 - 1].coord_xy[1];
605  }
606  }
607 
608  angleuv = DefineAngle(vk1[0], vk1[1], node3x, node3y);
609  velx = vk1[0] * cos(angleuv) - vk1[1] * sin(angleuv);
610  vely = vk1[0] * sin(angleuv) + vk1[1] * cos(angleuv);
611  angle = DefineAngle(velx, vely, node3x, node3y);
612 
613  if (angle > eps) {
614  velx = vk1[0] * cos(angleuv) + vk1[1] * sin(angleuv);
615  vely = -1 * vk1[0] * sin(angleuv) + vk1[1] * cos(angleuv);
616  angle = DefineAngle(velx, vely, node3x, node3y);
617  }
618 
619  if (angle < eps) {
620  node[k1 - 1].velocity[v1][0] = velx;
621  node[k1 - 1].velocity[v1][1] = vely;
622  }
623  }
624 
625  /* check those cells on boundaries and on intersection */
626  if ((vk1[0]*vk2[0] < 0) && (vk1[1]*vk2[1] < 0) && ((node[k1 - 1].typeN + node[k2 - 1].typeN) != 20)) {
627  // printf("corner2, %d %d %lf %lf\n", node[k2-1].fracture[0], k2, node[k2-1].coord_xy[0], node[k2-1].coord_xy[0]);
628  if (node[k1 - 1].typeN == 12) {
629  node[k1 - 1].velocity[v1][0] = (-1) * node[k1 - 1].velocity[v1][0];
630  node[k1 - 1].velocity[v1][1] = (-1) * node[k1 - 1].velocity[v1][1];
631  } else {
632  if (node[k2 - 1].typeN == 12) {
633  node[k2 - 1].velocity[v2][0] = (-1) * node[k2 - 1].velocity[v2][0];
634  node[k2 - 1].velocity[v2][1] = (-1) * node[k2 - 1].velocity[v2][1];
635  }
636  }
637  }
638  }
639  } //if k1*k2
640 
641  if((k1 * k2 != 0) && (kk3 != 0)) {
642  if (node[kk3 - 1].numneighb == 2) {
643  int notc;
644  notc = CornerVelocity(i, kk3, k2, k1, vv3, v2, v1);
645  }
646  }
647  }// loop on i
648 
649  return;
650 }
651 
653 
658 int CornerVelocity(int i, int m1, int m2, int m3, int s1, int s2, int s3) {
659  double v[2] = {0, 0}, u[2] = {0, 0}, vk1[2] = {0, 0}, vk2[2] = {0, 0}, vk3[2] = {0, 0}, ang_uv, ang_vk1u, ang_vk1v, ang_vk1vk2, ang_vk1vk3;
660  int notcorner;
661  // calculate the angle between two boundary edges of the cell i
662  // m1 - is a central node,
663  // u - is an edge (vector) between nodes m1 and m2
664  // v - is an edge between nodes m1 and m3
665  // vk1 - velocity of node m1 on cell i
666  // vk2 - velocity on node m2 on cell i
667  // vk3 - velocity on node m3 on cell i
668  // if angle is 0 or pi between velocity of central node and one of the edges and it's velocity, then velocity of central node is rotated on angle with the second edge
669  u[0] = node[m2 - 1].coord_xy[XindexC(m2, i)] - node[m1 - 1].coord_xy[XindexC(m1, i)];
670  u[1] = node[m2 - 1].coord_xy[YindexC(m2, i)] - node[m1 - 1].coord_xy[YindexC(m1, i)];
671  v[0] = node[m3 - 1].coord_xy[XindexC(m3, i)] - node[m1 - 1].coord_xy[XindexC(m1, i)];
672  v[1] = node[m3 - 1].coord_xy[YindexC(m3, i)] - node[m1 - 1].coord_xy[YindexC(m1, i)];
673  ang_uv = DefineAngle(u[0], u[1], v[0], v[1]);
674  /* if current cell is a corner cell then calculate angles of velocities and edges */
675  double eps = 0.001, angle;
676  double velx, vely;
677 
678  if (((ang_uv >= (pi - eps)) && (ang_uv <= (pi + eps))) || ((ang_uv >= (-eps)) && (ang_uv <= (+eps)))) {
679  notcorner = 1;
680  }
681 
682  if ((ang_uv < (pi - eps)) && (ang_uv > eps)) {
683  notcorner = 0;
684  vk1[0] = node[m1 - 1].velocity[s1][0];
685  vk1[1] = node[m1 - 1].velocity[s1][1];
686  vk2[0] = node[m2 - 1].velocity[s2][0];
687  vk2[1] = node[m2 - 1].velocity[s2][1];
688  vk3[0] = node[m3 - 1].velocity[s3][0];
689  vk3[1] = node[m3 - 1].velocity[s3][1];
690  ang_vk1v = DefineAngle(vk1[0], vk1[1], v[0], v[1]);
691  ang_vk1u = DefineAngle(vk1[0], vk1[1], u[0], u[1]);
692  ang_vk1vk2 = DefineAngle(vk1[0], vk1[1], vk2[0], vk2[1]);
693  ang_vk1vk3 = DefineAngle(vk1[0], vk1[1], vk3[0], vk3[1]);
694 
695  if ((ang_vk1u > (pi - eps)) && (ang_vk1u < (pi + eps)) && ((ang_vk1vk2 < eps) || (ang_vk1vk2 > (pi - eps)))) {
696  if ((node[m1 - 1].typeN < 300) && (node[m3 - 1].typeN < 300)) {
697  velx = node[m1 - 1].velocity[s1][0] * cos(ang_vk1v) - node[m1 - 1].velocity[s1][1] * sin(ang_vk1v);
698  vely = node[m1 - 1].velocity[s1][0] * sin(ang_vk1v) + node[m1 - 1].velocity[s1][1] * cos(ang_vk1v);
699  angle = DefineAngle(velx, vely, v[0], v[1]);
700 
701  if (angle > eps) {
702  velx = node[m1 - 1].velocity[s1][0] * cos(ang_vk1v) + node[m1 - 1].velocity[s1][1] * sin(ang_vk1v);
703  vely = -1 * node[m1 - 1].velocity[s1][0] * sin(ang_vk1v) + node[m1 - 1].velocity[s1][1] * cos(ang_vk1v);
704  angle = DefineAngle(velx, vely, v[0], v[1]);
705  }
706 
707  if (angle < eps) {
708  node[m1 - 1].velocity[s1][0] = velx;
709  node[m1 - 1].velocity[s1][1] = vely;
710  }
711  }
712  } else {
713  if ((ang_vk1v > (pi - eps)) && (ang_vk1v < (pi + eps)) && ((ang_vk1vk3 < eps) || (ang_vk1vk3 > (pi - eps)))) {
714  if ((node[m1 - 1].typeN < 300) && (node[m2 - 1].typeN < 300)) {
715  velx = node[m1 - 1].velocity[s1][0] * cos(ang_vk1u) - node[m1 - 1].velocity[s1][1] * sin(ang_vk1u);
716  vely = node[m1 - 1].velocity[s1][0] * sin(ang_vk1u) + node[m1 - 1].velocity[s1][1] * cos(ang_vk1u);
717  angle = DefineAngle(velx, vely, u[0], u[1]);
718 
719  if (angle > eps) {
720  velx = node[m1 - 1].velocity[s1][0] * cos(ang_vk1u) + node[m1 - 1].velocity[s1][1] * sin(ang_vk1u);
721  vely = -1 * node[m1 - 1].velocity[s1][0] * sin(ang_vk1u) + node[m1 - 1].velocity[s1][1] * cos(ang_vk1u);
722  angle = DefineAngle(velx, vely, u[0], u[1]);
723  }
724 
725  if (angle < eps) {
726  node[m1 - 1].velocity[s1][0] = velx;
727  node[m1 - 1].velocity[s1][1] = vely;
728  }
729  }
730  }
731  }
732 
733  /* if node is on inflow boundary:*/
734  double ang_u, ang_v;
735 
736  if (((node[m1 - 1].typeN == 310) || (node[m1 - 1].typeN == 312)) && (node[m2 - 1].typeN != 310) && (node[m2 - 1].typeN != 312)) {
737  ang_u = DefineAngle(vk1[0], vk1[1], u[0], u[1]);
738  velx = node[m1 - 1].velocity[s1][0] * cos(ang_u) - node[m1 - 1].velocity[s1][1] * sin(ang_u);
739  vely = node[m1 - 1].velocity[s1][0] * sin(ang_u) + node[m1 - 1].velocity[s1][1] * cos(ang_u);
740  angle = DefineAngle(velx, vely, u[0], u[1]);
741 
742  if (angle > eps) {
743  velx = node[m1 - 1].velocity[s1][0] * cos(ang_u) + node[m1 - 1].velocity[s1][1] * sin(ang_u);
744  vely = -1 * node[m1 - 1].velocity[s1][0] * sin(ang_u) + node[m1 - 1].velocity[s1][1] * cos(ang_u);
745  angle = DefineAngle(velx, vely, u[0], u[1]);
746  }
747 
748  if (angle < eps) {
749  node[m1 - 1].velocity[s1][0] = velx;
750  node[m1 - 1].velocity[s1][1] = vely;
751  }
752  } else {
753  if (((node[m1 - 1].typeN == 310) || (node[m1 - 1].typeN == 312)) && (node[m3 - 1].typeN != 310) && (node[m3 - 1].typeN != 312)) {
754  ang_v = DefineAngle(vk1[0], vk1[1], v[0], v[1]);
755  velx = node[m1 - 1].velocity[s1][0] * cos(ang_v) - node[m1 - 1].velocity[s1][1] * sin(ang_v);
756  vely = node[m1 - 1].velocity[s1][0] * sin(ang_v) + node[m1 - 1].velocity[s1][1] * cos(ang_v);
757  angle = DefineAngle(velx, vely, v[0], v[1]);
758 
759  if (angle > eps) {
760  velx = node[m1 - 1].velocity[s1][0] * cos(ang_v) + node[m1 - 1].velocity[s1][1] * sin(ang_v);
761  vely = -1 * node[m1 - 1].velocity[s1][0] * sin(ang_v) + node[m1 - 1].velocity[s1][1] * cos(ang_v);
762  angle = DefineAngle(velx, vely, v[0], v[1]);
763  }
764 
765  if (angle < eps) {
766  node[m1 - 1].velocity[s1][0] = velx;
767  node[m1 - 1].velocity[s1][1] = vely;
768  }
769  }
770  }
771  }
772 
773  return notcorner;
774 }
776 int XindexC(int nodenum, int ii) {
778  int xind = 0;
779 
780  if (node[nodenum - 1].fracture[0] == cell[ii].fracture) {
781  xind = 0;
782  }
783 
784  if (node[nodenum - 1].fracture[1] == cell[ii].fracture) {
785  xind = 3;
786  }
787 
788  return xind;
789 }
791 int YindexC(int nodenum, int ii) {
793  int yind = 0;
794 
795  if (node[nodenum - 1].fracture[0] == cell[ii].fracture) {
796  yind = 1;
797  }
798 
799  if (node[nodenum - 1].fracture[1] == cell[ii].fracture) {
800  yind = 4;
801  }
802 
803  return yind;
804 }
806 void HalfPolygonVelocity(int i, int k, int fractn, int indc, unsigned int fractj[max_neighb]) {
808  unsigned int subcell1f[max_neighb], subcell2f[max_neighb];
809  unsigned short int sc1f = 0, sc2f = 0, j;
810  double inters_v[2] = {0.0, 0.0}, pr, bx = 0, by = 0;
811  unsigned short int s1 = 0, ss = 0, kk, edge11 = 200, edge12 = 200;
812  double normxarea21[max_neighb - 1][2];
813  double normxarea22[max_neighb - 1][2];
814  double length = 1.0;
815  struct lb lbound;
816  int ic = 0;
817 
818  if (indc == 2) {
819  ic = 3;
820  }
821 
822  for (j = 0; j < k; j++) {
823  /* first, check if it is a boundary polygon */
824  ss = 0;
825 
826  for (kk = 0; kk < 4; kk++) {
827  if (node[i].cells[fractj[j]][kk] != 0) {
828  if (node[i].fracts[fractj[j]][kk] == fractn) {
829  ss = ss + 1;
830  }
831  }
832  }
833 
834  /* if this edge belongs to one cell in a fracture,
835  then it is a boundary edge and boundary angle will be defined */
836  if (ss == 1) {
837  if (edge11 == 200) {
838  edge11 = fractj[j];
839  } else {
840  edge12 = fractj[j];
841  }
842 
843  s1 = s1 + 1;
844  }
845 
846  /* define a vector of intersection edge (use it in vector cross product) */
847  if ((node[i].type[fractj[j]] == 2) || (node[i].type[fractj[j]] == 12)) {
848  inters_v[0] = node[i].coord_xy[ic] - node[node[i].indnodes[fractj[j]] - 1].coord_xy[ic];
849  inters_v[1] = node[i].coord_xy[ic + 1] - node[node[i].indnodes[fractj[j]] - 1].coord_xy[ic + 1];
850  }
851  } //loop j
852 
853  for (j = 0; j < k; j++) {
854  if (node[node[i].indnodes[fractj[j]] - 1].fracture[0] == fractn) {
855  bx = node[i].coord_xy[ic] - node[node[i].indnodes[fractj[j]] - 1].coord_xy[0];
856  by = node[i].coord_xy[ic + 1] - node[node[i].indnodes[fractj[j]] - 1].coord_xy[1];
857  length = sqrt(bx * bx + by * by);
858  }
859 
860  if (node[node[i].indnodes[fractj[j]] - 1].fracture[1] == fractn) {
861  bx = node[i].coord_xy[ic] - node[node[i].indnodes[fractj[j]] - 1].coord_xy[3];
862  by = node[i].coord_xy[ic + 1] - node[node[i].indnodes[fractj[j]] - 1].coord_xy[4];
863  length = sqrt(bx * bx + by * by);
864  }
865 
866  pr = bx * inters_v[1] - by * inters_v[0];
867 
868  /* sign of product defines the side of edges and devides polygon into
869  two subcells: opposite sides from intersection line */
870  if (pr <= 0.0) {
871  sc1f++;
872  subcell1f[sc1f - 1] = fractj[j];
873 
874  if (fehm == 1) {
875  normxarea21[sc1f - 1][0] = -1.*(bx * node[i].area[fractj[j]]);
876  normxarea21[sc1f - 1][1] = -1.*(by * node[i].area[fractj[j]]);
877  }
878 
879  if (pflotran == 1) {
880  normxarea21[sc1f - 1][0] = -1.*(bx * node[i].area[fractj[j]] / (length));
881  normxarea21[sc1f - 1][1] = -1.*(by * node[i].area[fractj[j]] / (length));
882  }
883 
884  if ((node[i].type[fractj[j]] != 2) && (node[i].type[fractj[j]] != 12)) {
885  kk = 0;
886 
887  /* defines indices of velocity vectors */
888  do {
889  if (cell[node[i].cells[fractj[j]][kk] - 1].node_ind[0] == i + 1) {
890  cell[node[i].cells[fractj[j]][kk] - 1].veloc_ind[0] = indc;
891  } else if (cell[node[i].cells[fractj[j]][kk] - 1].node_ind[1] == i + 1) {
892  cell[node[i].cells[fractj[j]][kk] - 1].veloc_ind[1] = indc;
893  } else {
894  cell[node[i].cells[fractj[j]][kk] - 1].veloc_ind[2] = indc;
895  }
896 
897  kk++;
898  } while (node[i].cells[fractj[j]][kk] != 0);
899  }
900  }
901 
902  if (pr >= 0.0) {
903  sc2f++;
904  subcell2f[sc2f - 1] = fractj[j];
905 
906  if (fehm == 1) {
907  normxarea22[sc2f - 1][0] = -1.*(bx * node[i].area[fractj[j]]);
908  normxarea22[sc2f - 1][1] = -1.*(by * node[i].area[fractj[j]]);
909  }
910 
911  if (pflotran == 1) {
912  normxarea22[sc2f - 1][0] = -1.*(bx * node[i].area[fractj[j]] / (length));
913  normxarea22[sc2f - 1][1] = -1.*(by * node[i].area[fractj[j]] / (length));
914  }
915 
916  if ((node[i].type[fractj[j]] != 2) && (node[i].type[fractj[j]] != 12)) {
917  kk = 0;
918 
919  /* defines indices of velocity vectors */
920  do {
921  if (cell[node[i].cells[fractj[j]][kk] - 1].node_ind[0] == i + 1) {
922  cell[node[i].cells[fractj[j]][kk] - 1].veloc_ind[0] = indc + 1;
923  } else if (cell[node[i].cells[fractj[j]][kk] - 1].node_ind[1] == i + 1) {
924  cell[node[i].cells[fractj[j]][kk] - 1].veloc_ind[1] = indc + 1;
925  } else {
926  cell[node[i].cells[fractj[j]][kk] - 1].veloc_ind[2] = indc + 1;
927  }
928 
929  kk++;
930  } while (node[i].cells[fractj[j]][kk] != 0);
931  }
932  }
933  }
934 
935  if ((s1 == 0) || (node[i].typeN == 302) || (node[i].typeN == 202) || (node[i].typeN == 312) || (node[i].typeN == 212)) {
936  if (sc1f != 0) {
937  VelocityInteriorNode (normxarea21, i, sc1f, subcell1f, indc);
938  }
939 
940  if (sc2f != 0) {
941  VelocityInteriorNode (normxarea22, i, sc2f, subcell2f, indc + 1);
942  }
943  } else {
944  if (s1 == 2) {
945  lbound = DefineBoundaryAngle (i, edge11, edge12, fractn, ic);
946 
947  if (sc1f != 0) {
948  VelocityExteriorNode (normxarea21, i, sc1f, subcell1f, lbound, indc );
949  }
950 
951  if (sc2f != 0) {
952  VelocityExteriorNode (normxarea22, i, sc2f, subcell2f, lbound, indc + 1 );
953  }
954  }
955  }
956 
957  return;
958 }
nnodes
unsigned int nnodes
Definition: main.c:29
pi
#define pi
Definition: FuncDef.h:2
cell
struct element * cell
Definition: main.c:42
porosity
double porosity
Definition: main.c:46
maindir
char maindir[125]
Definition: main.c:43
ncells
unsigned int ncells
Definition: main.c:30
timeunit
double timeunit
Definition: main.c:51
density
double density
Definition: main.c:47
pflotran
unsigned int pflotran
Definition: main.c:27
fehm
unsigned int fehm
Definition: main.c:28
max_neighb
unsigned int max_neighb
Definition: main.c:32
node
struct vertex * node
Definition: main.c:40
fracture
struct material * fracture
Definition: main.c:39
nfract
unsigned int nfract
Definition: main.c:31
HalfPolygonVelocity
void HalfPolygonVelocity(int i, int k, int fractn, int indc, unsigned int fractj[max_neighb])
Definition: VelocityReconstruction.c:806
YindexC
int YindexC(int nodenum, int ii)
Definition: VelocityReconstruction.c:791
DefineAngle
double DefineAngle(double u1, double u2, double v1, double v2)
Definition: VelocityReconstruction.c:347
VelocityInteriorNode
void VelocityInteriorNode(double normx_area[][2], int i, int number, unsigned int indj[max_neighb], int vi)
Definition: VelocityReconstruction.c:356
XindexC
int XindexC(int nodenum, int ii)
Definition: VelocityReconstruction.c:776
CornerVelocity
int CornerVelocity(int i, int m1, int m2, int m3, int s1, int s2, int s3)
Definition: VelocityReconstruction.c:658
BoundaryCells
void BoundaryCells()
Definition: VelocityReconstruction.c:447
DarcyVelocity
void DarcyVelocity()
Definition: VelocityReconstruction.c:29
MatrixProducts
struct matr MatrixProducts(double normxarea[][2], int number)
Definition: VelocityReconstruction.c:265
DefineBoundaryAngle
struct lb DefineBoundaryAngle(int i, unsigned int edge_1, unsigned int edge_2, int f1, int coorf)
Definition: VelocityReconstruction.c:304
VelocityExteriorNode
void VelocityExteriorNode(double norm_xarea[][2], int i, int number, unsigned int indj[max_neighb], struct lb lbound, int vi)
Definition: VelocityReconstruction.c:391
OutputVelocities
void OutputVelocities()
Definition: VelocityReconstruction.c:214
element::node_ind
unsigned int node_ind[3]
Definition: FuncDef.h:167
element::veloc_ind
unsigned int veloc_ind[3]
Definition: FuncDef.h:170
inpfile::filename
char filename[120]
Definition: InitialPartPositions.c:11
lb::angle
double angle
Definition: VelocityReconstruction.c:11
lb::norm
double norm[2]
Definition: VelocityReconstruction.c:12
lb::length_b
double length_b
Definition: VelocityReconstruction.c:10
matr::matrinvG
double matrinvG[2][40]
Definition: VelocityReconstruction.c:16
matr::matinv
double matinv[2][2]
Definition: VelocityReconstruction.c:17
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::indnodes
unsigned int * indnodes
Definition: FuncDef.h:136
vertex::velocity
double velocity[4][2]
Definition: FuncDef.h:130
vertex::flux
double * flux
Definition: FuncDef.h:148
vertex::fracture
unsigned int fracture[2]
Definition: FuncDef.h:112
vertex::typeN
unsigned int typeN
Definition: FuncDef.h:109
vertex::numneighb
unsigned int numneighb
Definition: FuncDef.h:127
vertex::type
unsigned int * type
Definition: FuncDef.h:145
vertex::area
double * area
Definition: FuncDef.h:151