OpenCAEPoro  v0.5.0
A simulator for multicomponent porous media flow
CornerGrid.cpp
Go to the documentation of this file.
1 
12 #include "CornerGrid.hpp"
13 
14 Point3D& Point3D::operator=(const Point3D& other)
15 {
16  x = other.x;
17  y = other.y;
18  z = other.z;
19  return *this;
20 }
21 
22 Point3D Point3D::operator+(const Point3D& other) const
23 {
24  return Point3D(x + other.x, y + other.y, z + other.z);
25 }
26 
27 Point3D Point3D::operator-(const Point3D& other) const
28 {
29  return Point3D(x - other.x, y - other.y, z - other.z);
30 }
31 
32 OCP_DBL Point3D::operator*(const Point3D& other) const
33 {
34  return x * other.x + y * other.y + z * other.z;
35 }
36 
37 Point3D& Point3D::operator+=(const Point3D& other)
38 {
39  x += other.x;
40  y += other.y;
41  z += other.z;
42  return *this;
43 }
44 
45 Point3D& Point3D::operator*=(const OCP_DBL& a)
46 {
47  x *= a;
48  y *= a;
49  z *= a;
50  return *this;
51 }
52 
53 Point3D& Point3D::operator/=(const OCP_DBL& a)
54 {
55  x /= a;
56  y /= a;
57  z /= a;
58  return *this;
59 }
60 
61 Point3D operator*(const Point3D& p, const OCP_DBL& a)
62 {
63  return Point3D(a * p.x, a * p.y, a * p.z);
64 }
65 
66 Point3D operator*(const OCP_DBL& a, const Point3D& p)
67 {
68  return Point3D(a * p.x, a * p.y, a * p.z);
69 }
70 
71 Point3D CrossProduct(const Point3D& p1, const Point3D& p2)
72 {
73  Point3D result;
74  result.x = p1.y * p2.z - p1.z * p2.y;
75  result.y = p1.z * p2.x - p1.x * p2.z;
76  result.z = p1.x * p2.y - p1.y * p2.x;
77  return result;
78 }
79 
80 Point3D Matrix3::operator*(const Point3D& v) const
81 {
82  Point3D result;
83  result.x = M[0][0] * v.x + M[0][1] * v.y + M[0][2] * v.z;
84  result.y = M[1][0] * v.x + M[1][1] * v.y + M[1][2] * v.z;
85  result.z = M[2][0] * v.x + M[2][1] * v.y + M[2][2] * v.z;
86  return result;
87 }
88 
89 OCP_DBL VolumHexahedron(const Hexahedron& h)
90 {
91  OCP_DBL result =
92  (h.p0.x * (h.p1.y * (-h.p2.z - h.p3.z + h.p4.z + h.p5.z) +
93  h.p2.y * (h.p1.z - h.p3.z) +
94  h.p3.y * (h.p1.z + h.p2.z - h.p4.z - h.p7.z) +
95  h.p4.y * (-h.p1.z + h.p3.z - h.p5.z + h.p7.z) +
96  h.p5.y * (-h.p1.z + h.p4.z) + h.p7.y * (h.p3.z - h.p4.z)) +
97  h.p1.x * (h.p0.y * (+h.p2.z + h.p3.z - h.p4.z - h.p5.z) +
98  h.p2.y * (-h.p0.z - h.p3.z + h.p5.z + h.p6.z) +
99  h.p3.y * (-h.p0.z + h.p2.z) + h.p4.y * (h.p0.z - h.p5.z) +
100  h.p5.y * (h.p0.z - h.p2.z + h.p4.z - h.p6.z) +
101  h.p6.y * (-h.p2.z + h.p5.z)) +
102  h.p2.x * (h.p0.y * (-h.p1.z + h.p3.z) +
103  h.p1.y * (h.p0.z + h.p3.z - h.p5.z - h.p6.z) +
104  h.p3.y * (-h.p0.z - h.p1.z + h.p6.z + h.p7.z) +
105  h.p5.y * (h.p1.z - h.p6.z) +
106  h.p6.y * (h.p1.z - h.p3.z + h.p5.z - h.p7.z) +
107  h.p7.y * (-h.p3.z + h.p6.z)) +
108  h.p3.x * (h.p0.y * (-h.p1.z - h.p2.z + h.p4.z + h.p7.z) +
109  h.p1.y * (h.p0.z - h.p2.z) +
110  h.p2.y * (h.p0.z + h.p1.z - h.p6.z - h.p7.z) +
111  h.p4.y * (-h.p0.z + h.p7.z) + h.p6.y * (h.p2.z - h.p7.z) +
112  h.p7.y * (-h.p0.z + h.p2.z - h.p4.z + h.p6.z)) +
113  h.p4.x * (h.p0.y * (h.p1.z - h.p3.z + h.p5.z - h.p7.z) +
114  h.p1.y * (-h.p0.z + h.p5.z) + h.p3.y * (h.p0.z - h.p7.z) +
115  h.p5.y * (-h.p0.z - h.p1.z + h.p6.z + h.p7.z) +
116  h.p6.y * (-h.p5.z + h.p7.z) +
117  h.p7.y * (h.p0.z + h.p3.z - h.p5.z - h.p6.z)) +
118  h.p5.x * (h.p0.y * (h.p1.z - h.p4.z) +
119  h.p1.y * (-h.p0.z + h.p2.z - h.p4.z + h.p6.z) +
120  h.p2.y * (-h.p1.z + h.p6.z) +
121  h.p4.y * (h.p0.z + h.p1.z - h.p6.z - h.p7.z) +
122  h.p6.y * (-h.p1.z - h.p2.z + h.p4.z + h.p7.z) +
123  h.p7.y * (h.p4.z - h.p6.z)) +
124  h.p6.x * (h.p1.y * (h.p2.z - h.p5.z) +
125  h.p2.y * (-h.p1.z + h.p3.z - h.p5.z + h.p7.z) +
126  h.p3.y * (-h.p2.z + h.p7.z) + h.p4.y * (h.p5.z - h.p7.z) +
127  h.p5.y * (h.p1.z + h.p2.z - h.p4.z - h.p7.z) +
128  h.p7.y * (-h.p2.z - h.p3.z + h.p4.z + h.p5.z)) +
129  h.p7.x * (h.p0.y * (-h.p3.z + h.p4.z) + h.p2.y * (h.p3.z - h.p6.z) +
130  h.p3.y * (h.p0.z - h.p2.z + h.p4.z - h.p6.z) +
131  h.p4.y * (-h.p0.z - h.p3.z + h.p5.z + h.p6.z) +
132  h.p5.y * (-h.p4.z + h.p6.z) +
133  h.p6.y * (h.p2.z + h.p3.z - h.p4.z - h.p5.z))) /
134  12;
135  return result;
136 }
137 
138 Point3D CenterHexahedron(const Hexahedron& h)
139 {
140  OCP_DBL r = 1.0 / 8.0;
141  Point3D result = r * (h.p0 + h.p1 + h.p2 + h.p3 + h.p4 + h.p5 + h.p6 + h.p7);
142  return result;
143 }
144 
145 Point3D VectorFace(const HexahedronFace& f)
146 {
147  Point3D p0, p1, p2, p3;
148  p0 = f.p2 - f.p1;
149  p1 = f.p0 - f.p1;
150  p2 = f.p0 - f.p3;
151  p3 = f.p2 - f.p3;
152  Point3D result = 0.5 * (CrossProduct(p0, p1) + CrossProduct(p2, p3));
153  return result;
154 }
155 
156 Point3D CenterFace(const HexahedronFace& f)
157 {
158  OCP_DBL r = 1.0 / 4.0;
159  Point3D result = r * (f.p0 + f.p1 + f.p2 + f.p3);
160  return result;
161 }
162 
163 Point2D CalCrossingPoint(const Point2D Line1[2], const Point2D Line2[2])
164 {
165  Point2D crosspoint;
166  //
167  // LOCALS
168  //
169  OCP_DBL a11, a12, a21, a22, b1, b2, detA, detX, detY;
170  //
171  // assume x = crosspoint.x
172  // y = crosspoint.y
173  // calculate x and y with equations in the following
174  //
175  // a11 a12 x b1
176  // [ ] ( ) = ( )
177  // a21 a22 y b2
178  //
179  a11 = Line1[1].y - Line1[0].y;
180  a12 = Line1[0].x - Line1[1].x;
181  a21 = Line2[1].y - Line2[0].y;
182  a22 = Line2[0].x - Line2[1].x;
183  b1 = a11 * Line1[0].x + a12 * Line1[0].y;
184  b2 = a21 * Line2[0].x + a22 * Line2[0].y;
185  detA = a11 * a22 - a12 * a21;
186 
187  if (fabs(detA) > SMALL_REAL) {
188  detX = b1 * a22 - b2 * a12;
189  detY = a11 * b2 - a21 * b1;
190  crosspoint.x = detX / detA;
191  crosspoint.y = detY / detA;
192  } else {
193  crosspoint = Line1[0];
194  }
195  return crosspoint;
196 }
197 
198 OCP_DBL CalAreaNotQuadr(const HexahedronFace& FACE1, const HexahedronFace& FACE2)
199 {
200  // Attention! Only for non quadrilateral!!! ---- Lishizhe
201  //
202  // This function calculate the common area of two quadrilaterals FACE1, FACE2.
203  //
204  // Order of points of Face follows
205  // 1 --- 0 0 --- 1
206  // | | or | |
207  // 2 --- 3 3 --- 2
208  // p0, p1 are upper, p2, p3 are lower
209  // y must be depth!!!
210  //
211  OCP_DBL CalAreaNotQuadr;
212  //
213  // LOCALS
214  //
215  USI iret;
216  Point2D crosspoint[4];
217  Point2D Line1[2], Line2[2];
218  HexahedronFace FACEtmp1, FACEtmp2;
219  Point3D area, point1, point2, point3;
220  //
221  CalAreaNotQuadr = 0;
222  iret = 0;
223  //
224  // the crossing relations of 4 lines:
225  // Line1 : point0 and point1 of face1
226  // Line2 : point2 and point3 of face1
227  // Line3 : point0 and point1 of face2
228  // Line4 : point2 and point3 of face2
229  //
230  // Line1 & Line3
231  //
232  Line1[0] = Point2D(FACE1.p0.x, FACE1.p0.y);
233  Line1[1] = Point2D(FACE1.p1.x, FACE1.p1.y);
234  Line2[0] = Point2D(FACE2.p0.x, FACE2.p0.y);
235  Line2[1] = Point2D(FACE2.p1.x, FACE2.p1.y);
236  crosspoint[0] = CalCrossingPoint(Line1, Line2);
237  if ((crosspoint[0].x - Line1[0].x) * (crosspoint[0].x - Line1[1].x) < 0)
238  iret = iret + 1;
239  //
240  // Line2 & Line3
241  //
242  Line1[0] = Point2D(FACE1.p2.x, FACE1.p2.y);
243  Line1[1] = Point2D(FACE1.p3.x, FACE1.p3.y);
244  Line2[0] = Point2D(FACE2.p0.x, FACE2.p0.y);
245  Line2[1] = Point2D(FACE2.p1.x, FACE2.p1.y);
246  crosspoint[1] = CalCrossingPoint(Line1, Line2);
247  if ((crosspoint[1].x - Line1[0].x) * (crosspoint[1].x - Line1[1].x) < 0)
248  iret = iret + 2;
249  //
250  // Line1 & Line4
251  //
252  Line1[0] = Point2D(FACE1.p0.x, FACE1.p0.y);
253  Line1[1] = Point2D(FACE1.p1.x, FACE1.p1.y);
254  Line2[0] = Point2D(FACE2.p2.x, FACE2.p2.y);
255  Line2[1] = Point2D(FACE2.p3.x, FACE2.p3.y);
256  crosspoint[2] = CalCrossingPoint(Line1, Line2);
257  if ((crosspoint[2].x - Line1[0].x) * (crosspoint[2].x - Line1[1].x) < 0)
258  iret = iret + 4;
259  //
260  // Line2 & Line4
261  //
262  Line1[0] = Point2D(FACE1.p2.x, FACE1.p2.y);
263  Line1[1] = Point2D(FACE1.p3.x, FACE1.p3.y);
264  Line2[0] = Point2D(FACE2.p2.x, FACE2.p2.y);
265  Line2[1] = Point2D(FACE2.p3.x, FACE2.p3.y);
266  crosspoint[3] = CalCrossingPoint(Line1, Line2);
267  if ((crosspoint[3].x - Line1[0].x) * (crosspoint[3].x - Line1[1].x) < 0)
268  iret = iret + 8;
269  //
270  // consider 12 cases of crossing relation combinations
271  //
272  switch (iret) {
273  case 1:
274  //
275  // Line1 & Line3 only
276  //
277  FACEtmp1.p1 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
278  FACEtmp2.p0 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
279 
280  if (FACE1.p0.y > FACE2.p0.y) {
281  FACEtmp1.p0 = FACE1.p0;
282  } else {
283  FACEtmp1.p0 = FACE2.p0;
284  }
285 
286  if (FACE1.p1.y > FACE2.p1.y) {
287  FACEtmp2.p1 = FACE1.p1;
288  } else {
289  FACEtmp2.p1 = FACE2.p1;
290  }
291 
292  if (FACE1.p3.y > FACE2.p3.y) {
293  FACEtmp1.p3 = FACE2.p3;
294  } else {
295  FACEtmp1.p3 = FACE1.p3;
296  }
297 
298  if (FACE1.p2.y > FACE2.p2.y) {
299  FACEtmp2.p2 = FACE2.p2;
300  } else {
301  FACEtmp2.p2 = FACE1.p2;
302  }
303 
304  FACEtmp1.p2 = Point3D(0.5 * (FACEtmp1.p3.x + FACEtmp2.p2.x),
305  0.5 * (FACEtmp1.p3.y + FACEtmp2.p2.y), 0);
306  FACEtmp2.p3 = FACEtmp1.p2;
307 
308  area = VectorFace(FACEtmp1);
309  CalAreaNotQuadr = fabs(area.z);
310  area = VectorFace(FACEtmp2);
311  CalAreaNotQuadr = CalAreaNotQuadr + fabs(area.z);
312  break;
313  case 2:
314  //
315  // Line2 & Line3 only
316  //
317  if (FACE1.p3.y > FACE2.p0.y) {
318  point1 = FACE1.p3;
319  point2 = FACE2.p0;
320  } else {
321  point1 = FACE1.p2;
322  point2 = FACE2.p1;
323  }
324  point3 = Point3D(crosspoint[1].x, crosspoint[1].y, 0);
325  area = CrossProduct(point1 - point3, point2 - point3);
326  CalAreaNotQuadr = fabs(area.z) * 0.5;
327  break;
328  case 3:
329  //
330  // Line1 & Line3
331  // Line2 & Line3
332  //
333  FACEtmp1.p0 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
334  FACEtmp1.p1 = Point3D(crosspoint[1].x, crosspoint[1].y, 0);
335  if (FACE1.p0.y < FACE2.p0.y) {
336  FACEtmp1.p2 = FACE1.p2;
337  FACEtmp1.p3 = FACE1.p1;
338  } else {
339  FACEtmp1.p2 = FACE1.p3;
340  FACEtmp1.p3 = FACE1.p0;
341  }
342  area = VectorFace(FACEtmp1);
343  CalAreaNotQuadr = fabs(area.z);
344  break;
345  case 4:
346  //
347  // Line1 & Line4 only
348  //
349  if (FACE1.p0.y < FACE2.p3.y) {
350  point1 = FACE1.p0;
351  point2 = FACE2.p3;
352  } else {
353  point1 = FACE1.p1;
354  point2 = FACE2.p2;
355  }
356  point3 = Point3D(crosspoint[2].x, crosspoint[2].y, 0);
357  area = CrossProduct(point1 - point3, point2 - point3);
358  CalAreaNotQuadr = fabs(area.z) * 0.5;
359  break;
360  case 5:
361  //
362  // Line1 & Line3
363  // Line1 & Line4
364  //
365  FACEtmp1.p0 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
366  FACEtmp1.p1 = Point3D(crosspoint[2].x, crosspoint[2].y, 0);
367  if (FACE2.p3.y > FACE1.p0.y) {
368  FACEtmp1.p2 = FACE2.p3;
369  FACEtmp1.p3 = FACE2.p0;
370  } else {
371  FACEtmp1.p2 = FACE2.p2;
372  FACEtmp1.p3 = FACE2.p1;
373  }
374  area = VectorFace(FACEtmp1);
375  CalAreaNotQuadr = fabs(area.z);
376  break;
377  case 8:
378  //
379  // Line2 & Line4 only
380  //
381  FACEtmp1.p2 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
382  FACEtmp2.p3 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
383 
384  if (FACE1.p0.y > FACE2.p0.y) {
385  FACEtmp1.p0 = FACE1.p0;
386  } else {
387  FACEtmp1.p0 = FACE2.p0;
388  }
389 
390  if (FACE1.p1.y > FACE2.p1.y) {
391  FACEtmp2.p1 = FACE1.p1;
392  } else {
393  FACEtmp2.p1 = FACE2.p1;
394  }
395 
396  if (FACE1.p3.y > FACE2.p3.y) {
397  FACEtmp1.p3 = FACE2.p3;
398  } else {
399  FACEtmp1.p3 = FACE1.p3;
400  }
401 
402  if (FACE1.p2.y > FACE2.p2.y) {
403  FACEtmp2.p2 = FACE2.p2;
404  } else {
405  FACEtmp2.p2 = FACE1.p2;
406  }
407 
408  FACEtmp1.p1 = Point3D(0.5 * (FACEtmp1.p0.x + FACEtmp2.p1.x),
409  0.5 * (FACEtmp1.p0.y + FACEtmp2.p1.y), 0);
410  FACEtmp2.p0 = FACEtmp1.p1;
411 
412  area = VectorFace(FACEtmp1);
413  CalAreaNotQuadr = fabs(area.z);
414  area = VectorFace(FACEtmp2);
415  CalAreaNotQuadr = CalAreaNotQuadr + fabs(area.z);
416  break;
417  case 9:
418  //
419  // Line1 & Line3
420  // Line2 & Line4
421  //
422  FACEtmp1.p1 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
423  FACEtmp2.p0 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
424  FACEtmp1.p2 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
425  FACEtmp2.p3 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
426 
427  if (FACE1.p0.y > FACE2.p0.y) {
428  FACEtmp1.p0 = FACE1.p0;
429  } else {
430  FACEtmp1.p0 = FACE2.p0;
431  }
432 
433  if (FACE1.p1.y > FACE2.p1.y) {
434  FACEtmp2.p1 = FACE1.p1;
435  } else {
436  FACEtmp2.p1 = FACE2.p1;
437  }
438 
439  if (FACE1.p3.y > FACE2.p3.y) {
440  FACEtmp1.p3 = FACE2.p3;
441  } else {
442  FACEtmp1.p3 = FACE1.p3;
443  }
444 
445  if (FACE1.p2.y > FACE2.p2.y) {
446  FACEtmp2.p2 = FACE2.p2;
447  } else {
448  FACEtmp2.p2 = FACE1.p2;
449  }
450 
451  area = VectorFace(FACEtmp1);
452  CalAreaNotQuadr = fabs(area.z);
453  area = VectorFace(FACEtmp2);
454  CalAreaNotQuadr = CalAreaNotQuadr + fabs(area.z);
455  break;
456  case 10:
457  //
458  // Line2 & Line3
459  // Line2 & Line4
460  //
461  FACEtmp1.p0 = Point3D(crosspoint[1].x, crosspoint[1].y, 0);
462  FACEtmp1.p1 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
463  if (FACE1.p2.y > FACE2.p1.y) {
464  FACEtmp1.p2 = FACE2.p1;
465  FACEtmp1.p3 = FACE2.p2;
466  } else {
467  FACEtmp1.p2 = FACE2.p3;
468  FACEtmp1.p3 = FACE2.p0;
469  }
470  area = VectorFace(FACEtmp1);
471  CalAreaNotQuadr = fabs(area.z);
472  break;
473  case 11:
474  //
475  // Line1 & Line3
476  // Line2 & Line3
477  // Line2 & Line4
478  //
479  FACEtmp1.p0 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
480  FACEtmp1.p2 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
481  FACEtmp1.p3 = Point3D(crosspoint[1].x, crosspoint[1].y, 0);
482  FACEtmp2.p3 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
483 
484  if (FACE1.p0.y < FACE2.p0.y) {
485  FACEtmp2.p1 = FACE1.p1;
486  FACEtmp2.p2 = FACE2.p2;
487  } else {
488  FACEtmp2.p1 = FACE1.p0;
489  FACEtmp2.p2 = FACE2.p3;
490  }
491 
492  FACEtmp1.p1 = Point3D(0.5 * (FACEtmp1.p0.x + FACEtmp2.p1.x),
493  0.5 * (FACEtmp1.p0.y + FACEtmp2.p1.y), 0);
494  FACEtmp2.p0 = FACEtmp1.p1;
495 
496  area = VectorFace(FACEtmp1);
497  CalAreaNotQuadr = fabs(area.z);
498  area = VectorFace(FACEtmp2);
499  CalAreaNotQuadr = CalAreaNotQuadr + fabs(area.z);
500  break;
501  case 12:
502  //
503  // Line1 & Line4
504  // Line2 & Line4
505  //
506  FACEtmp1.p0 = Point3D(crosspoint[2].x, crosspoint[2].y, 0);
507  FACEtmp1.p1 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
508  if (FACE1.p2.y > FACE2.p2.y) {
509  FACEtmp1.p2 = FACE1.p3;
510  FACEtmp1.p3 = FACE1.p0;
511  } else {
512  FACEtmp1.p2 = FACE1.p2;
513  FACEtmp1.p3 = FACE1.p1;
514  }
515  area = VectorFace(FACEtmp1);
516  CalAreaNotQuadr = fabs(area.z);
517  break;
518  case 13:
519  //
520  // Line1 & Line3
521  // Line1 & Line4
522  // Line2 & Line4
523  //
524  FACEtmp1.p2 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
525  FACEtmp2.p1 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
526  FACEtmp2.p2 = Point3D(crosspoint[2].x, crosspoint[2].y, 0);
527  FACEtmp2.p3 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
528 
529  if (FACE1.p2.y > FACE2.p2.y) {
530  FACEtmp1.p0 = FACE2.p0;
531  FACEtmp1.p3 = FACE1.p3;
532  } else {
533  FACEtmp1.p0 = FACE2.p1;
534  FACEtmp1.p3 = FACE1.p2;
535  }
536 
537  FACEtmp1.p1 = Point3D(0.5 * (FACEtmp1.p0.x + FACEtmp2.p1.x),
538  0.5 * (FACEtmp1.p0.y + FACEtmp2.p1.y), 0);
539  FACEtmp2.p0 = FACEtmp1.p1;
540 
541  area = VectorFace(FACEtmp1);
542  CalAreaNotQuadr = fabs(area.z);
543  area = VectorFace(FACEtmp2);
544  CalAreaNotQuadr = CalAreaNotQuadr + fabs(area.z);
545  break;
546  case 15:
547  //
548  // Line1 & Line3
549  // Line2 & Line3
550  // Line1 & Line4
551  // Line2 & Line4
552  //
553  FACEtmp1.p0 = Point3D(crosspoint[0].x, crosspoint[0].y, 0);
554  FACEtmp1.p1 = Point3D(crosspoint[2].x, crosspoint[2].y, 0);
555  FACEtmp1.p2 = Point3D(crosspoint[3].x, crosspoint[3].y, 0);
556  FACEtmp1.p3 = Point3D(crosspoint[1].x, crosspoint[1].y, 0);
557  area = VectorFace(FACEtmp1);
558  CalAreaNotQuadr = fabs(area.z);
559  break;
560  default:
561  CalAreaNotQuadr = 0;
562  break;
563  }
564  return CalAreaNotQuadr;
565 }
566 
567 void ConnGrid::Allocate(const USI& max_neighbor)
568 {
569  nConn = 0;
570  maxConn = max_neighbor;
571  halfConn.resize(maxConn);
572 }
573 
574 void ConnGrid::AddHalfConn(const OCP_USI& n,
575  const Point3D& area,
576  const Point3D& d,
577  const USI& direction,
578  const OCP_DBL& flag)
579 {
580  if (nConn >= maxConn) {
581  maxConn *= 2;
582  halfConn.resize(maxConn);
583  // get larger space
584  if (maxConn > MAX_NEIGHBOR) {
585  // OCP_ABORT("Too many Neighbors!");
586  }
587  }
588  halfConn[nConn].Ad_dd = area * d / (d * d) * flag;
589  halfConn[nConn].d = d;
590  halfConn[nConn].neigh = n;
591  halfConn[nConn].directionType = direction;
592  nConn++;
593 
594  // cout << n << " " << direction << " " << halfConn[nConn-1].Ad_dd << endl;
595 }
596 
597 void OCP_COORD::Allocate(const USI& Nx, const USI& Ny, const USI& Nz)
598 {
599  nx = Nx;
600  ny = Ny;
601  nz = Nz;
602  numGrid = nx * ny * nz;
603 
604  COORDDATA = new OCP_DBL**[3];
605  for (USI i = 0; i < 3; i++) {
606  COORDDATA[i] = new OCP_DBL*[2];
607  for (USI j = 0; j < 2; j++) {
608  COORDDATA[i][j] = new OCP_DBL[(nx + 1) * (ny + 1)];
609  }
610  }
611 
612  ZCORNDATA = new OCP_DBL***[nx];
613  for (USI i = 0; i < nx; i++) {
614  ZCORNDATA[i] = new OCP_DBL**[ny];
615  for (USI j = 0; j < ny; j++) {
616  ZCORNDATA[i][j] = new OCP_DBL*[nz];
617  for (USI k = 0; k < nz; k++) {
618  ZCORNDATA[i][j][k] = new OCP_DBL[8];
619  }
620  }
621  }
622 
623  cornerPoints = new Hexahedron**[nx];
624  for (USI i = 0; i < nx; i++) {
625  cornerPoints[i] = new Hexahedron*[ny];
626  for (USI j = 0; j < ny; j++) {
627  cornerPoints[i][j] = new Hexahedron[nz];
628  }
629  }
630 
631  v.resize(numGrid);
632  depth.resize(numGrid);
633  dx.resize(numGrid);
634  dy.resize(numGrid);
635  dz.resize(numGrid);
636  center.resize(numGrid);
637 }
638 
639 void OCP_COORD::InputData(const vector<OCP_DBL>& coord, const vector<OCP_DBL>& zcorn)
640 {
641  if (coord.empty() || !InputCOORDDATA(coord)) {
642  OCP_ABORT("ERROR COORD!");
643  }
644  if (zcorn.empty() || !InputZCORNDATA(zcorn)) {
645  OCP_ABORT("ERROR ZCORN!");
646  }
647 }
648 
649 OCP_BOOL OCP_COORD::InputCOORDDATA(const vector<OCP_DBL>& coord)
650 {
651  // See Eclipse -- COORD
652  OCP_BOOL flag = OCP_FALSE;
653  OCP_USI iter = 0;
654 
655  for (USI J = 0; J < ny + 1; J++) {
656  for (USI I = 0; I < nx + 1; I++) {
657  // top
658  for (USI i = 0; i < 3; i++) {
659  COORDDATA[i][0][J * (nx + 1) + I] = coord[iter];
660  iter++;
661  }
662  // bottom
663  for (USI i = 0; i < 3; i++) {
664  COORDDATA[i][1][J * (nx + 1) + I] = coord[iter];
665  iter++;
666  }
667  }
668  }
669 
670  flag = OCP_TRUE;
671  return flag;
672 }
673 
674 OCP_BOOL OCP_COORD::InputZCORNDATA(const vector<OCP_DBL>& zcorn)
675 {
676  // See Eclipse -- ZCORN
677  OCP_BOOL flag = OCP_FALSE;
678  OCP_USI iter = 0;
679 
680  for (USI K = 0; K < nz; K++) {
681  for (USI J = 0; J < ny; J++) {
682  for (USI I = 0; I < nx; I++) {
683  ZCORNDATA[I][J][K][0] = zcorn[iter];
684  iter++;
685  ZCORNDATA[I][J][K][1] = zcorn[iter];
686  iter++;
687  }
688  for (USI I = 0; I < nx; I++) {
689  ZCORNDATA[I][J][K][3] = zcorn[iter];
690  iter++;
691  ZCORNDATA[I][J][K][2] = zcorn[iter];
692  iter++;
693  }
694  }
695  for (USI J = 0; J < ny; J++) {
696  for (USI I = 0; I < nx; I++) {
697  ZCORNDATA[I][J][K][4] = zcorn[iter];
698  iter++;
699  ZCORNDATA[I][J][K][5] = zcorn[iter];
700  iter++;
701  }
702  for (USI I = 0; I < nx; I++) {
703  ZCORNDATA[I][J][K][7] = zcorn[iter];
704  iter++;
705  ZCORNDATA[I][J][K][6] = zcorn[iter];
706  iter++;
707  }
708  }
709  }
710 
711  flag = OCP_TRUE;
712  return flag;
713 }
714 
715 void OCP_COORD::SetAllFlags(const HexahedronFace& oFace, const HexahedronFace& Face)
716 {
717  tmpFace = Face;
718 
719  if (oFace.p0.z > Face.p0.z + TEENY) {
720  tmpFace.p0 = oFace.p0;
721  flagp0 = 1;
722  } else if (oFace.p0.z < Face.p0.z - TEENY)
723  flagp0 = -1;
724  else
725  flagp0 = 0;
726 
727  if (oFace.p1.z > Face.p1.z + TEENY)
728  flagp1 = 1;
729  else if (oFace.p1.z < Face.p1.z - TEENY) {
730  tmpFace.p1 = oFace.p1;
731  flagp1 = -1;
732  } else
733  flagp1 = 0;
734 
735  if (oFace.p2.z > Face.p2.z + TEENY)
736  flagp2 = 1;
737  else if (oFace.p2.z < Face.p2.z - TEENY) {
738  tmpFace.p2 = oFace.p2;
739  flagp2 = -1;
740  } else
741  flagp2 = 0;
742 
743  if (oFace.p3.z > Face.p3.z + TEENY) {
744  tmpFace.p3 = oFace.p3;
745  flagp3 = 1;
746  } else if (oFace.p3.z < Face.p3.z - TEENY)
747  flagp3 = -1;
748  else
749  flagp3 = 0;
750 
751  // check if interface is empty set
752  // check if interface is quadrilateral
753  // check if the one contains the other one
754 
755  if (((oFace.p1.z <= Face.p0.z) && (oFace.p2.z <= Face.p3.z)) ||
756  ((oFace.p0.z >= Face.p1.z) && (oFace.p3.z >= Face.p2.z))) {
757  flagJump = OCP_TRUE;
758  } else {
759  flagJump = OCP_FALSE;
760  if ((flagp0 * flagp3 >= 0) && (oFace.p0.z <= Face.p1.z) &&
761  (oFace.p3.z <= Face.p2.z) && (flagp1 * flagp2 >= 0) &&
762  (oFace.p1.z >= Face.p0.z) && (oFace.p2.z >= Face.p3.z)) {
763  flagQuad = OCP_TRUE;
764  } else {
765  flagQuad = OCP_FALSE;
766  }
767  }
768 }
769 
770 void OCP_COORD::SetupCornerPoints()
771 {
772  OCP_USI cindex, oindex; // current block index and the other block index
773  OCP_USI nxny = nx * ny;
774 
775  // allocate memoery for connections
776  vector<ConnGrid> blockconn(numGrid);
777  for (OCP_USI iloop = 0; iloop < numGrid; iloop++) {
778  blockconn[iloop].Allocate(10);
779  }
780 
781  // setup each block including coordinates of points, center, depth, and volume
782  OCP_DBL xtop, ytop, ztop, xbottom, ybottom, zbottom, xvalue, yvalue, zvalue;
783  for (USI k = 0; k < nz; k++) {
784  for (USI j = 0; j < ny; j++) {
785  for (USI i = 0; i < nx; i++) {
786  //
787  // corner point 0 and 4
788  //
789  xtop = COORDDATA[0][0][j * (nx + 1) + i];
790  ytop = COORDDATA[1][0][j * (nx + 1) + i];
791  ztop = COORDDATA[2][0][j * (nx + 1) + i];
792  xbottom = COORDDATA[0][1][j * (nx + 1) + i];
793  ybottom = COORDDATA[1][1][j * (nx + 1) + i];
794  zbottom = COORDDATA[2][1][j * (nx + 1) + i];
795 
796  zvalue = ZCORNDATA[i][j][k][0];
797  xvalue =
798  xbottom - (zbottom - zvalue) / (zbottom - ztop) * (xbottom - xtop);
799  yvalue =
800  ybottom - (zbottom - zvalue) / (zbottom - ztop) * (ybottom - ytop);
801  cornerPoints[i][j][k].p0 = Point3D(xvalue, yvalue, zvalue);
802 
803  zvalue = ZCORNDATA[i][j][k][4];
804  xvalue =
805  xbottom - (zbottom - zvalue) / (zbottom - ztop) * (xbottom - xtop);
806  yvalue =
807  ybottom - (zbottom - zvalue) / (zbottom - ztop) * (ybottom - ytop);
808  cornerPoints[i][j][k].p4 = Point3D(xvalue, yvalue, zvalue);
809  //
810  // corner point 1 and 5
811  //
812  xtop = COORDDATA[0][0][j * (nx + 1) + i + 1];
813  ytop = COORDDATA[1][0][j * (nx + 1) + i + 1];
814  ztop = COORDDATA[2][0][j * (nx + 1) + i + 1];
815  xbottom = COORDDATA[0][1][j * (nx + 1) + i + 1];
816  ybottom = COORDDATA[1][1][j * (nx + 1) + i + 1];
817  zbottom = COORDDATA[2][1][j * (nx + 1) + i + 1];
818 
819  zvalue = ZCORNDATA[i][j][k][1];
820  xvalue =
821  xbottom - (zbottom - zvalue) / (zbottom - ztop) * (xbottom - xtop);
822  yvalue =
823  ybottom - (zbottom - zvalue) / (zbottom - ztop) * (ybottom - ytop);
824  cornerPoints[i][j][k].p1 = Point3D(xvalue, yvalue, zvalue);
825 
826  zvalue = ZCORNDATA[i][j][k][5];
827  xvalue =
828  xbottom - (zbottom - zvalue) / (zbottom - ztop) * (xbottom - xtop);
829  yvalue =
830  ybottom - (zbottom - zvalue) / (zbottom - ztop) * (ybottom - ytop);
831  cornerPoints[i][j][k].p5 = Point3D(xvalue, yvalue, zvalue);
832  //
833  // corner point 2 and 6
834  //
835  xtop = COORDDATA[0][0][(j + 1) * (nx + 1) + i + 1];
836  ytop = COORDDATA[1][0][(j + 1) * (nx + 1) + i + 1];
837  ztop = COORDDATA[2][0][(j + 1) * (nx + 1) + i + 1];
838  xbottom = COORDDATA[0][1][(j + 1) * (nx + 1) + i + 1];
839  ybottom = COORDDATA[1][1][(j + 1) * (nx + 1) + i + 1];
840  zbottom = COORDDATA[2][1][(j + 1) * (nx + 1) + i + 1];
841 
842  zvalue = ZCORNDATA[i][j][k][2];
843  xvalue =
844  xbottom - (zbottom - zvalue) / (zbottom - ztop) * (xbottom - xtop);
845  yvalue =
846  ybottom - (zbottom - zvalue) / (zbottom - ztop) * (ybottom - ytop);
847  cornerPoints[i][j][k].p2 = Point3D(xvalue, yvalue, zvalue);
848 
849  zvalue = ZCORNDATA[i][j][k][6];
850  xvalue =
851  xbottom - (zbottom - zvalue) / (zbottom - ztop) * (xbottom - xtop);
852  yvalue =
853  ybottom - (zbottom - zvalue) / (zbottom - ztop) * (ybottom - ytop);
854  cornerPoints[i][j][k].p6 = Point3D(xvalue, yvalue, zvalue);
855  //
856  // corner point 3 and 7
857  //
858  xtop = COORDDATA[0][0][(j + 1) * (nx + 1) + i];
859  ytop = COORDDATA[1][0][(j + 1) * (nx + 1) + i];
860  ztop = COORDDATA[2][0][(j + 1) * (nx + 1) + i];
861  xbottom = COORDDATA[0][1][(j + 1) * (nx + 1) + i];
862  ybottom = COORDDATA[1][1][(j + 1) * (nx + 1) + i];
863  zbottom = COORDDATA[2][1][(j + 1) * (nx + 1) + i];
864 
865  zvalue = ZCORNDATA[i][j][k][3];
866  xvalue =
867  xbottom - (zbottom - zvalue) / (zbottom - ztop) * (xbottom - xtop);
868  yvalue =
869  ybottom - (zbottom - zvalue) / (zbottom - ztop) * (ybottom - ytop);
870  cornerPoints[i][j][k].p3 = Point3D(xvalue, yvalue, zvalue);
871 
872  zvalue = ZCORNDATA[i][j][k][7];
873  xvalue =
874  xbottom - (zbottom - zvalue) / (zbottom - ztop) * (xbottom - xtop);
875  yvalue =
876  ybottom - (zbottom - zvalue) / (zbottom - ztop) * (ybottom - ytop);
877  cornerPoints[i][j][k].p7 = Point3D(xvalue, yvalue, zvalue);
878 
879  // calculate volumes and pore volumes
880  cindex = k * nxny + j * nx + i;
881  //
882  // NOTE: if there are several points not well ordered, the calculated
883  // volume will be negative.
884  //
885  v[cindex] = VolumHexahedron(cornerPoints[i][j][k]); // NTG
886  v[cindex] = fabs(v[cindex]);
887  center[cindex] = CenterHexahedron(cornerPoints[i][j][k]);
888  depth[cindex] = center[cindex].z;
889  }
890  }
891  }
892 
893  // find neighbor and calculate transmissibility
894  OCP_USI num_conn = 0; // record the num of connection, a->b & b->a are both included
895  Point3D Pcenter, Pface, Pc2f; // center of Hexahedron
896  HexahedronFace Face, oFace; // current face, the other face
897  HexahedronFace FaceP, oFaceP; // Projection of Face and the other face
898  Point3D areaV; // area vector of interface
899  OCP_DBL areaP; // area of projection of interface
900  OCP_INT iznnc;
901  Point3D dxpoint, dypoint, dzpoint;
902 
903  // test
904  // cornerPoints[13][1][72].p0; cornerPoints[13][1][72].p1;
905  // cornerPoints[13][1][72].p2; cornerPoints[13][1][72].p3;
906  // cornerPoints[13][1][72].p4; cornerPoints[13][1][72].p5;
907  // cornerPoints[13][1][72].p6; cornerPoints[13][1][72].p7;
908 
909  // cornerPoints[13][2][74].p0; cornerPoints[13][2][74].p1;
910  // cornerPoints[13][2][74].p2; cornerPoints[13][2][74].p3;
911  // cornerPoints[13][2][74].p4; cornerPoints[13][2][74].p5;
912  // cornerPoints[13][2][74].p6; cornerPoints[13][2][74].p7;
913 
915  // Attention that The coordinate axis follows the right-hand rule ! //
917  //
918  // o----> x
919  // /|
920  // y z
921  // For a face, p0 and p3 are the points of upper edge of quadrilateral,
922  // p1 and p2 are the points of lower edge
923  // p0 ---- p3
924  // | |
925  // | |
926  // p1 ---- p2
927 
928  // Determine flagForward
929  if (COORDDATA[1][0][nx + 1] > COORDDATA[1][0][0])
930  flagForward = 1.0;
931  else
932  flagForward = -1.0;
933 
934  for (USI k = 0; k < nz; k++) {
935  for (USI j = 0; j < ny; j++) {
936  for (USI i = 0; i < nx; i++) {
937  // begin from each block
938  const Hexahedron& block = cornerPoints[i][j][k];
939  cindex = k * nxny + j * nx + i;
940  Pcenter = center[cindex];
941 
942  // cout << "============= " << cindex << " =============" << endl;
943  //
944  // (x-) direction
945  //
946 
947  Face.p0 = block.p0;
948  Face.p1 = block.p4;
949  Face.p2 = block.p7;
950  Face.p3 = block.p3;
951  Pface = CenterFace(Face);
952  Pc2f = Pface - Pcenter;
953  dxpoint = Pc2f;
954 
955  if (i == 0) {
956  // nothing to do
957  } else {
958 
959  const Hexahedron& leftblock = cornerPoints[i - 1][j][k];
960  oindex = k * nxny + j * nx + i - 1;
961 
962  oFace.p0 = leftblock.p1;
963  oFace.p1 = leftblock.p5;
964  oFace.p2 = leftblock.p6;
965  oFace.p3 = leftblock.p2;
966 
967  SetAllFlags(oFace, Face);
968 
969  // calculate the interface of two face
970  if (flagJump) {
971  // nothing to do
972  } else {
973  if (flagQuad) {
974  areaV = VectorFace(tmpFace);
975  } else {
976  FaceP.p0 = Point3D(Face.p3.y, Face.p3.z, 0);
977  FaceP.p1 = Point3D(Face.p0.y, Face.p0.z, 0);
978  FaceP.p2 = Point3D(Face.p1.y, Face.p1.z, 0);
979  FaceP.p3 = Point3D(Face.p2.y, Face.p2.z, 0);
980  oFaceP.p0 = Point3D(oFace.p3.y, oFace.p3.z, 0);
981  oFaceP.p1 = Point3D(oFace.p0.y, oFace.p0.z, 0);
982  oFaceP.p2 = Point3D(oFace.p1.y, oFace.p1.z, 0);
983  oFaceP.p3 = Point3D(oFace.p2.y, oFace.p2.z, 0);
984  areaP = CalAreaNotQuadr(FaceP, oFaceP);
985  // attention the direction of vector
986  areaV = VectorFace(Face);
987  // correct
988  if (fabs(areaV.x) < 1E-6) {
989  OCP_WARNING("x is too small");
990  } else {
991  areaV.y = areaV.y / fabs(areaV.x) * areaP;
992  areaV.z = areaV.z / fabs(areaV.x) * areaP;
993  areaV.x = OCP_SIGN(areaV.x) * areaP;
994  }
995  }
996  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 1,
997  flagForward);
998  num_conn++;
999  }
1000 
1001  // then find all NNC for current block
1002  // check if upNNC and downNNC exist
1003  if ((flagp0 > 0) || (flagp3 > 0))
1004  upNNC = OCP_TRUE;
1005  else
1006  upNNC = OCP_FALSE;
1007  if ((flagp1 < 0) || (flagp2 < 0))
1008  downNNC = OCP_TRUE;
1009  else
1010  downNNC = OCP_FALSE;
1011 
1012  iznnc = -1;
1013  while (upNNC) {
1014  // if (-iznnc > k) break;
1015  if (-iznnc - static_cast<OCP_INT>(k) > 0) break;
1016  // find object block
1017  const Hexahedron& leftblock = cornerPoints[i - 1][j][k + iznnc];
1018  oindex = (k + iznnc) * nxny + j * nx + i - 1;
1019  oFace.p0 = leftblock.p1;
1020  oFace.p1 = leftblock.p5;
1021  oFace.p2 = leftblock.p6;
1022  oFace.p3 = leftblock.p2;
1023 
1024  SetAllFlags(oFace, Face);
1025 
1026  // calculate the interface of two face
1027  if (flagJump) {
1028  // nothing to do
1029  } else {
1030  if (flagQuad) {
1031  areaV = VectorFace(tmpFace);
1032  } else {
1033  FaceP.p0 = Point3D(Face.p3.y, Face.p3.z, 0);
1034  FaceP.p1 = Point3D(Face.p0.y, Face.p0.z, 0);
1035  FaceP.p2 = Point3D(Face.p1.y, Face.p1.z, 0);
1036  FaceP.p3 = Point3D(Face.p2.y, Face.p2.z, 0);
1037  oFaceP.p0 = Point3D(oFace.p3.y, oFace.p3.z, 0);
1038  oFaceP.p1 = Point3D(oFace.p0.y, oFace.p0.z, 0);
1039  oFaceP.p2 = Point3D(oFace.p1.y, oFace.p1.z, 0);
1040  oFaceP.p3 = Point3D(oFace.p2.y, oFace.p2.z, 0);
1041  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1042  // attention the direction of vector
1043  areaV = VectorFace(Face);
1044  // correct
1045  if (fabs(areaV.x) < 1E-6) {
1046  OCP_WARNING("x is too small");
1047  } else {
1048  areaV.y = areaV.y / fabs(areaV.x) * areaP;
1049  areaV.z = areaV.z / fabs(areaV.x) * areaP;
1050  areaV.x = OCP_SIGN(areaV.x) * areaP;
1051  }
1052  }
1053  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 1,
1054  flagForward);
1055  num_conn++;
1056  }
1057  iznnc--;
1058  if ((flagp0 > 0) || (flagp3 > 0))
1059  upNNC = OCP_TRUE;
1060  else
1061  upNNC = OCP_FALSE;
1062  }
1063 
1064  iznnc = 1;
1065  while (downNNC) {
1066  if (k + iznnc > nz - 1) break;
1067  // find object block
1068  const Hexahedron& leftblock = cornerPoints[i - 1][j][k + iznnc];
1069  oindex = (k + iznnc) * nxny + j * nx + i - 1;
1070  oFace.p0 = leftblock.p1;
1071  oFace.p1 = leftblock.p5;
1072  oFace.p2 = leftblock.p6;
1073  oFace.p3 = leftblock.p2;
1074 
1075  SetAllFlags(oFace, Face);
1076 
1077  // calculate the interface of two face
1078  if (flagJump) {
1079  // nothing to do
1080  } else {
1081  if (flagQuad) {
1082  areaV = VectorFace(tmpFace);
1083  } else {
1084  FaceP.p0 = Point3D(Face.p3.y, Face.p3.z, 0);
1085  FaceP.p1 = Point3D(Face.p0.y, Face.p0.z, 0);
1086  FaceP.p2 = Point3D(Face.p1.y, Face.p1.z, 0);
1087  FaceP.p3 = Point3D(Face.p2.y, Face.p2.z, 0);
1088  oFaceP.p0 = Point3D(oFace.p3.y, oFace.p3.z, 0);
1089  oFaceP.p1 = Point3D(oFace.p0.y, oFace.p0.z, 0);
1090  oFaceP.p2 = Point3D(oFace.p1.y, oFace.p1.z, 0);
1091  oFaceP.p3 = Point3D(oFace.p2.y, oFace.p2.z, 0);
1092  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1093  // attention the direction of vector
1094  areaV = VectorFace(Face);
1095  // correct
1096  if (fabs(areaV.x) < 1E-6) {
1097  OCP_WARNING("x is too small");
1098  } else {
1099  areaV.y = areaV.y / fabs(areaV.x) * areaP;
1100  areaV.z = areaV.z / fabs(areaV.x) * areaP;
1101  areaV.x = OCP_SIGN(areaV.x) * areaP;
1102  }
1103  }
1104  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 1,
1105  flagForward);
1106  num_conn++;
1107  }
1108  iznnc++;
1109 
1110  if ((flagp1 < 0) || (flagp2 < 0))
1111  downNNC = OCP_TRUE;
1112  else
1113  downNNC = OCP_FALSE;
1114  }
1115  }
1116 
1117  //
1118  // (x+) direction
1119  //
1120  Face.p0 = block.p2;
1121  Face.p1 = block.p6;
1122  Face.p2 = block.p5;
1123  Face.p3 = block.p1;
1124  Pface = CenterFace(Face);
1125  Pc2f = Pface - Pcenter;
1126  dxpoint = Pc2f - dxpoint;
1127 
1128  if (i == nx - 1) {
1129  // nothing to do
1130  } else {
1131 
1132  const Hexahedron& rightblock = cornerPoints[i + 1][j][k];
1133  oindex = k * nxny + j * nx + i + 1;
1134 
1135  oFace.p0 = rightblock.p3;
1136  oFace.p1 = rightblock.p7;
1137  oFace.p2 = rightblock.p4;
1138  oFace.p3 = rightblock.p0;
1139 
1140  SetAllFlags(oFace, Face);
1141 
1142  // calculate the interface of two face
1143  if (flagJump) {
1144  // nothing to do
1145  } else {
1146  if (flagQuad) {
1147  areaV = VectorFace(tmpFace);
1148  } else {
1149  FaceP.p0 = Point3D(Face.p3.y, Face.p3.z, 0);
1150  FaceP.p1 = Point3D(Face.p0.y, Face.p0.z, 0);
1151  FaceP.p2 = Point3D(Face.p1.y, Face.p1.z, 0);
1152  FaceP.p3 = Point3D(Face.p2.y, Face.p2.z, 0);
1153  oFaceP.p0 = Point3D(oFace.p3.y, oFace.p3.z, 0);
1154  oFaceP.p1 = Point3D(oFace.p0.y, oFace.p0.z, 0);
1155  oFaceP.p2 = Point3D(oFace.p1.y, oFace.p1.z, 0);
1156  oFaceP.p3 = Point3D(oFace.p2.y, oFace.p2.z, 0);
1157  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1158  // attention the direction of vector
1159  areaV = VectorFace(Face);
1160  // correct
1161  if (fabs(areaV.x) < 1E-6) {
1162  OCP_WARNING("x is too small");
1163  } else {
1164  areaV.y = areaV.y / fabs(areaV.x) * areaP;
1165  areaV.z = areaV.z / fabs(areaV.x) * areaP;
1166  areaV.x = OCP_SIGN(areaV.x) * areaP;
1167  }
1168  }
1169  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 1,
1170  flagForward);
1171  num_conn++;
1172  }
1173 
1174  // then find all NNC for current block
1175  if ((flagp0 > 0) || (flagp3 > 0))
1176  upNNC = OCP_TRUE;
1177  else
1178  upNNC = OCP_FALSE;
1179  if ((flagp1 < 0) || (flagp2 < 0))
1180  downNNC = OCP_TRUE;
1181  else
1182  downNNC = OCP_FALSE;
1183 
1184  iznnc = -1;
1185  while (upNNC) {
1186  // if (-iznnc > k) break;
1187  if (-iznnc - static_cast<OCP_INT>(k) > 0) break;
1188  // find object block
1189  const Hexahedron& rightblock =
1190  cornerPoints[i + 1][j][k + iznnc];
1191  oindex = (k + iznnc) * nxny + j * nx + i + 1;
1192  oFace.p0 = rightblock.p3;
1193  oFace.p1 = rightblock.p7;
1194  oFace.p2 = rightblock.p4;
1195  oFace.p3 = rightblock.p0;
1196 
1197  SetAllFlags(oFace, Face);
1198 
1199  // calculate the interface of two face
1200  if (flagJump) {
1201  // nothing to do
1202  } else {
1203  if (flagQuad) {
1204  areaV = VectorFace(tmpFace);
1205  } else {
1206  FaceP.p0 = Point3D(Face.p3.y, Face.p3.z, 0);
1207  FaceP.p1 = Point3D(Face.p0.y, Face.p0.z, 0);
1208  FaceP.p2 = Point3D(Face.p1.y, Face.p1.z, 0);
1209  FaceP.p3 = Point3D(Face.p2.y, Face.p2.z, 0);
1210  oFaceP.p0 = Point3D(oFace.p3.y, oFace.p3.z, 0);
1211  oFaceP.p1 = Point3D(oFace.p0.y, oFace.p0.z, 0);
1212  oFaceP.p2 = Point3D(oFace.p1.y, oFace.p1.z, 0);
1213  oFaceP.p3 = Point3D(oFace.p2.y, oFace.p2.z, 0);
1214  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1215  // attention the direction of vector
1216  areaV = VectorFace(Face);
1217  // correct
1218  if (fabs(areaV.x) < 1E-6) {
1219  OCP_WARNING("x is too small");
1220  } else {
1221  areaV.y = areaV.y / fabs(areaV.x) * areaP;
1222  areaV.z = areaV.z / fabs(areaV.x) * areaP;
1223  areaV.x = OCP_SIGN(areaV.x) * areaP;
1224  }
1225  }
1226  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 1,
1227  flagForward);
1228  num_conn++;
1229  }
1230  iznnc--;
1231 
1232  if ((flagp0 > 0) || (flagp3 > 0))
1233  upNNC = OCP_TRUE;
1234  else
1235  upNNC = OCP_FALSE;
1236  }
1237 
1238  iznnc = 1;
1239  while (downNNC) {
1240  if (k + iznnc > nz - 1) break;
1241  // find object block
1242  const Hexahedron& rightblock =
1243  cornerPoints[i + 1][j][k + iznnc];
1244  oindex = (k + iznnc) * nxny + j * nx + i + 1;
1245  oFace.p0 = rightblock.p3;
1246  oFace.p1 = rightblock.p7;
1247  oFace.p2 = rightblock.p4;
1248  oFace.p3 = rightblock.p0;
1249 
1250  SetAllFlags(oFace, Face);
1251 
1252  // calculate the interface of two face
1253  if (flagJump) {
1254  // nothing to do
1255  } else {
1256  if (flagQuad) {
1257  areaV = VectorFace(tmpFace);
1258  } else {
1259  FaceP.p0 = Point3D(Face.p3.y, Face.p3.z, 0);
1260  FaceP.p1 = Point3D(Face.p0.y, Face.p0.z, 0);
1261  FaceP.p2 = Point3D(Face.p1.y, Face.p1.z, 0);
1262  FaceP.p3 = Point3D(Face.p2.y, Face.p2.z, 0);
1263  oFaceP.p0 = Point3D(oFace.p3.y, oFace.p3.z, 0);
1264  oFaceP.p1 = Point3D(oFace.p0.y, oFace.p0.z, 0);
1265  oFaceP.p2 = Point3D(oFace.p1.y, oFace.p1.z, 0);
1266  oFaceP.p3 = Point3D(oFace.p2.y, oFace.p2.z, 0);
1267  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1268  // attention the direction of vector
1269  areaV = VectorFace(Face);
1270  // correct
1271  if (fabs(areaV.x) < 1E-6) {
1272  OCP_WARNING("x is too small");
1273  } else {
1274  areaV.y = areaV.y / fabs(areaV.x) * areaP;
1275  areaV.z = areaV.z / fabs(areaV.x) * areaP;
1276  areaV.x = OCP_SIGN(areaV.x) * areaP;
1277  }
1278  }
1279  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 1,
1280  flagForward);
1281  num_conn++;
1282  }
1283  iznnc++;
1284 
1285  if ((flagp1 < 0) || (flagp2 < 0))
1286  downNNC = OCP_TRUE;
1287  else
1288  downNNC = OCP_FALSE;
1289  }
1290  }
1291 
1292  //
1293  // (y-) direction
1294  //
1295  Face.p0 = block.p1;
1296  Face.p1 = block.p5;
1297  Face.p2 = block.p4;
1298  Face.p3 = block.p0;
1299  Pface = CenterFace(Face);
1300  Pc2f = Pface - Pcenter;
1301  dypoint = Pc2f;
1302 
1303  if (j == 0) {
1304  // nothing to do
1305  } else {
1306 
1307  const Hexahedron& backblock = cornerPoints[i][j - 1][k];
1308  oindex = k * nxny + (j - 1) * nx + i;
1309 
1310  oFace.p0 = backblock.p2;
1311  oFace.p1 = backblock.p6;
1312  oFace.p2 = backblock.p7;
1313  oFace.p3 = backblock.p3;
1314 
1315  SetAllFlags(oFace, Face);
1316 
1317  // calculate the interface of two face
1318  if (flagJump) {
1319  // nothing to do
1320  } else {
1321  if (flagQuad) {
1322  areaV = VectorFace(tmpFace);
1323  } else {
1324  FaceP.p0 = Point3D(Face.p0.x, Face.p0.z, 0);
1325  FaceP.p1 = Point3D(Face.p3.x, Face.p3.z, 0);
1326  FaceP.p2 = Point3D(Face.p2.x, Face.p2.z, 0);
1327  FaceP.p3 = Point3D(Face.p1.x, Face.p1.z, 0);
1328  oFaceP.p0 = Point3D(oFace.p0.x, oFace.p0.z, 0);
1329  oFaceP.p1 = Point3D(oFace.p3.x, oFace.p3.z, 0);
1330  oFaceP.p2 = Point3D(oFace.p2.x, oFace.p2.z, 0);
1331  oFaceP.p3 = Point3D(oFace.p1.x, oFace.p1.z, 0);
1332  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1333  // attention the direction of vector
1334  areaV = VectorFace(Face);
1335  // correct
1336  if (fabs(areaV.y) < 1E-6) {
1337  OCP_WARNING("y is too small");
1338  } else {
1339  areaV.x = areaV.x / fabs(areaV.y) * areaP;
1340  areaV.z = areaV.z / fabs(areaV.y) * areaP;
1341  areaV.y = OCP_SIGN(areaV.y) * areaP;
1342  }
1343  }
1344  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 2,
1345  flagForward);
1346  num_conn++;
1347  }
1348 
1349  // then find all NNC for current block
1350  if ((flagp0 > 0) || (flagp3 > 0))
1351  upNNC = OCP_TRUE;
1352  else
1353  upNNC = OCP_FALSE;
1354  if ((flagp1 < 0) || (flagp2 < 0))
1355  downNNC = OCP_TRUE;
1356  else
1357  downNNC = OCP_FALSE;
1358 
1359  iznnc = -1;
1360  while (upNNC) {
1361  // if (-iznnc > k) break;
1362  if (-iznnc - static_cast<OCP_INT>(k) > 0) break;
1363  // find object block
1364  const Hexahedron& backblock = cornerPoints[i][j - 1][k + iznnc];
1365  oindex = (k + iznnc) * nxny + (j - 1) * nx + i;
1366  oFace.p0 = backblock.p2;
1367  oFace.p1 = backblock.p6;
1368  oFace.p2 = backblock.p7;
1369  oFace.p3 = backblock.p3;
1370 
1371  SetAllFlags(oFace, Face);
1372 
1373  // calculate the interface of two face
1374  if (flagJump) {
1375  // nothing to do
1376  } else {
1377  if (flagQuad) {
1378  areaV = VectorFace(tmpFace);
1379  } else {
1380  FaceP.p0 = Point3D(Face.p0.x, Face.p0.z, 0);
1381  FaceP.p1 = Point3D(Face.p3.x, Face.p3.z, 0);
1382  FaceP.p2 = Point3D(Face.p2.x, Face.p2.z, 0);
1383  FaceP.p3 = Point3D(Face.p1.x, Face.p1.z, 0);
1384  oFaceP.p0 = Point3D(oFace.p0.x, oFace.p0.z, 0);
1385  oFaceP.p1 = Point3D(oFace.p3.x, oFace.p3.z, 0);
1386  oFaceP.p2 = Point3D(oFace.p2.x, oFace.p2.z, 0);
1387  oFaceP.p3 = Point3D(oFace.p1.x, oFace.p1.z, 0);
1388  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1389  // attention the direction of vector
1390  areaV = VectorFace(Face);
1391  // correct
1392  if (fabs(areaV.y) < 1E-6) {
1393  OCP_WARNING("y is too small");
1394  } else {
1395  areaV.x = areaV.x / fabs(areaV.y) * areaP;
1396  areaV.z = areaV.z / fabs(areaV.y) * areaP;
1397  areaV.y = OCP_SIGN(areaV.y) * areaP;
1398  }
1399  }
1400  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 2,
1401  flagForward);
1402  num_conn++;
1403  }
1404  iznnc--;
1405 
1406  if ((flagp0 > 0) || (flagp3 > 0))
1407  upNNC = OCP_TRUE;
1408  else
1409  upNNC = OCP_FALSE;
1410  }
1411 
1412  iznnc = 1;
1413  while (downNNC) {
1414  if (k + iznnc > nz - 1) break;
1415  // find object block
1416  const Hexahedron& backblock = cornerPoints[i][j - 1][k + iznnc];
1417  oindex = (k + iznnc) * nxny + (j - 1) * nx + i;
1418  oFace.p0 = backblock.p2;
1419  oFace.p1 = backblock.p6;
1420  oFace.p2 = backblock.p7;
1421  oFace.p3 = backblock.p3;
1422 
1423  SetAllFlags(oFace, Face);
1424 
1425  // calculate the interface of two face
1426  if (flagJump) {
1427  // nothing to do
1428  } else {
1429  if (flagQuad) {
1430  areaV = VectorFace(tmpFace);
1431  } else {
1432  FaceP.p0 = Point3D(Face.p0.x, Face.p0.z, 0);
1433  FaceP.p1 = Point3D(Face.p3.x, Face.p3.z, 0);
1434  FaceP.p2 = Point3D(Face.p2.x, Face.p2.z, 0);
1435  FaceP.p3 = Point3D(Face.p1.x, Face.p1.z, 0);
1436  oFaceP.p0 = Point3D(oFace.p0.x, oFace.p0.z, 0);
1437  oFaceP.p1 = Point3D(oFace.p3.x, oFace.p3.z, 0);
1438  oFaceP.p2 = Point3D(oFace.p2.x, oFace.p2.z, 0);
1439  oFaceP.p3 = Point3D(oFace.p1.x, oFace.p1.z, 0);
1440  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1441  // attention the direction of vector
1442  areaV = VectorFace(Face);
1443  // correct
1444  if (fabs(areaV.y) < 1E-6) {
1445  OCP_WARNING("y is too small");
1446  } else {
1447  areaV.x = areaV.x / fabs(areaV.y) * areaP;
1448  areaV.z = areaV.z / fabs(areaV.y) * areaP;
1449  areaV.y = OCP_SIGN(areaV.y) * areaP;
1450  }
1451  }
1452  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 2,
1453  flagForward);
1454  num_conn++;
1455  }
1456  iznnc++;
1457 
1458  if ((flagp1 < 0) || (flagp2 < 0))
1459  downNNC = OCP_TRUE;
1460  else
1461  downNNC = OCP_FALSE;
1462  }
1463  }
1464 
1465  //
1466  // (y+) direction
1467  //
1468  Face.p0 = block.p3;
1469  Face.p1 = block.p7;
1470  Face.p2 = block.p6;
1471  Face.p3 = block.p2;
1472  Pface = CenterFace(Face);
1473  Pc2f = Pface - Pcenter;
1474  dypoint = Pc2f - dypoint;
1475 
1476  if (j == ny - 1) {
1477  // nothing to do
1478  } else {
1479 
1480  const Hexahedron& frontblock = cornerPoints[i][j + 1][k];
1481  oindex = k * nxny + (j + 1) * nx + i;
1482 
1483  oFace.p0 = frontblock.p0;
1484  oFace.p1 = frontblock.p4;
1485  oFace.p2 = frontblock.p5;
1486  oFace.p3 = frontblock.p1;
1487 
1488  SetAllFlags(oFace, Face);
1489 
1490  // calculate the interface of two face
1491  if (flagJump) {
1492  // nothing to do
1493  } else {
1494  if (flagQuad) {
1495  areaV = VectorFace(tmpFace);
1496  } else {
1497  FaceP.p0 = Point3D(Face.p0.x, Face.p0.z, 0);
1498  FaceP.p1 = Point3D(Face.p3.x, Face.p3.z, 0);
1499  FaceP.p2 = Point3D(Face.p2.x, Face.p2.z, 0);
1500  FaceP.p3 = Point3D(Face.p1.x, Face.p1.z, 0);
1501  oFaceP.p0 = Point3D(oFace.p0.x, oFace.p0.z, 0);
1502  oFaceP.p1 = Point3D(oFace.p3.x, oFace.p3.z, 0);
1503  oFaceP.p2 = Point3D(oFace.p2.x, oFace.p2.z, 0);
1504  oFaceP.p3 = Point3D(oFace.p1.x, oFace.p1.z, 0);
1505  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1506  // attention the direction of vector
1507  areaV = VectorFace(Face);
1508  // correct
1509  if (fabs(areaV.y) < 1E-6) {
1510  OCP_WARNING("y is too small");
1511  } else {
1512  areaV.x = areaV.x / fabs(areaV.y) * areaP;
1513  areaV.z = areaV.z / fabs(areaV.y) * areaP;
1514  areaV.y = OCP_SIGN(areaV.y) * areaP;
1515  }
1516  }
1517  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 2,
1518  flagForward);
1519  num_conn++;
1520  }
1521 
1522  // then find all NNC for current block
1523  if ((flagp0 > 0) || (flagp3 > 0))
1524  upNNC = OCP_TRUE;
1525  else
1526  upNNC = OCP_FALSE;
1527  if ((flagp1 < 0) || (flagp2 < 0))
1528  downNNC = OCP_TRUE;
1529  else
1530  downNNC = OCP_FALSE;
1531 
1532  iznnc = -1;
1533  while (upNNC) {
1534  // if (-iznnc > k) break;
1535  if (-iznnc - static_cast<OCP_INT>(k) > 0) break;
1536  // find object block
1537  const Hexahedron& frontblock =
1538  cornerPoints[i][j + 1][k + iznnc];
1539  oindex = (k + iznnc) * nxny + (j + 1) * nx + i;
1540  oFace.p0 = frontblock.p0;
1541  oFace.p1 = frontblock.p4;
1542  oFace.p2 = frontblock.p5;
1543  oFace.p3 = frontblock.p1;
1544 
1545  SetAllFlags(oFace, Face);
1546 
1547  // calculate the interface of two face
1548  if (flagJump) {
1549  // nothing to do
1550  } else {
1551  if (flagQuad) {
1552  areaV = VectorFace(tmpFace);
1553  } else {
1554  FaceP.p0 = Point3D(Face.p0.x, Face.p0.z, 0);
1555  FaceP.p1 = Point3D(Face.p3.x, Face.p3.z, 0);
1556  FaceP.p2 = Point3D(Face.p2.x, Face.p2.z, 0);
1557  FaceP.p3 = Point3D(Face.p1.x, Face.p1.z, 0);
1558  oFaceP.p0 = Point3D(oFace.p0.x, oFace.p0.z, 0);
1559  oFaceP.p1 = Point3D(oFace.p3.x, oFace.p3.z, 0);
1560  oFaceP.p2 = Point3D(oFace.p2.x, oFace.p2.z, 0);
1561  oFaceP.p3 = Point3D(oFace.p1.x, oFace.p1.z, 0);
1562  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1563  // attention the direction of vector
1564  areaV = VectorFace(Face);
1565  // correct
1566  if (fabs(areaV.y) < 1E-6) {
1567  OCP_WARNING("y is too small");
1568  } else {
1569  areaV.x = areaV.x / fabs(areaV.y) * areaP;
1570  areaV.z = areaV.z / fabs(areaV.y) * areaP;
1571  areaV.y = OCP_SIGN(areaV.y) * areaP;
1572  }
1573  }
1574  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 2,
1575  flagForward);
1576  num_conn++;
1577  }
1578  iznnc--;
1579 
1580  if ((flagp0 > 0) || (flagp3 > 0))
1581  upNNC = OCP_TRUE;
1582  else
1583  upNNC = OCP_FALSE;
1584  }
1585 
1586  iznnc = 1;
1587  while (downNNC) {
1588  if (k + iznnc > nz - 1) break;
1589  // find object block
1590  const Hexahedron& frontblock =
1591  cornerPoints[i][j + 1][k + iznnc];
1592  oindex = (k + iznnc) * nxny + (j + 1) * nx + i;
1593  oFace.p0 = frontblock.p0;
1594  oFace.p1 = frontblock.p4;
1595  oFace.p2 = frontblock.p5;
1596  oFace.p3 = frontblock.p1;
1597 
1598  SetAllFlags(oFace, Face);
1599 
1600  // calculate the interface of two face
1601  if (flagJump) {
1602  // nothing to do
1603  } else {
1604  if (flagQuad) {
1605  areaV = VectorFace(tmpFace);
1606  } else {
1607  FaceP.p0 = Point3D(Face.p0.x, Face.p0.z, 0);
1608  FaceP.p1 = Point3D(Face.p3.x, Face.p3.z, 0);
1609  FaceP.p2 = Point3D(Face.p2.x, Face.p2.z, 0);
1610  FaceP.p3 = Point3D(Face.p1.x, Face.p1.z, 0);
1611  oFaceP.p0 = Point3D(oFace.p0.x, oFace.p0.z, 0);
1612  oFaceP.p1 = Point3D(oFace.p3.x, oFace.p3.z, 0);
1613  oFaceP.p2 = Point3D(oFace.p2.x, oFace.p2.z, 0);
1614  oFaceP.p3 = Point3D(oFace.p1.x, oFace.p1.z, 0);
1615  areaP = CalAreaNotQuadr(FaceP, oFaceP);
1616  // attention the direction of vector
1617  areaV = VectorFace(Face);
1618  // correct
1619  if (fabs(areaV.y) < 1E-6) {
1620  OCP_WARNING("y is too small");
1621  } else {
1622  areaV.x = areaV.x / fabs(areaV.y) * areaP;
1623  areaV.z = areaV.z / fabs(areaV.y) * areaP;
1624  areaV.y = OCP_SIGN(areaV.y) * areaP;
1625  }
1626  }
1627  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 2,
1628  flagForward);
1629  num_conn++;
1630  }
1631  iznnc++;
1632 
1633  if ((flagp1 < 0) || (flagp2 < 0))
1634  downNNC = OCP_TRUE;
1635  else
1636  downNNC = OCP_FALSE;
1637  }
1638  }
1639 
1640  //
1641  // (z-) direction
1642  //
1643  Face.p0 = block.p0;
1644  Face.p1 = block.p3;
1645  Face.p2 = block.p2;
1646  Face.p3 = block.p1;
1647  Pface = CenterFace(Face);
1648  Pc2f = Pface - Pcenter;
1649  dzpoint = Pc2f;
1650  if (k == 0) {
1651  // nothing to do
1652  } else {
1653  // upblock
1654  oindex = (k - 1) * nxny + j * nx + i;
1655 
1656  tmpFace = Face;
1657  areaV = VectorFace(tmpFace);
1658  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 3, flagForward);
1659  num_conn++;
1660  }
1661 
1662  //
1663  // (z+) direction
1664  //
1665  Face.p0 = block.p5;
1666  Face.p1 = block.p6;
1667  Face.p2 = block.p7;
1668  Face.p3 = block.p4;
1669  Pface = CenterFace(Face);
1670  Pc2f = Pface - Pcenter;
1671  dzpoint = Pc2f - dzpoint;
1672 
1673  if (k == nz - 1) {
1674  // nothing to do
1675  } else {
1676  // downblock
1677  oindex = (k + 1) * nxny + j * nx + i;
1678 
1679  tmpFace = Face;
1680  areaV = VectorFace(tmpFace);
1681  blockconn[cindex].AddHalfConn(oindex, areaV, Pc2f, 3, flagForward);
1682  num_conn++;
1683  }
1684 
1685  // calculate dx,dy,dz
1686  dx[cindex] = sqrt(dxpoint.x * dxpoint.x + dxpoint.y * dxpoint.y +
1687  dxpoint.z * dxpoint.z);
1688  dy[cindex] = sqrt(dypoint.x * dypoint.x + dypoint.y * dypoint.y +
1689  dypoint.z * dypoint.z);
1690  dz[cindex] = sqrt(dzpoint.x * dzpoint.x + dzpoint.y * dzpoint.y +
1691  dzpoint.z * dzpoint.z);
1692 
1693  OCP_ASSERT(isfinite(dx[cindex]), "Wrong dx!");
1694  OCP_ASSERT(isfinite(dy[cindex]), "Wrong dy!");
1695  OCP_ASSERT(isfinite(dz[cindex]), "Wrong dz!");
1696  }
1697  }
1698  }
1699 
1700  OCP_ASSERT(num_conn % 2 == 0, "Wrong Conn!");
1701  numConnMax = num_conn / 2;
1702  connect.resize(numConnMax);
1703  //
1704  // calculate the x,y,z direction transmissibilities of each block and save them
1705  //
1706  // make the connections
1707  OCP_USI iter_conn = 0;
1708  for (OCP_USI n = 0; n < numGrid; n++) {
1709  for (USI j = 0; j < blockconn[n].nConn; j++) {
1710  OCP_USI nn = blockconn[n].halfConn[j].neigh;
1711  if (nn < n) continue;
1712  USI jj;
1713  for (jj = 0; jj < blockconn[nn].nConn; jj++) {
1714  if (blockconn[nn].halfConn[jj].neigh == n) {
1715  break;
1716  }
1717  }
1718  if (jj == blockconn[nn].nConn) {
1719  continue;
1720  }
1721  if (blockconn[n].halfConn[j].Ad_dd <= 0 ||
1722  blockconn[nn].halfConn[jj].Ad_dd <= 0) {
1723  // OCP_FALSE connection
1724  continue;
1725  }
1726 
1727  //
1728  // now, blockconn[n].halfConn[j]
1729  // blockconn[nn].halfConn[jj]
1730  // are a pair of connections
1731  connect[iter_conn].begin = n;
1732  connect[iter_conn].Ad_dd_begin = blockconn[n].halfConn[j].Ad_dd;
1733  connect[iter_conn].end = nn;
1734  connect[iter_conn].Ad_dd_end = blockconn[nn].halfConn[jj].Ad_dd;
1735  connect[iter_conn].directionType = blockconn[n].halfConn[j].directionType;
1736  iter_conn++;
1737  }
1738  }
1739  numConn = iter_conn;
1740 }
1741 
1742 /*----------------------------------------------------------------------------*/
1743 /* Brief Change History of This File */
1744 /*----------------------------------------------------------------------------*/
1745 /* Author Date Actions */
1746 /*----------------------------------------------------------------------------*/
1747 /* Shizhe Li Nov/19/2021 Create file */
1748 /*----------------------------------------------------------------------------*/
CalAreaNotQuadr
OCP_DBL CalAreaNotQuadr(const HexahedronFace &FACE1, const HexahedronFace &FACE2)
???
Definition: CornerGrid.cpp:198
VectorFace
Point3D VectorFace(const HexahedronFace &f)
Find the normal vector of a face.
Definition: CornerGrid.cpp:145
CrossProduct
Point3D CrossProduct(const Point3D &p1, const Point3D &p2)
Cross product.
Definition: CornerGrid.cpp:71
CalCrossingPoint
Point2D CalCrossingPoint(const Point2D Line1[2], const Point2D Line2[2])
???
Definition: CornerGrid.cpp:163
CenterHexahedron
Point3D CenterHexahedron(const Hexahedron &h)
Find the center of a hexahedron.
Definition: CornerGrid.cpp:138
CenterFace
Point3D CenterFace(const HexahedronFace &f)
Find the center of a face.
Definition: CornerGrid.cpp:156
operator*
Point3D operator*(const Point3D &p, const OCP_DBL &a)
Point * a.
Definition: CornerGrid.cpp:61
VolumHexahedron
OCP_DBL VolumHexahedron(const Hexahedron &h)
Get the volume of a hexahedron.
Definition: CornerGrid.cpp:89
CornerGrid.hpp
Declaration of classes related to the corner grid.
TEENY
const OCP_DBL TEENY
Used for checking distance b/w center to face.
Definition: CornerGrid.hpp:27
SMALL_REAL
const OCP_DBL SMALL_REAL
Used for checking determinate of a small matrix.
Definition: CornerGrid.hpp:26
MAX_NEIGHBOR
const USI MAX_NEIGHBOR
Max number of neighbors allowed.
Definition: CornerGrid.hpp:29
USI
unsigned int USI
Generic unsigned integer.
Definition: OCPConst.hpp:23
OCP_DBL
double OCP_DBL
Double precision.
Definition: OCPConst.hpp:27
OCP_USI
unsigned int OCP_USI
Long unsigned integer.
Definition: OCPConst.hpp:25
OCP_INT
int OCP_INT
Long integer.
Definition: OCPConst.hpp:26
OCP_BOOL
unsigned int OCP_BOOL
OCP_BOOL in OCP.
Definition: OCPConst.hpp:29
OCP_WARNING
#define OCP_WARNING(msg)
Log warning messages.
Definition: UtilError.hpp:39
OCP_ASSERT
#define OCP_ASSERT(cond, msg)
Assert condition and log user messages in DEBUG mode.
Definition: UtilError.hpp:58
OCP_ABORT
#define OCP_ABORT(msg)
Abort if critical error happens.
Definition: UtilError.hpp:47
HexahedronFace
A face of a hexahedron cell.
Definition: CornerGrid.hpp:88
Hexahedron
A hexahedron cell.
Definition: CornerGrid.hpp:81
Point2D
A point in 2D.
Definition: CornerGrid.hpp:33
Point3D
A point in 3D.
Definition: CornerGrid.hpp:47
Point3D::operator-
Point3D operator-(const Point3D &other) const
Subtraction.
Definition: CornerGrid.cpp:27
Point3D::operator+
Point3D operator+(const Point3D &other) const
Addition.
Definition: CornerGrid.cpp:22
Point3D::operator=
Point3D & operator=(const Point3D &other)
equal
Definition: CornerGrid.cpp:14
Point3D::operator*
OCP_DBL operator*(const Point3D &other) const
Multiplication.
Definition: CornerGrid.cpp:32
Generated on Tue Feb 7 2023 01:49:50 for OpenCAEPoro by doxygen 1.9.1