48 inline static scalar COMBLK(
const scalar a,
const scalar
b)
63 inline static bool obsHasMinSize(
const vector& span,
const PDRparams& tol)
90 int *cfmin,
int *cfmax
106 <<
"The overlap scratch array is too small, has "
113 if (xmax <= grid.
first() || grid.
last() <= xmin)
127 xmin = grid.
clip(xmin);
128 xmax = grid.
clip(xmax);
134 for (label ix = *cmin; ix <= *cmax; ++ix)
142 olap[*cmax] = (xmax - xmin) / grid.
width(*cmax);
146 if (grid[*cmin] < xmin)
148 olap[*cmin] = (grid[*cmin+1] - xmin) / grid.
width(*cmin);
151 if (xmax < grid[*cmax+1])
153 olap[*cmax] = (xmax - grid[*cmax]) / grid.
width(*cmax);
156 assert(olap[*cmax] >= 0.0);
194 for (label ia = amin; ia <= amax; ++ia)
196 for (label ib = bmin; ib <= bmax; ++ib)
198 ab_olap(ia,ib) = a_olap[ia] * b_olap[ib];
208 scalar ac, scalar bc, scalar dia,
209 scalar theta, scalar wa, scalar wb,
246 scalar
count, a_lblk, b_lblk, perim, dummy;
259 for (label ia = amin; ia <= amax; ++ia)
262 const scalar a1 = agrid[ia];
263 const scalar a2 = agrid[ia+1];
266 const scalar af1 = agrid.
C(ia-1);
267 const scalar af2 = agrid.
C(ia);
269 for (label ib = bmin; ib <= bmax; ++ib)
272 const scalar b1 = bgrid[ib];
273 const scalar b2 = bgrid[ib+1];
276 const scalar bf1 = bgrid.
C(ib-1);
277 const scalar bf2 = bgrid.
C(ib);
284 ac, bc, 0.5*dia, a1, a2, b1, b2, &perim,
285 &dummy, &dummy, &b_lblk, &a_lblk
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);
298 scalar
area = (a2 - a1) * (b2 - b1);
299 scalar rat = dia * dia /
area - 1.5;
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;
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;
319 ac_lblock(ia,ib) = a_lblk;
320 bc_lblock(ia,ib) = b_lblk;
321 ab_perim(ia,ib) = perim;
326 if (ac >= af1 && ac < af2)
332 af1, af2, b1, b2, &count, &dummy, &dummy
339 &perim, &count, &dummy, &dummy, &dummy
346 ac, bc, theta, wa, wb, af1, af2, b1, b2,
347 &count, vdrag, &dummy, &a_lblk, &b_lblk, &dummy, &dummy
349 a_lblock(ia,ib) = a_lblk;
355 if (bc >= bf1 && bc < bf2)
360 bc, ac, 0.5*dia, bf1, bf2, a1, a2,
361 &count, &(vdrag.yy()), &dummy
369 &perim, &dummy, &count, &dummy, &dummy
376 ac, bc, theta, wa, wb,
378 &count, vdrag, &dummy, &a_lblk, &b_lblk, &dummy, &dummy
380 b_lblock(ia,ib) = b_lblk;
394 const scalar multiplier
398 const label nBlock =
range.size();
401 scalar totVolume = 0;
403 if (nBlock < 2)
return 0;
413 for (label i1 = 0; i1 < nBlock-1; ++i1)
424 for (label i2 = i1 + 1; i2 < nBlock; ++i2)
431 if (max1.x() <= blk2.
x())
439 || max1.z() <= blk2.
z()
440 || max2.y() <= blk1.
y()
441 || max2.z() <= blk1.
z()
468 over.
typeId = 81 + int(15 * multiplier);
470 if (obsHasMinSize(over.
span, pars))
473 totVolume -= over.
volume();
481 blocks.
append(std::move(newBlocks));
499 const label nBlock =
range.size();
500 const label nCyl = cylinders.
size();
503 scalar totVolume = 0;
505 if (!nBlock || !nCyl)
return 0;
507 scalar
area, a_lblk, b_lblk, dummy, a_centre, b_centre;
521 for (label i1 = 0; i1 < nBlock; i1++)
532 while (i2 < nCyl-1 && cylinders[cylOrder[i2]] < blk1)
537 for (; i2 < nCyl; ++i2)
551 const scalar zm2 = cyl2.
z() + cyl2.
len();
552 if (blk1.
z() > zm2 || cyl2.
z() > max1.z())
continue;
554 if ( cyl2.
dia() == 0.0 )
558 cyl2.
x(), cyl2.
y(), cyl2.
theta(), cyl2.
wa, cyl2.
wb,
561 &dummy, dum2, &dummy, &a_lblk, &b_lblk,
569 cyl2.
x(), cyl2.
y(), 0.5*cyl2.
dia(),
572 &dummy, &dummy, &dummy, &dummy, &dummy
576 cyl2.
x(), cyl2.
y(), 0.5*cyl2.
dia(),
579 &dummy, &dummy, &b_centre
583 cyl2.
y(), cyl2.
x(), 0.5*cyl2.
dia(),
586 &dummy, &dummy, &a_centre
589 if (
equal(area, 0))
continue;
597 const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
599 a_lblk *= blk1.
span.
x() * ratio;
600 b_lblk *= blk1.
span.
y() * ratio;
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());
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);
618 const scalar ym2 = cyl2.
y() + cyl2.
len();
619 if (blk1.
y() > ym2 || cyl2.
y() > max1.y())
continue;
621 if ( cyl2.
dia() == 0.0 )
625 cyl2.
z(), cyl2.
x(), cyl2.
theta(), cyl2.
wa, cyl2.
wb,
628 &dummy, dum2, &dummy, &a_lblk, &b_lblk,
636 cyl2.
z(), cyl2.
x(), 0.5*cyl2.
dia(),
639 &dummy, &dummy, &dummy, &dummy, &dummy
644 cyl2.
z(), cyl2.
x(), 0.5*cyl2.
dia(),
647 &dummy, &dummy, &b_centre
652 cyl2.
x(), cyl2.
z(), 0.5*cyl2.
dia(),
655 &dummy, &dummy, &a_centre
659 if (
equal(area, 0))
continue;
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;
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());
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();
681 const scalar xm2 = cyl2.
x() + cyl2.
len();
682 if (blk1.
x() > xm2 || cyl2.
x() > max1.x())
continue;
684 if ( cyl2.
dia() == 0.0 )
688 cyl2.
y(), cyl2.
z(), cyl2.
theta(), cyl2.
wa, cyl2.
wb,
691 &dummy, dum2, &dummy, &a_lblk, &b_lblk,
699 cyl2.
y(), cyl2.
z(), 0.5*cyl2.
dia(),
702 &dummy, &dummy, &dummy, &dummy, &dummy
707 cyl2.
y(), cyl2.
z(), 0.5*cyl2.
dia(),
710 &dummy, &dummy, &b_centre
715 cyl2.
z(), cyl2.
y(), 0.5*cyl2.
dia(),
718 &dummy, &dummy, &a_centre
723 if (
equal(area, 0))
continue;
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;
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());
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();
746 over.
typeId = PDRobstacle::IGNORE;
751 assert(over.
x() > -10000.0);
753 if (obsHasMinSize(over.
span, pars))
756 totVolume -= over.
volume();
763 blocks.
append(std::move(newBlocks));
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.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Grid locations in an axis direction.
scalar C(const label i) const
Cell centre at element position.
label findCell(const scalar p) const
Find the cell index enclosing this location.
label nCells() const
The number of cells in this direction.
label nPoints() const
The number of points in this direction.
scalar width(const label i) const
Cell size at element position.
const scalar & clip(const scalar &val) const
Obstacle definitions for PDR.
scalar volume() const
Volume of the obstacle.
point pt
The obstacle location.
direction orient
The x/y/z orientation (0,1,2)
scalar x() const
Obstacle position accessors.
vector span
The obstacle dimensions (for boxes)
int typeId
The obstacle type-id.
Parameters for PDRsetFields.
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
A templated (2 x 2) symmetric tensor of objects of <T>, effectively containing 3 elements,...
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
T & first()
Return the first element of the list.
SubList< T > slice(const label pos, label len=-1)
Return SubList slice (non-const access) - no range checking.
void size(const label n)
Older name for setAddressableSize.
T & last()
Return the last element of the list.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
A range or interval of labels defined by a start and a size.
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.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
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)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Cmpt cmptProduct(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)