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