PDRutilsOverlap.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-2021 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 // A sign-corrected multiply
47 // This is used for porosity of obstacle intersections
48 inline static scalar COMBLK(const scalar a, const scalar b)
49 {
50 if (a < 0)
51 {
52 return -a * b;
53 }
54
55 return a * b;
56 }
57
58
59 // Obstacle satisfies some minimum size checks.
60 // A volume check misses thin plates, so use area.
61 // Thin sheet overlaps can be produced by touching objects
62 // if the obs_extend parameter is > 0.
63 inline static bool obsHasMinSize(const vector& span, const PDRparams& tol)
64 {
65 return
66 (
67 (cmptProduct(span) > tol.min_overlap_vol)
68 &&
69 (
70 (span.x() * span.y() > tol.min_overlap_area)
71 || (span.y() * span.z() > tol.min_overlap_area)
72 || (span.z() * span.x() > tol.min_overlap_area)
73 )
74 );
75 }
76
77} // End namespace Foam
78
79
80// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
81
82
84(
85 scalar xmin,
86 scalar xmax,
87 const PDRblock::location& grid,
88 List<scalar>& olap,
89 int *cmin, int *cmax,
90 int *cfmin, int *cfmax
91)
92{
93 // Looking at one coordinate direction, called x here, for something
94 // that extends from xmin to xmax, calculate how much it overlaps
95 // each cell in this direction. Result returned in 'olap' List is
96 // the proportion of the grid step overlapped, i.e dimensionless.
97 // First and last steps overlapped given by *cmin, *cmax
98 // Ditto for shifted grid given by *cfmin, *cfmax.
99
100 // Initially zero everywhere
101 olap = Zero;
102
103 if (olap.size() < grid.nPoints())
104 {
106 << "The overlap scratch array is too small, has "
107 << olap.size() << " but needs " << grid.nPoints() << nl
108 << exit(FatalError);
109 }
110
111
112 // No intersection with the box
113 if (xmax <= grid.first() || grid.last() <= xmin)
114 {
115 // Mark as bad range, cannot iterate
116 *cmin = 0;
117 *cmax = -1;
118
119 // Another bad range (cannot iterate) but for extra safety ensure
120 // that (cfmin -> cmin) and (cmax -> cfmax) cannot iterate either
121 *cfmin = 1;
122 *cfmax = -2;
123 return;
124 }
125
126 // Ensure search is within the (point) bounds
127 xmin = grid.clip(xmin);
128 xmax = grid.clip(xmax);
129
130 // The begin/end of the obstacle
131 *cmin = grid.findCell(xmin);
132 *cmax = grid.findCell(xmax);
133
134 for (label ix = *cmin; ix <= *cmax; ++ix)
135 {
136 olap[ix] = 1.0;
137 }
138
139 // Fixup ends
140 if (*cmax == *cmin)
141 {
142 olap[*cmax] = (xmax - xmin) / grid.width(*cmax);
143 }
144 else
145 {
146 if (grid[*cmin] < xmin)
147 {
148 olap[*cmin] = (grid[*cmin+1] - xmin) / grid.width(*cmin);
149 }
150
151 if (xmax < grid[*cmax+1])
152 {
153 olap[*cmax] = (xmax - grid[*cmax]) / grid.width(*cmax);
154 }
155 }
156 assert(olap[*cmax] >= 0.0);
157
158
159 // Is xmin below/above the cell-centre (for virtual staggered-grid) ?
160 *cfmin =
161 (
162 xmin < grid.C(*cmin)
163 ? *cmin
164 : Foam::min(*cmin+1, grid.nCells()-1)
165 );
166
167 // Is xmax below/above the cell-centre (for virtual staggered-grid) ?
168 *cfmax =
169 (
170 xmax < grid.C(*cmax)
171 ? *cmax
172 : Foam::min(*cmax+1, grid.nCells()-1)
173 );
174}
175
176
177/**************************************************************************************************/
178
180(
181 const UList<scalar>& a_olap, label amin, label amax,
182 const UList<scalar>& b_olap, label bmin, label bmax,
183 SquareMatrix<scalar>& ab_olap
184)
185{
186 // We go one over the relevant min/max limits since these values might be
187 // used. If not, they would have been zeroed in one_d_overlap
188
189 amin = Foam::max(0, amin-1);
190 bmin = Foam::max(0, bmin-1);
191 amax = Foam::min(a_olap.size()-1, amax+1);
192 bmax = Foam::min(b_olap.size()-1, bmax+1);
193
194 for (label ia = amin; ia <= amax; ++ia)
195 {
196 for (label ib = bmin; ib <= bmax; ++ib)
197 {
198 ab_olap(ia,ib) = a_olap[ia] * b_olap[ib];
199 }
200 }
201}
202
203
204/**************************************************************************************************/
205
207(
208 scalar ac, scalar bc, scalar dia,
209 scalar theta, scalar wa, scalar wb,
210 const PDRblock::location& agrid, label amin, label amax,
211 const PDRblock::location& bgrid, label bmin, label bmax,
212 SquareMatrix<scalar>& ab_olap,
213 SquareMatrix<scalar>& ab_perim,
214 SquareMatrix<scalar>& a_lblock,
215 SquareMatrix<scalar>& ac_lblock,
216 SquareMatrix<scalar>& c_count,
218 SquareMatrix<scalar>& b_lblock,
219 SquareMatrix<scalar>& bc_lblock
220)
221{
222 /* This routine calculates the proportion of each (two-dimensional) grid cell
223 overlapped by the circle or angled rectangle. Coordinates are labelled a and b.
224 On entry:
225 ac, bc coordinates of centre of circle or rectangle
226 dia diameter of circle (zero for rectangle)
227 theta, wa, wb parameters for rectangle
228 agrid[] locations of grid lines of a-grid
229 amin, amax first and last cells in a-grid overlapped by object
230 (similarly for b)
231 On exit:
232 abolap 2-D array of (proportionate) area blockage by grid cell
233 a_lblock 2-D array of (proportionate) blockage to a-direction flow
234 (This will be area blockage when extruded in the third coordinate).
235 a_count (2-D array)The contribution of this object to the count of obstacles blocking
236 a-direction flow. This is only non-zero if the object is inside the
237 lateral boundaries of the cell. It is large negative if the cell is
238 totally blocked in this direction.
239 (similarly for b)
240 c_drag 2-D array of tensor that will give tensor drag in each cell (when multiplied
241 Cd, cylinder length, and 0.5 rho*U^2) Dimension: L.
242
243 Note that this routine does not zero array elements outside the amin to amax, bmin to bmax area.
244 */
245
246 scalar count, a_lblk, b_lblk, perim, dummy;
247
248 symmTensor2D vdrag(Zero);
249
250 // Prevent stepping outside of the array when the obstacle is on the
251 // upper boundary
252
253 // Upper limit of inclusive range is nCells-1
254 amin = Foam::max(0, amin);
255 bmin = Foam::max(0, bmin);
256 amax = Foam::min(amax, agrid.nCells()-1);
257 bmax = Foam::min(bmax, bgrid.nCells()-1);
258
259 for (label ia = amin; ia <= amax; ++ia)
260 {
261 // Cell-centred grid
262 const scalar a1 = agrid[ia];
263 const scalar a2 = agrid[ia+1];
264
265 // Left-shifted staggered face grid (-1 addressing is OK)
266 const scalar af1 = agrid.C(ia-1);
267 const scalar af2 = agrid.C(ia);
268
269 for (label ib = bmin; ib <= bmax; ++ib)
270 {
271 // Cell-centred grid
272 const scalar b1 = bgrid[ib];
273 const scalar b2 = bgrid[ib+1];
274
275 // Left-shifted staggered face grid (-1 addressing is OK)
276 const scalar bf1 = bgrid.C(ib-1);
277 const scalar bf2 = bgrid.C(ib);
278
279 // Do the centred cell
280 if ( dia > 0.0 )
281 {
282 ab_olap(ia,ib) = inters_cy
283 (
284 ac, bc, 0.5*dia, a1, a2, b1, b2, &perim,
285 &dummy, &dummy, &b_lblk, &a_lblk
286 );
287/* The last two arguments of the above call appear to be reversed, but the inters_cy routine returns
288 the amount of overlap in the a and b direcvtions, which are the blockage to the b and a directions. */
289
290/* abolap * cell area is area of cylinder in this cell. Divide by PI%D^2/4 to get proportion of cylinder in cell
291 For whole cylinger c_drag should be = D, so multiply by D. */
292
293 c_drag(ia,ib).xx() = c_drag(ia,ib).yy() = 4.0 * ab_olap(ia,ib) * (a2 - a1) * (b2 - b1) / dia / mathematical::pi;
294 c_drag(ia,ib).xy() = Zero;
295 c_count(ia,ib) = perim / (mathematical::pi * dia);
296
297//******?????
298 scalar area = (a2 - a1) * (b2 - b1);
299 scalar rat = dia * dia / area - 1.5;
300 if (rat > 0.0)
301 {
302 scalar da = ac - 0.5 * (a1 + a2);
303 scalar db = bc - 0.5 * (b1 + b2);
304 scalar dc = std::hypot(da, db);
305 scalar rat1 = min(max((dc / sqrt(area) - 0.3) * 1.4, 0), 1);
306 scalar drg0 = c_drag(ia,ib).xx();
307 scalar drg1 = c_drag(ia,ib).yy();
308 scalar drg = std::hypot(drg0, drg1);
309 c_drag(ia,ib).xx() = drg * ( 1.0 - rat1 ) + drg * da*da/dc/dc * rat1;
310 c_drag(ia,ib).yy() = drg * ( 1.0 - rat1 ) + drg * db*db/dc/dc * rat1;
311 c_drag(ia,ib).xy() = drg * da*db/dc/dc *rat1;
312 }
313 }
314 else
315 {
316 ab_olap(ia,ib) = inters_db( ac, bc, theta, wa, wb, a1, a2, b1, b2, &count, c_drag(ia,ib), &perim, &a_lblk, &b_lblk, &dummy, &dummy );
317 c_count(ia,ib) = perim / ( wa + wb ) * 0.5;
318 }
319 ac_lblock(ia,ib) = a_lblk;
320 bc_lblock(ia,ib) = b_lblk;
321 ab_perim(ia,ib) = perim;
322
323 // Do the a-shifted cell
324 if ( dia > 0.0 ) // I.e. a cylinder, not a d.b.
325 {
326 if (ac >= af1 && ac < af2)
327 {
328 // Only want to block one layer of faces
329 a_lblock(ia,ib) = l_blockage
330 (
331 ac, bc, 0.5*dia,
332 af1, af2, b1, b2, &count, &dummy, &dummy
333 );
334 }
336 (
337 ac, bc, 0.5*dia,
338 af1, af2, b1, b2,
339 &perim, &count, &dummy, &dummy, &dummy
340 );
341 }
342 else
343 {
345 (
346 ac, bc, theta, wa, wb, af1, af2, b1, b2,
347 &count, vdrag, &dummy, &a_lblk, &b_lblk, &dummy, &dummy
348 );
349 a_lblock(ia,ib) = a_lblk;
350 }
351
352 // Do the b-shifted cell
353 if ( dia > 0.0 )
354 {
355 if (bc >= bf1 && bc < bf2)
356 {
357 // Only want to block one layer of faces
358 b_lblock(ia,ib) = l_blockage
359 (
360 bc, ac, 0.5*dia, bf1, bf2, a1, a2,
361 &count, &(vdrag.yy()), &dummy
362 );
363 }
364
366 (
367 ac, bc, 0.5*dia,
368 a1, a2, bf1, bf2,
369 &perim, &dummy, &count, &dummy, &dummy
370 );
371 }
372 else
373 {
375 (
376 ac, bc, theta, wa, wb,
377 a1, a2, bf1, bf2,
378 &count, vdrag, &dummy, &a_lblk, &b_lblk, &dummy, &dummy
379 );
380 b_lblock(ia,ib) = b_lblk;
381 }
382 }
383 }
384
385} // End circle_overlap
386
387
388/**************************************************************************************************/
389
390scalar block_overlap
391(
393 const labelRange& range,
394 const scalar multiplier
395)
396{
397 // Size information
398 const label nBlock = range.size();
399
400 // The return value
401 scalar totVolume = 0;
402
403 if (nBlock < 2) return 0;
404
405
406 // Sort blocks by their x-position (with sortBias)
407 labelList blkOrder;
408 sortedOrder(blocks.slice(range), blkOrder);
409
410 DynamicList<PDRobstacle> newBlocks;
411
412 // Work through the sorted blocks
413 for (label i1 = 0; i1 < nBlock-1; ++i1)
414 {
415 const PDRobstacle& blk1 = blocks[range[blkOrder[i1]]];
416
417 // Upper coordinates
418 const vector max1 = blk1.pt + blk1.span;
419
420 // For second block start with the next one on the list, and
421 // stop when we find the first one whose biased x-position
422 // is beyond the end of the block1
423
424 for (label i2 = i1 + 1; i2 < nBlock; ++i2)
425 {
426 const PDRobstacle& blk2 = blocks[range[blkOrder[i2]]];
427
428 // Upper coordinates
429 const vector max2 = blk2.pt + blk2.span;
430
431 if (max1.x() <= blk2.x())
432 {
433 break;
434 }
435
436 if
437 (
438 max1.y() <= blk2.y()
439 || max1.z() <= blk2.z()
440 || max2.y() <= blk1.y()
441 || max2.z() <= blk1.z()
442 || (blk1.vbkge * blk2.vbkge <= 0)
443 )
444 {
445 continue;
446 }
447
448
449 {
450 PDRobstacle over;
451
452 over.pt = max(blk1.pt, blk2.pt);
453 over.span = min(max1, max2) - over.pt;
454
455 assert(cmptProduct(over.span) > 0.0);
456
457 // This routine should only have been called for all +ve o r all -ve obstacles
458 assert(blk1.vbkge * blk2.vbkge > 0);
459 /* At the first level of intersection, we create an obstacle of blockage -1 (if both objects solid)
460 to cancel out the double counting. (multiplier is 1).
461 ?? COMBLK does a (sign corrected) multiply; is this corrrect for porous obstacles?
462 Depends on how blockages were summed in the first place. In fact this -ve obstacle
463 concept only works if the blockages are summed??*/
464 over.vbkge = - COMBLK( blk1.vbkge, blk2.vbkge ) * multiplier;
465 over.xbkge = - COMBLK( blk1.xbkge, blk2.xbkge ) * multiplier;
466 over.ybkge = - COMBLK( blk1.ybkge, blk2.ybkge ) * multiplier;
467 over.zbkge = - COMBLK( blk1.zbkge, blk2.zbkge ) * multiplier;
468 over.typeId = 81 + int(15 * multiplier); // Not subsequently used
469
470 if (obsHasMinSize(over.span, pars))
471 {
472 // Obstacle satisfies some minimum size checks
473 totVolume -= over.volume();
474
475 newBlocks.append(over);
476 }
477 }
478 }
479 }
480
481 blocks.append(std::move(newBlocks));
482
483 return totVolume;
484}
485
486
487/**************************************************************************************************/
488
489using namespace Foam::PDRutils;
490
492(
494 const labelRange& range,
495 const UList<PDRobstacle>& cylinders
496)
497{
498 // Size information
499 const label nBlock = range.size();
500 const label nCyl = cylinders.size();
501
502 // The return value
503 scalar totVolume = 0;
504
505 if (!nBlock || !nCyl) return 0;
506
507 scalar area, a_lblk, b_lblk, dummy, a_centre, b_centre;
508 symmTensor2D dum2;
509
510
511 // Sort blocks and cylinders by their x-position (with sortBias)
512 labelList blkOrder;
513 sortedOrder(blocks.slice(range), blkOrder);
514
515 labelList cylOrder;
516 sortedOrder(cylinders, cylOrder);
517
518 DynamicList<PDRobstacle> newBlocks;
519
520 // Work through the sorted blocks
521 for (label i1 = 0; i1 < nBlock; i1++)
522 {
523 const PDRobstacle& blk1 = blocks[range[blkOrder[i1]]];
524
525 // Upper coordinates
526 const vector max1 = blk1.pt + blk1.span;
527
528 // Cyls whose end is before start of this block no longer
529 // need to be considered
530
531 label i2 = 0;
532 while (i2 < nCyl-1 && cylinders[cylOrder[i2]] < blk1)
533 {
534 ++i2;
535 }
536
537 for (/*nil*/; i2 < nCyl; ++i2)
538 {
539 const PDRobstacle& cyl2 = cylinders[cylOrder[i2]];
540
541 // Calculate overlap in axis direction; if zero continue.
542 // Calculate 2-d overlap and c 0f g; if area zero continue.
543
544 PDRobstacle over;
545
546
547 switch (cyl2.orient)
548 {
549 case vector::Z:
550 {
551 const scalar zm2 = cyl2.z() + cyl2.len();
552 if (blk1.z() > zm2 || cyl2.z() > max1.z()) continue;
553
554 if ( cyl2.dia() == 0.0 )
555 {
557 (
558 cyl2.x(), cyl2.y(), cyl2.theta(), cyl2.wa, cyl2.wb,
559 blk1.x(), max1.x(),
560 blk1.y(), max1.y(),
561 &dummy, dum2, &dummy, &a_lblk, &b_lblk,
562 &a_centre, &b_centre
563 );
564 }
565 else
566 {
568 (
569 cyl2.x(), cyl2.y(), 0.5*cyl2.dia(),
570 blk1.x(), max1.x(),
571 blk1.y(), max1.y(),
572 &dummy, &dummy, &dummy, &dummy, &dummy
573 );
574 b_lblk = l_blockage
575 (
576 cyl2.x(), cyl2.y(), 0.5*cyl2.dia(),
577 blk1.x(), max1.x(),
578 blk1.y(), max1.y(),
579 &dummy, &dummy, &b_centre
580 );
581 a_lblk = l_blockage
582 (
583 cyl2.y(), cyl2.x(), 0.5*cyl2.dia(),
584 blk1.y(), max1.y(),
585 blk1.x(), max1.x(),
586 &dummy, &dummy, &a_centre
587 );
588 }
589 if (equal(area, 0)) continue;
590 assert(a_lblk >0.0);
591 assert(b_lblk >0.0);
592
593 // The intersection between a circle and a rectangle can be an odd shape.
594 // We have its area. a_lblk and b_lblk are dimensions of enclosing rectangle
595 // and a_centre and b_centre its centre. We scale this rectangle down to
596 // the correct areacorrect area, as a rectangular approximation to the intersection.
597 const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
598
599 a_lblk *= blk1.span.x() * ratio;
600 b_lblk *= blk1.span.y() * ratio;
601 assert(b_lblk >0.0);
602 assert(a_lblk >0.0);
603
604 over.x() = a_centre - 0.5 * a_lblk;
605 over.y() = b_centre - 0.5 * b_lblk;
606 over.z() = max(blk1.z(), cyl2.z());
607
608 over.span.x() = a_lblk;
609 over.span.y() = b_lblk;
610 over.span.z() = min(max1.z(), cyl2.z() + cyl2.len()) - over.z();
611 assert(over.x() > -200.0);
612 assert(over.x() < 2000.0);
613 }
614 break;
615
616 case vector::Y:
617 {
618 const scalar ym2 = cyl2.y() + cyl2.len();
619 if (blk1.y() > ym2 || cyl2.y() > max1.y()) continue;
620
621 if ( cyl2.dia() == 0.0 )
622 {
624 (
625 cyl2.z(), cyl2.x(), cyl2.theta(), cyl2.wa, cyl2.wb,
626 blk1.z(), max1.z(),
627 blk1.x(), max1.x(),
628 &dummy, dum2, &dummy, &a_lblk, &b_lblk,
629 &a_centre, &b_centre
630 );
631 }
632 else
633 {
635 (
636 cyl2.z(), cyl2.x(), 0.5*cyl2.dia(),
637 blk1.z(), max1.z(),
638 blk1.x(), max1.x(),
639 &dummy, &dummy, &dummy, &dummy, &dummy
640 );
641
642 b_lblk = l_blockage
643 (
644 cyl2.z(), cyl2.x(), 0.5*cyl2.dia(),
645 blk1.z(), max1.z(),
646 blk1.x(), max1.x(),
647 &dummy, &dummy, &b_centre
648 );
649
650 a_lblk = l_blockage
651 (
652 cyl2.x(), cyl2.z(), 0.5*cyl2.dia(),
653 blk1.x(), max1.x(),
654 blk1.z(), max1.z(),
655 &dummy, &dummy, &a_centre
656 );
657 }
658
659 if (equal(area, 0)) continue;
660 assert(a_lblk >0.0);
661 assert(b_lblk >0.0);
662
663 // a_lblk and b_lblk are dimensions of enclosing rectangle.
664 // Need to scale to correct area
665 const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
666 a_lblk *= blk1.span.z() * ratio;
667 b_lblk *= blk1.span.x() * ratio;
668
669 over.z() = a_centre - a_lblk * 0.5;
670 over.x() = b_centre - b_lblk * 0.5;
671 over.y() = max(blk1.y(), cyl2.y());
672
673 over.span.z() = a_lblk;
674 over.span.x() = b_lblk;
675 over.span.y() = min(max1.y(), cyl2.y() + cyl2.len()) - over.y();
676 }
677 break;
678
679 case vector::X:
680 {
681 const scalar xm2 = cyl2.x() + cyl2.len();
682 if (blk1.x() > xm2 || cyl2.x() > max1.x()) continue;
683
684 if ( cyl2.dia() == 0.0 )
685 {
687 (
688 cyl2.y(), cyl2.z(), cyl2.theta(), cyl2.wa, cyl2.wb,
689 blk1.y(), max1.y(),
690 blk1.z(), max1.z(),
691 &dummy, dum2, &dummy, &a_lblk, &b_lblk,
692 &a_centre, &b_centre
693 );
694 }
695 else
696 {
698 (
699 cyl2.y(), cyl2.z(), 0.5*cyl2.dia(),
700 blk1.y(), max1.y(),
701 blk1.z(), max1.z(),
702 &dummy, &dummy, &dummy, &dummy, &dummy
703 );
704
705 b_lblk = l_blockage
706 (
707 cyl2.y(), cyl2.z(), 0.5*cyl2.dia(),
708 blk1.y(), max1.y(),
709 blk1.z(), max1.z(),
710 &dummy, &dummy, &b_centre
711 );
712
713 a_lblk = l_blockage
714 (
715 cyl2.z(), cyl2.y(), 0.5*cyl2.dia(),
716 blk1.z(), max1.z(),
717 blk1.y(), max1.y(),
718 &dummy, &dummy, &a_centre
719 );
720
721 }
722
723 if (equal(area, 0)) continue;
724 assert(a_lblk >0.0);
725 assert(b_lblk >0.0);
726
727 // a_lblk and b_lblk are dimensions of enclosing rectangle.
728 // Need to scale to correct area
729 const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
730 assert(ratio >-10000.0);
731 assert(ratio <10000.0);
732 a_lblk *= blk1.span.y() * ratio;
733 b_lblk *= blk1.span.z() * ratio;
734
735 over.y() = a_centre - a_lblk * 0.5;
736 over.z() = b_centre - b_lblk * 0.5;
737 over.x() = max(blk1.x(), cyl2.x());
738
739 over.span.y() = a_lblk;
740 over.span.z() = b_lblk;
741 over.span.x() = min(max1.x(), cyl2.x() + cyl2.len()) - over.x();
742 }
743 break;
744 }
745 over.vbkge = over.xbkge = over.ybkge = over.zbkge = -1.0;
746 over.typeId = PDRobstacle::IGNORE;
747
748 assert(cmptProduct(over.span) > 0.0);
749 assert(b_lblk >0.0);
750 assert(a_lblk >0.0);
751 assert(over.x() > -10000.0);
752
753 if (obsHasMinSize(over.span, pars))
754 {
755 // Obstacle satisfies some minimum size checks
756 totVolume -= over.volume();
757
758 newBlocks.append(over);
759 }
760 }
761 }
762
763 blocks.append(std::move(newBlocks));
764
765 return totVolume;
766}
767
768
769// ************************************************************************* //
scalar range
Preparation of fields for PDRFoam.
scalar block_overlap(DynamicList< PDRobstacle > &blocks, const labelRange &range, const scalar multiplier=1.0)
Calculate block/block overlaps.
scalar block_cylinder_overlap(DynamicList< PDRobstacle > &blocks, const labelRange &range, const UList< PDRobstacle > &cylinders)
Calculate block/cylinder overlaps.
Y[inertIndex] max(0.0)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
Grid locations in an axis direction.
Definition: PDRblock.H:181
scalar C(const label i) const
Cell centre at element position.
Definition: PDRblockI.H:101
label findCell(const scalar p) const
Find the cell index enclosing this location.
label nCells() const
The number of cells in this direction.
Definition: PDRblockI.H:37
label nPoints() const
The number of points in this direction.
Definition: PDRblockI.H:43
scalar width(const label i) const
Cell size at element position.
Definition: PDRblockI.H:91
const scalar & clip(const scalar &val) const
Definition: PDRblockI.H:127
Obstacle definitions for PDR.
Definition: PDRobstacle.H:75
scalar volume() const
Volume of the obstacle.
point pt
The obstacle location.
Definition: PDRobstacle.H:123
direction orient
The x/y/z orientation (0,1,2)
Definition: PDRobstacle.H:116
scalar x() const
Obstacle position accessors.
Definition: PDRobstacle.H:225
scalar y() const
Definition: PDRobstacle.H:226
scalar z() const
Definition: PDRobstacle.H:227
scalar theta() const
Definition: PDRobstacle.H:131
scalar len() const
Definition: PDRobstacle.H:132
scalar dia() const
Definition: PDRobstacle.H:130
vector span
The obstacle dimensions (for boxes)
Definition: PDRobstacle.H:126
int typeId
The obstacle type-id.
Definition: PDRobstacle.H:113
Parameters for PDRsetFields.
Definition: PDRparams.H:56
scalar min_overlap_vol
Definition: PDRparams.H:109
scalar min_overlap_area
Definition: PDRparams.H:110
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
Definition: SquareMatrix.H:66
A templated (2 x 2) symmetric tensor of objects of <T>, effectively containing 3 elements,...
Definition: SymmTensor2D.H:61
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
T & first()
Return the first element of the list.
Definition: UListI.H:202
SubList< T > slice(const label pos, label len=-1)
Return SubList slice (non-const access) - no range checking.
Definition: SubList.H:165
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
A range or interval of labels defined by a start and a size.
Definition: labelRange.H:58
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
Definition: BitOps.H:78
bool equal(const Matrix< Form1, Type > &A, const Matrix< Form2, Type > &B, const bool verbose=false, const label maxDiffs=10, const scalar relTol=1e-5, const scalar absTol=1e-8)
Compare matrix elements for absolute or relative equality.
Definition: MatrixTools.C:34
Utilities for PDR (eg, for setFields)
void two_d_overlap(const UList< scalar > &a_olap, label amin, label amax, const UList< scalar > &b_olap, label bmin, label bmax, SquareMatrix< scalar > &ab_olap)
Combine two 1D overlaps.
void one_d_overlap(scalar xmin, scalar xmax, const PDRblock::location &grid, List< scalar > &olap, int *cmin, int *cmax, int *cfmin, int *cfmax)
Determine 1-D overlap locations for a geometric entity.
void circle_overlap(scalar ac, scalar bc, scalar dia, scalar theta, scalar wa, scalar wb, const PDRblock::location &agrid, label amin, label amax, const PDRblock::location &bgrid, label bmin, label bmax, SquareMatrix< scalar > &ab_olap, SquareMatrix< scalar > &ab_perim, SquareMatrix< scalar > &a_lblock, SquareMatrix< scalar > &ac_lblock, SquareMatrix< scalar > &c_count, SquareMatrix< symmTensor2D > &c_drag, SquareMatrix< scalar > &b_lblock, SquareMatrix< scalar > &bc_lblock)
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.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Cmpt cmptProduct(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:598
dimensionedScalar sqrt(const dimensionedScalar &ds)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & b
Definition: createFields.H:27