PDRutilsIntersect.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2016 Shell Research Ltd.
9 Copyright (C) 2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "PDRsetFields.H"
30#include "PDRutilsInternal.H"
32
33#ifndef FULLDEBUG
34#ifndef NDEBUG
35#define NDEBUG
36#endif
37#endif
38#include <cassert>
39
40using namespace Foam::constant;
41
42// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
43
44namespace Foam
45{
46
47// Calculate the area of the sector of a circle whose ends are at
48// (dxa, dya) and (dxb, dyb) relative to the centre. radsqu is radius
49// squared.
50//
51// (We trust that this is consistent with the other parameters..)
52inline static scalar sector
53(
54 scalar dxa, scalar dya,
55 scalar dxb, scalar dyb
56)
57{
58 scalar angle = (::atan2(dyb, dxb) - ::atan2(dya, dxa));
59
60 if (angle < -1.0E-10)
61 {
62 angle += mathematical::twoPi;
63 }
64
65 return angle;
66}
67
68} // End namespace Foam
69
70
71// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
72
74(
75 double xc, double yc, double rad,
76 double x1, double x2,
77 double y1, double y2,
78 scalar* perim_p,
79 scalar* x_proj_edge_p, scalar* y_proj_edge_p,
80 scalar* x_overlap_p, scalar* y_overlap_p
81)
82{
83 double angle, area, del;
84 double x_int[6][2], y_int[6][2]; // Coordinates of intersections between the circle and sides of the rectangle.
85 double x_arc[6][2], y_arc[6][2]; // Coordinates of end orc (within cell) for each quadrant
86 double dx[6], dy[6];
87
88 double x_olap_min = GREAT;
89 double x_olap_max = -GREAT;
90 double y_olap_min = GREAT;
91 double y_olap_max = -GREAT;
92 int n_vert, n_oppv;
93 int no_intersection;
94
95 const double dx1 = (x1 - xc);
96 const double dx2 = (x2 - xc);
97 const double dy1 = (y1 - yc);
98 const double dy2 = (y2 - yc);
99
100 const double dx1squ = dx1 * dx1;
101 const double dx2squ = dx2 * dx2;
102 const double dy1squ = dy1 * dy1;
103 const double dy2squ = dy2 * dy2;
104 const double radsqu = rad * rad;
105
106
107 /* Going clockwise from (x1, y1), vertices are labelled 1,2,3,4, with 0 the same as 4
108 and 5 the same as 1 (so that we can find the vertices on either side of any of 1 to 4).*/
109
110 dx[1] = dx1; dy[1] = dy1;
111 dx[2] = dx1; dy[2] = dy2;
112 dx[3] = dx2; dy[3] = dy2;
113 dx[4] = dx2; dy[4] = dy1;
114 dx[0] = dx2; dy[0] = dy1;
115 dx[5] = dx1; dy[5] = dy1;
116
117 // The positions of the ends of the arcs, if these points are
118 // inside the cell, they will be changed, if necessary, below.
119
120 x_arc[2][0] = x_arc[1][1] = -rad; y_arc[2][0] = y_arc[1][1] = 0.0;
121 x_arc[3][0] = x_arc[2][1] = 0.0; y_arc[3][0] = y_arc[2][1] = rad;
122 x_arc[4][0] = x_arc[3][1] = rad; y_arc[4][0] = y_arc[3][1] = 0.0;
123 x_arc[1][0] = x_arc[4][1] = 0.0; y_arc[1][0] = y_arc[4][1] = -rad;
124
125 // We catch arcs that are entirely inside the rectangle
126 // Note: this is wrong for a circle completely outside, but that
127 // will be dealt with separately
128
129 int arc_in[6] = { /* zero-initialied */ };
130 arc_in[1] = (dx1 < -rad && dy1 < -rad) ? 1 : 0;
131 arc_in[2] = (dx1 < -rad && dy2 > rad) ? 1 : 0;
132 arc_in[3] = (dx2 > rad && dy2 > rad) ? 1 : 0;
133 arc_in[4] = (dx2 > rad && dy1 < -rad) ? 1 : 0;
134
135 // Work out which vertices are in the circle
136
137 int vert_in[6];
138 vert_in[1] = (dx1squ + dy1squ <= radsqu);
139 vert_in[2] = (dx1squ + dy2squ <= radsqu);
140 vert_in[3] = (dx2squ + dy2squ <= radsqu);
141 vert_in[4] = (dx2squ + dy1squ <= radsqu);
142 vert_in[0] = vert_in[4];
143 vert_in[5] = vert_in[1];
144
145 int n_in = 0;
146 if (vert_in[1]) ++n_in;
147 if (vert_in[2]) ++n_in;
148 if (vert_in[3]) ++n_in;
149 if (vert_in[4]) ++n_in;
150
151
152 /* We now calculate the points of intersection of the circle with, successively,
153 x=x1, y=y2, x=x2. y=y1.
154
155 Where there are two intersections with one side, need to be careful about
156 the order of the two points (i.e. clockwise round the rectangle) so that
157 later on we get the right sector (short or long way round the circumference) */
158
159 int n_int[6] = { /* zero-initialied */ };
160 n_int[1] = 0;
161 if ( dx1squ <= radsqu)
162 {
163 del = std::sqrt( radsqu - dx1squ);
164 if ( ( ( -del ) <= dy2 ) && ( del >= dy1 ) )
165 {
166 x_int[1][0] = x_int[1][1] = dx1;
167 if ( (-del ) > dy1 )
168 {
169 y_int[1][0] = -del;
170 n_int[1]++;
171 // This intersection will be an end of the 3rd- or 4th-quadrant arc
172 if ( dx1 > 0.0 ) { x_arc[4][1] = dx1; y_arc[4][1] = -del; arc_in[4] = 1; }
173 else { x_arc[1][1] = dx1; y_arc[1][1] = -del; arc_in[1] = 1; }
174 }
175 if ( ( del ) < dy2 )
176 {
177 y_int[1][n_int[1]] = del;
178 n_int[1]++;
179 if ( dx1 > 0.0 ) { x_arc[3][0] = dx1; y_arc[3][0] = del; arc_in[3] = 1; }
180 else { x_arc[2][0] = dx1; y_arc[2][0] = del; arc_in[2] = 1; }
181 }
182 }
183 }
184
185 n_int[2] = 0;
186 if ( dy2squ <= radsqu)
187 {
188 del = std::sqrt( radsqu - dy2squ);
189 if ( ( ( -del ) <= dx2 ) && ( del >= dx1 ) )
190 {
191 y_int[2][0] = y_int[2][1] = dy2;
192 if ( (-del ) > dx1 )
193 {
194 x_int[2][0] = -del;
195 n_int[2]++;
196 if ( dy2 > 0.0 ) { x_arc[2][1] = -del; y_arc[2][1] = dy2; arc_in[2] = 1; }
197 else { x_arc[1][1] = -del; y_arc[1][1] = dy2; arc_in[1] = 1; }
198 }
199 if ( ( del ) < dx2 )
200 {
201 x_int[2][n_int[2]] = del;
202 n_int[2]++;
203 if ( dy2 > 0.0 ) { x_arc[3][0] = del; y_arc[3][0] = dy2; arc_in[3] = 1; }
204 else { x_arc[4][0] = del; y_arc[4][0] = dy2; arc_in[4] = 1; }
205 }
206 }
207 }
208
209 n_int[3] = 0;
210 if ( dx2squ <= radsqu)
211 {
212 del = std::sqrt( radsqu - dx2squ);
213 if ( ( ( -del ) <= dy2 ) && ( del >= dy1 ) )
214 {
215 x_int[3][0] = x_int[3][1] = dx2;
216 if ( ( del ) < dy2 )
217 {
218 y_int[3][0] = del;
219 n_int[3]++;
220 if ( dx2 > 0.0 ) { x_arc[3][1] = dx2; y_arc[3][1] = del; arc_in[3] = 1; }
221 else { x_arc[2][1] = dx2; y_arc[2][1] = del; arc_in[2] = 1; }
222 }
223 if ( (-del ) > dy1 )
224 {
225 y_int[3][n_int[3]] = -del;
226 n_int[3]++;
227 if ( dx2 > 0.0 ) { x_arc[4][0] = dx2; y_arc[4][0] = -del; arc_in[4] = 1; }
228 else { x_arc[1][0] = dx2; y_arc[1][0] = -del; arc_in[1] = 1; }
229 }
230 }
231 }
232
233 n_int[4] = 0;
234 if ( dy1squ <= radsqu)
235 {
236 del = std::sqrt( radsqu - dy1squ);
237 if ( ( ( -del ) <= dx2 ) && ( del >= dx1 ) )
238 {
239 y_int[4][0] = y_int[4][1] = dy1;
240 if ( ( del ) < dx2 )
241 {
242 x_int[4][0] = del;
243 n_int[4]++;
244 if ( dy1 > 0.0 ) { x_arc[3][1] = del; y_arc[3][1] = dy1; arc_in[3] = 1; }
245 else { x_arc[4][1] = del; y_arc[4][1] = dy1; arc_in[4] = 1; }
246 }
247 if ( (-del ) > dx1 )
248 {
249 x_int[4][n_int[4]] = -del;
250 n_int[4]++;
251 if ( dy1 > 0.0 ) { x_arc[2][0] = -del; y_arc[2][0] = dy1; arc_in[2] = 1; }
252 else { x_arc[1][0] = -del; y_arc[1][0] = dy1; arc_in[1] = 1; }
253 }
254 }
255 }
256
257 n_int[0] = n_int[4];
258 n_int[5] = n_int[1];
259
260 y_int[0][0] = y_int[0][1] = dy1;
261 x_int[0][0] = x_int[4][0];
262 x_int[0][1] = x_int[4][1];
263 x_int[5][0] = x_int[5][1] = dx1;
264 y_int[5][0] = y_int[1][0];
265 y_int[5][1] = y_int[1][1];
266
267 /* There are five separate cases, depending of the number of vertices inside the circle */
268 switch ( n_in )
269 {
270 case 0:
271 {
272 /* We start with the whole area of the circle, and then subtract any bits that stick out. */
273 area = radsqu * mathematical::pi;
274 *perim_p = mathematical::twoPi * rad;
275 no_intersection = true;
276 for (n_vert = 1; n_vert < 5; n_vert++)
277 {
278 assert(n_int[n_vert] != 1);
279 if (n_int[n_vert] == 2)
280 {
281 /* The area of the bit to be subtracted is a sector minus a triangle. */
282 no_intersection = false;
283 angle = sector( x_int[n_vert][1], y_int[n_vert][1], x_int[n_vert][0], y_int[n_vert][0]);
284 area -= angle * radsqu * 0.5;
285 *perim_p -= angle * rad;
286 /* Two trinagles specified here, but one has zero area. */
287 area += ( - ( y_int[n_vert][1] - y_int[n_vert][0] ) * x_int[n_vert][0]
288 + ( x_int[n_vert][1] - x_int[n_vert][0] ) * y_int[n_vert][0] ) / 2.0;
289 }
290 }
291 /* Need to allow for when the circle is completely out side the rectanglle
292 by checking if the centre is outside the rectangle */
293 if ( no_intersection )
294 {
295 if ( (dx1>0) ||(dx2<0) || (dy1>0) || (dy2<0) )
296 {
297 *perim_p = *x_proj_edge_p = *y_proj_edge_p = 0.0;
298 area = *x_overlap_p = *y_overlap_p = 0.0;
299 return area;
300 }
301 }
302
303 break;
304 }
305
306 case 1:
307 {
308 /* Find which vertex is inside */
309 n_vert = 1;
310 while ( !vert_in[n_vert] ) { n_vert++; assert( n_vert < 5 ); }
311 assert( n_int[n_vert-1] == 1 );
312 if ( n_int[n_vert] != 1 )
313 {
314 assert( n_int[n_vert] == 1 );
315 }
316 angle = sector( x_int[n_vert-1][0], y_int[n_vert-1][0], x_int[n_vert][0], y_int[n_vert][0]);
317 area = angle * radsqu * 0.5;
318 *perim_p = angle * rad;
319 /* We subtract (or add) two triangles; the other two evaluate to zero */
320 area -= ( - ( x_int[n_vert][0] - dx[n_vert] ) * dy[n_vert]
321 + ( x_int[n_vert-1][0] - dx[n_vert] ) * dy[n_vert]
322 + ( y_int[n_vert][0] - dy[n_vert] ) * dx[n_vert]
323 - ( y_int[n_vert-1][0] - dy[n_vert] ) * dx[n_vert] ) / 2.0;
324
325 break;
326 }
327
328 case 2:
329 {
330 /* This time n_vert is the number of the side which is completely inside the circle */
331 n_vert = 1;
332 while ( !(vert_in[n_vert] && vert_in[n_vert+1]) ) { n_vert++; assert( n_vert < 5 ); }
333 assert( n_int[n_vert-1] == 1 );
334 assert( n_int[n_vert+1] == 1 );
335 angle = sector( x_int[n_vert-1][0], y_int[n_vert-1][0], x_int[n_vert+1][0], y_int[n_vert+1][0]);
336 area = angle * radsqu * 0.5;
337 *perim_p = angle * rad;
338 /* We subtract (or add) three triangles; the other three evaluate to zero */
339 area += ( ( x_int[n_vert+1][0] - dx[n_vert+1] ) * dy[n_vert+1]
340 - ( x_int[n_vert-1][0] - dx[n_vert] ) * dy[n_vert]
341 - ( y_int[n_vert+1][0] - dy[n_vert+1] ) * dx[n_vert+1]
342 + ( y_int[n_vert-1][0] - dy[n_vert] ) * dx[n_vert]
343 + ( dx[n_vert+1] -dx[n_vert] ) * dy[n_vert]
344 - ( dy[n_vert+1] -dy[n_vert] ) * dx[n_vert] ) / 2.0;
345
346 switch ( n_vert )
347 {
348 case 1: x_olap_min = dx1; break;
349 case 2: y_olap_max = dy2; break;
350 case 3: x_olap_max = dx2; break;
351 case 4: y_olap_min = dy1; break;
352 }
353
354 break;
355 }
356
357 case 3:
358 {
359 /* Find which vertex is NOT inside */
360 n_vert = 1;
361 while ( vert_in[n_vert] ) { n_vert++; assert( n_vert < 5 ); }
362 assert( n_int[n_vert-1] == 1 );
363 assert( n_int[n_vert] == 1 );
364 n_oppv = (n_vert + 2) % 4;
365 angle = sector( x_int[n_vert][0], y_int[n_vert][0], x_int[n_vert-1][0], y_int[n_vert-1][0]);
366 area = angle * radsqu * 0.5;
367 *perim_p = angle * rad;
368 /* We subtract (or add) four triangles; the other four evaluate to zero */
369 area += ( - ( x_int[n_vert][0] - dx[n_vert+1] ) * dy[n_vert+1]
370 + ( x_int[n_vert-1][0] - dx[n_vert-1] ) * dy[n_vert-1]
371 + ( y_int[n_vert][0] - dy[n_vert+1] ) * dx[n_vert+1]
372 - ( y_int[n_vert-1][0] - dy[n_vert-1] ) * dx[n_vert-1]
373 + ( dx[n_oppv] -dx[n_vert+1] ) * dy[n_oppv]
374 - ( dx[n_oppv] -dx[n_vert-1] ) * dy[n_oppv]
375 - ( dy[n_oppv] -dy[n_vert+1] ) * dx[n_oppv]
376 + ( dy[n_oppv] -dy[n_vert-1] ) * dx[n_oppv] ) / 2.0;
377
378 x_olap_min = dx1;
379 y_olap_max = dy2;
380 x_olap_max = dx2;
381 y_olap_min = dy1;
382
383 break;
384 }
385
386 case 4:
387 {
388 /* Easy! We have the whole rectangle. */
389 area = *x_overlap_p = *y_overlap_p = 1.0; // Normalised
390 *perim_p = *x_proj_edge_p = *y_proj_edge_p = 0.0;
391 return area;
392
393 break;
394 }
395 }
396
397 // The area may be very small negative by rounding errors
398 assert(area >=-1.0E-4);
399 if (area < 0.0) area = 0.0;
400 /* Return the overlap as a fraction of the rectangle's area. */
401 area /= ( (x2 - x1 ) * ( y2 - y1 ) );
402
403 // Sum the parts of the circumference that are inside the circle, projected onto axes
404 *x_proj_edge_p =
405 (
406 (y_arc[1][1] - y_arc[1][0]) * arc_in[1]
407 + (y_arc[2][1] - y_arc[2][0]) * arc_in[2]
408 + (y_arc[3][0] - y_arc[3][1]) * arc_in[3]
409 + (y_arc[4][0] - y_arc[4][1]) * arc_in[4]
410 );
411
412 *y_proj_edge_p =
413 (
414 (x_arc[1][0] - x_arc[1][1]) * arc_in[1]
415 + (x_arc[2][1] - x_arc[2][0]) * arc_in[2]
416 + (x_arc[3][1] - x_arc[3][0]) * arc_in[3]
417 + (x_arc[4][0] - x_arc[4][1]) * arc_in[4]
418 );
419
420 if (arc_in[1])
421 {
422 x_olap_min = min(x_olap_min, x_arc[1][1]);
423 x_olap_max = max(x_olap_max, x_arc[1][0]);
424 y_olap_min = min(y_olap_min, y_arc[1][0]);
425 y_olap_max = max(y_olap_max, y_arc[1][1]);
426 }
427 if (arc_in[2])
428 {
429 x_olap_min = min(x_olap_min, x_arc[2][0]);
430 x_olap_max = max(x_olap_max, x_arc[2][1]);
431 y_olap_min = min(y_olap_min, y_arc[2][0]);
432 y_olap_max = max(y_olap_max, y_arc[2][1]);
433 }
434 if (arc_in[3])
435 {
436 x_olap_min = min(x_olap_min, x_arc[3][0]);
437 x_olap_max = max(x_olap_max, x_arc[3][1]);
438 y_olap_min = min(y_olap_min, y_arc[3][1]);
439 y_olap_max = max(y_olap_max, y_arc[3][0]);
440 }
441 if (arc_in[4])
442 {
443 x_olap_min = min(x_olap_min, x_arc[4][1]);
444 x_olap_max = max(x_olap_max, x_arc[4][0]);
445 y_olap_min = min(y_olap_min, y_arc[4][1]);
446 y_olap_max = max(y_olap_max, y_arc[4][0]);
447 }
448
449 *x_overlap_p = ( x_olap_max - x_olap_min ) / ( x2 - x1 );
450 *y_overlap_p = ( y_olap_max - y_olap_min ) / ( y2 - y1 );
451 assert ( *x_overlap_p >= -floatSMALL );
452 assert ( *y_overlap_p >= -floatSMALL );
453
454 return area;
455} // End intersect
456
457
458// ************************************************************************* //
459
461(
462 double xc, double yc, double rad,
463 double x1, double x2,
464 double y1, double y2,
465 scalar* count_p, scalar* drag_p, scalar* centre_p
466)
467{
468 double xi = 0.0, lb, lb1, lb2, del;
469 bool within = true; // Indicates that the the intersection does not overlap the ends of the line
470
471 /* xi is the side we need to calc. intersections with */
472 if ( xc < x1 ) { xi = x1; }
473 else if ( xc > x2 ) { xi = x2; }
474
475 if ( xi == 0.0 )
476 {
477 del = rad; // The relevant lowest ( or highest) point is at end of vertical radius
478 }
479 else // The relevant lowest ( or highest) point at intersection with x = xi
480 {
481 del = rad*rad - ( xi - xc ) * ( xi - xc );
482 if ( del < 0.0 ) { del = 0.0; } // No intersection
483 else { del = std::sqrt(del); }
484 }
485
486 if ( ( yc + del ) > y2 ) { lb2 = y2; within = false; } else { lb2 = yc + del; }
487 if ( ( yc - del ) < y1 ) { lb1 = y1; within = false; } else { lb1 = yc - del; }
488
489 lb = (lb2 - lb1) / (y2 - y1);
490 *centre_p = (lb2 + lb1) * 0.5;
491
492 if ( lb < 0.0 ) lb = 0.0;
493
494 /* *count_p is 0 if the circle overlaps either y-side of the rectangle,
495 1 if the circle is entirely inside the rectangle
496 reduced if it overlaps x-sides.
497 A negative value indicates total blockage*/
498 if ( within && (lb > 0.0) )
499 {
500 *count_p = 1.0;
501 if ( ( xc - rad ) < x1 ) *count_p -= 0.5;
502 if ( ( xc + rad ) > x2 ) *count_p -= 0.5;
503 }
504 else
505 {
506 *count_p = 0.0;
507 }
508 *drag_p = lb * 1.2; //*drag_p = lb * CD_ROUND;
509 if ( lb > 0.99 ) { *count_p = -1000.0; *drag_p = 1000.0; }
510 assert(lb >-100.0);
511 return lb;
512}// End l_blockage
513
514
515// ************************************************************************* //
516
518(
519 double xc, double yc, double theta,
520 double wa, double wb,
521 double x1, double x2,
522 double y1, double y2,
523 scalar* count_p,
524 symmTensor2D& vdrag, scalar* perim_p,
525 scalar* x_lblk_p, scalar* y_lblk_p,
526 scalar* x_centre_p, scalar* y_centre_p
527)
528{
529 double x_int[6][2], y_int[6][2]; // Coordinates of intersections between the circle and sides of the rectangle.
530 double area, lpa, lpb, len;
531
532 double m = ::tan( theta );
533 double cth = ::cos( theta );
534 double sth = ::sin( theta );
535
536 double was = wa * sth * 0.5;
537 double wac = wa * cth * 0.5;
538 double wbs = wb * sth * 0.5;
539 double wbc = wb * cth * 0.5;
540 double waos = wa / sth * 0.5;
541 double waoc = wa / cth * 0.5;
542 double wbos = wb / sth * 0.5;
543 double wboc = wb / cth * 0.5;
544
545 double xb[6], yb[6], xp1, xp2, yp1, yp2;
546
547 double dx1 = (x1 - xc);
548 double dx2 = (x2 - xc);
549 double dy1 = (y1 - yc);
550 double dy2 = (y2 - yc);
551
552 *count_p = 0;
553
554// The vertices of the rectangle (all coordinates relative to centre of rectangle)
555 xb[1] = -wac - wbs;
556 yb[1] = -was + wbc;
557 xb[3] = wac + wbs;
558 yb[3] = was - wbc;
559 xb[2] = wac - wbs;
560 yb[2] = was + wbc;
561 xb[4] = -wac + wbs;
562 yb[4] = -was - wbc;
563
564 // First parameter of x_int or y_int determines which side of the cell we intersecting with
565 // Second parameter 0 is first intersection, 1 is second, going clockwise
566
567 if ( xb[1] < dx1 ) // i.e. if left corner of block is to the left of x1
568 {
569 // Where one of lower sides of block intersects with x=x1
570 // Innermost max determines which intersection is the genuine one
571 // (not if whole block is to left of x1)
572 y_int[1][0] = min(max(max(dx1 * m - wboc, -dx1 / m - waos), dy1), dy2);
573 // Upper intersection
574 y_int[1][1] = min(max(min(dx1 * m + wboc, -dx1 / m + waos), dy1), dy2);
575 }
576 else
577 {
578 y_int[1][1] = dy1;
579 y_int[1][0] = dy2;
580 // We add a quarter to count for each vertex inside the cell
581 if ( (yb[1] > dy1) && (yb[1] < dy2) ) // ?? Seems inefficient ??
582 { *count_p += 0.25; }
583 }
584 if ( xb[3] > dx2 )
585 {
586 y_int[3][1] = min(max(max(dx2 * m - wboc, -dx2 / m - waos), dy1), dy2);
587 y_int[3][0] = min(max(min(dx2 * m + wboc, -dx2 / m + waos), dy1), dy2);
588 }
589 else
590 {
591 y_int[3][0] = dy1;
592 y_int[3][1] = dy2;
593 if (yb[3] > dy1 && yb[3] < dy2)
594 {
595 *count_p += 0.25;
596 }
597 }
598 if (yb[2] > dy2)
599 {
600 x_int[2][0] = min(max(max(dy2 / m - wbos, -dy2 * m - waoc), dx1), dx2);
601 x_int[2][1] = min(max(min(dy2 / m + wbos, -dy2 * m + waoc), dx1), dx2);
602 }
603 else
604 {
605 x_int[2][0] = dx2;
606 x_int[2][1] = dx1;
607 if ( (xb[2] > dx1) && (xb[2] < dx2) )
608 { *count_p += 0.25; }
609 }
610 if ( yb[4] < dy1 )
611 {
612 x_int[4][1] = min(max(max(dy1 / m - wbos, -dy1 * m - waoc ), dx1), dx2);
613 x_int[4][0] = min(max(min(dy1 / m + wbos, -dy1 * m + waoc ), dx1), dx2);
614 }
615 else
616 {
617 x_int[4][1] = dx2;
618 x_int[4][0] = dx1;
619 if ( (xb[4] > dx1) && (xb[4] < dx2) )
620 { *count_p += 0.25; }
621 }
622
623 y_int[0][0] = y_int[0][1] = dy1;
624 x_int[0][0] = x_int[4][0];
625 x_int[0][1] = x_int[4][1];
626 x_int[5][0] = x_int[5][1] = dx1;
627 y_int[5][0] = y_int[1][0];
628 y_int[5][1] = y_int[1][1];
629
630
631// We can now define a smaller enclosing rectangle
632
633 xp1 = min(x_int[2][0], x_int[4][1]); // Leftmost of the intersections with top and bottom of cell
634 if ( yb[1] > dy1 && yb[1] < dy2 ) xp1 = min(xp1, xb[1] ); // left corner of block
635 xp1 = max(xp1, dx1); // Make sure it is not to the left of x1
636
637 yp2 = max(y_int[1][1], y_int[3][0] );
638 if ( xb[2] > dx1 && xb[2] < dx2 ) yp2 = max(yp2, yb[2] );
639 yp2 = min(yp2, dy2);
640
641 xp2 = max(x_int[2][1], x_int[4][0] );
642 if ( yb[3] > dy1 && yb[3] < dy2 ) xp2 = max(xp2, xb[3] );
643 xp2 = min(xp2, dx2);
644
645 yp1 = min(y_int[1][0], y_int[3][1]);
646 if ( xb[4] > dx1 && xb[4] < dx2 ) yp1 = min(yp1, yb[4] );
647 yp1 = max(yp1, dy1 );
648
649 // Conveniently, the dimensions of the enclosing rectangle give us the line blockages
650 *x_lblk_p = (xp2 - xp1 ) / (x2 - x1 );
651 if ( *x_lblk_p < 0.0 ) { *x_lblk_p = 0.0; *count_p = 0.0; }; // ?? Better to trap no intersection earlier??
652 *y_lblk_p = (yp2 - yp1 ) / (y2 - y1 );
653 if ( *y_lblk_p < 0.0 ) { *y_lblk_p = 0.0; *count_p = 0.0; };
654
655 *x_centre_p = xc + (xp2 + xp1 ) * 0.5;
656 *y_centre_p = yc + (yp2 + yp1 ) * 0.5;
657
658 *perim_p = lpa = lpb = 0.0;;
659 area = (xp2 - xp1 ) * ( yp2 - yp1 );
660 {
661 double dxx, dyy;
662 // Lower left
663 dyy = max(0.0, min(yb[1], y_int[1][0]) - yp1);
664 dxx = min(xb[4], x_int[0][1] ) - xp1;
665
666 if ( ( dxx * dyy) > 0.0 )
667 {
668 area -= dxx * dyy * 0.5;
669 len = std::hypot(dxx, dyy);
670 lpa += len * 0.5;
671 *perim_p += len;
672 }
673 // Upper left
674 dxx = max(0.0, min(xb[2], x_int[2][0]) - xp1);
675 dyy = yp2 - max(yb[1], y_int[1][1] );
676 if ( ( dxx * dyy) > 0.0 )
677 {
678 area -= dxx * dyy * 0.5;
679 len = std::hypot(dxx, dyy);
680 lpb += len * 0.5;
681 *perim_p += len;
682 }
683 // Upper right
684 dyy = max(0.0, yp2 - max(yb[3], y_int[3][0]));
685 dxx = xp2 - max(xb[2], x_int[2][1] );
686 if ( ( dxx * dyy) > 0.0 )
687 {
688 area -= dxx * dyy * 0.5;
689 len = std::hypot(dxx, dyy);
690 lpa += len * 0.5;
691 *perim_p += len;
692 }
693 // Lower right
694 dxx = max(0.0, xp2 - max(xb[4], x_int[4][0]));
695 dyy = min(yb[3], y_int[3][1] ) - yp1;
696 if ( ( dxx * dyy) > 0.0 )
697 {
698 area -= dxx * dyy * 0.5;
699 len = std::hypot(dxx, dyy);
700 lpb += len * 0.5;
701 *perim_p += len;
702 }
703
704 }
705
706 vdrag.xx() = lpa * cth * cth + lpb * sth * sth;
707 vdrag.xy() = lpa * cth * sth - lpb * sth * cth;
708 vdrag.yy() = lpa * sth * sth + lpb * cth * cth;
709
710 return area / ( (x2 - x1 ) * ( y2 - y1 ) );
711} // End inters_db
712
713
714// ************************************************************************* //
Preparation of fields for PDRFoam.
#define floatSMALL
Definition: PDRsetFields.H:54
Y[inertIndex] max(0.0)
radiation::radiationModel & rad
A templated (2 x 2) symmetric tensor of objects of <T>, effectively containing 3 elements,...
Definition: SymmTensor2D.H:61
const Cmpt & xx() const
Definition: SymmTensor2DI.H:80
const Cmpt & xy() const
Definition: SymmTensor2DI.H:86
const Cmpt & yy() const
Definition: SymmTensor2DI.H:98
double l_blockage(double xc, double yc, double rad, double x1, double x2, double y1, double y2, scalar *count_p, scalar *drag_p, scalar *centre_p)
double inters_db(double xc, double yc, double theta, double wa, double wb, double x1, double x2, double y1, double y2, scalar *count_p, symmTensor2D &vdrag, scalar *perim_p, scalar *x_lblk, scalar *y_lblk, scalar *x_centre_p, scalar *y_centre_p)
The area overlap in the plane of a diagonal block and a cell.
double inters_cy(double xc, double yc, double rad, double x1, double x2, double y1, double y2, scalar *perim_p, scalar *x_proj_edge_p, scalar *y_proj_edge_p, scalar *x_overlap_p, scalar *y_overlap_p)
Area of intersection between circle and rectangle.
Different types of constants.
const wordList area
Standard area field types (scalar, vector, tensor, etc)
Namespace for OpenFOAM.
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
dimensionedScalar y1(const dimensionedScalar &ds)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33