mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-12-28 03:37:59 +00:00
GIT: rename directory PDRsetFields -> PDR
This commit is contained in:
Mark Olesen
committed by
Andrew Heather
Andrew Heather
parent
649b1b1971
commit
e2d7ad5c60
@ -0,0 +1,769 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2016 Shell Research Ltd.
|
||||
Copyright (C) 2019 OpenCFD Ltd.
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "PDRsetFields.H"
|
||||
#include "PDRutilsInternal.H"
|
||||
#include "mathematicalConstants.H"
|
||||
|
||||
#ifndef FULLDEBUG
|
||||
#ifndef NDEBUG
|
||||
#define NDEBUG
|
||||
#endif
|
||||
#endif
|
||||
#include <cassert>
|
||||
|
||||
using namespace Foam::constant;
|
||||
|
||||
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
// A sign-corrected multiply
|
||||
// This is used for porosity of obstacle intersections
|
||||
inline static scalar COMBLK(const scalar a, const scalar b)
|
||||
{
|
||||
if (a < 0)
|
||||
{
|
||||
return -a * b;
|
||||
}
|
||||
|
||||
return a * b;
|
||||
}
|
||||
|
||||
|
||||
// Obstacle satisfies some minimum size checks.
|
||||
// A volume check misses thin plates, so use area.
|
||||
// Thin sheet overlaps can be produced by touching objects
|
||||
// if the obs_extend parameter is > 0.
|
||||
inline static bool obsHasMinSize(const vector& span, const PDRparams& tol)
|
||||
{
|
||||
return
|
||||
(
|
||||
(cmptProduct(span) > tol.min_overlap_vol)
|
||||
&&
|
||||
(
|
||||
(span.x() * span.y() > tol.min_overlap_area)
|
||||
|| (span.y() * span.z() > tol.min_overlap_area)
|
||||
|| (span.z() * span.x() > tol.min_overlap_area)
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
|
||||
void Foam::PDRutils::one_d_overlap
|
||||
(
|
||||
scalar xmin,
|
||||
scalar xmax,
|
||||
const PDRblock::location& grid,
|
||||
List<scalar>& olap,
|
||||
int *cmin, int *cmax,
|
||||
int *cfmin, int *cfmax
|
||||
)
|
||||
{
|
||||
// Looking at one coordinate direction, called x here, for something
|
||||
// that extends from xmin to xmax, calculate how much it overlaps
|
||||
// each cell in this direction. Result returned in 'olap' List is
|
||||
// the proportion of the grid step overlapped, i.e dimensionless.
|
||||
// First and last steps overlapped given by *cmin, *cmax
|
||||
// Ditto for shifted grid given by *cfmin, *cfmax.
|
||||
|
||||
// Initially zero everywhere
|
||||
olap = Zero;
|
||||
|
||||
if (olap.size() < grid.nPoints())
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "The overlap scratch array is too small, has "
|
||||
<< olap.size() << " but needs " << grid.nPoints() << nl
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
|
||||
// No intersection with the box
|
||||
if (xmax <= grid.first() || grid.last() <= xmin)
|
||||
{
|
||||
// Mark as bad range, cannot iterate
|
||||
*cmin = 0;
|
||||
*cmax = -1;
|
||||
|
||||
// Another bad range (cannot iterate) but for extra safety ensure
|
||||
// that (cfmin -> cmin) and (cmax -> cfmax) cannot iterate either
|
||||
*cfmin = 1;
|
||||
*cfmax = -2;
|
||||
return;
|
||||
}
|
||||
|
||||
// Ensure search is within the (point) bounds
|
||||
xmin = grid.clip(xmin);
|
||||
xmax = grid.clip(xmax);
|
||||
|
||||
// The begin/end of the obstacle
|
||||
*cmin = grid.findCell(xmin);
|
||||
*cmax = grid.findCell(xmax);
|
||||
|
||||
for (label ix = *cmin; ix <= *cmax; ++ix)
|
||||
{
|
||||
olap[ix] = 1.0;
|
||||
}
|
||||
|
||||
// Fixup ends
|
||||
if (*cmax == *cmin)
|
||||
{
|
||||
olap[*cmax] = (xmax - xmin) / grid.width(*cmax);
|
||||
}
|
||||
else
|
||||
{
|
||||
if (grid[*cmin] < xmin)
|
||||
{
|
||||
olap[*cmin] = (grid[*cmin+1] - xmin) / grid.width(*cmin);
|
||||
}
|
||||
|
||||
if (xmax < grid[*cmax+1])
|
||||
{
|
||||
olap[*cmax] = (xmax - grid[*cmax]) / grid.width(*cmax);
|
||||
}
|
||||
}
|
||||
assert(olap[*cmax] >= 0.0);
|
||||
|
||||
|
||||
// Is xmin below/above the cell-centre (for virtual staggered-grid) ?
|
||||
*cfmin =
|
||||
(
|
||||
xmin < grid.C(*cmin)
|
||||
? *cmin
|
||||
: Foam::min(*cmin+1, grid.nCells()-1)
|
||||
);
|
||||
|
||||
// Is xmax below/above the cell-centre (for virtual staggered-grid) ?
|
||||
*cfmax =
|
||||
(
|
||||
xmax < grid.C(*cmax)
|
||||
? *cmax
|
||||
: Foam::min(*cmax+1, grid.nCells()-1)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
/**************************************************************************************************/
|
||||
|
||||
void Foam::PDRutils::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
|
||||
)
|
||||
{
|
||||
// We go one over the relevant min/max limits since these values might be
|
||||
// used. If not, they would have been zeroed in one_d_overlap
|
||||
|
||||
amin = Foam::max(0, amin-1);
|
||||
bmin = Foam::max(0, bmin-1);
|
||||
amax = Foam::min(a_olap.size()-1, amax+1);
|
||||
bmax = Foam::min(b_olap.size()-1, bmax+1);
|
||||
|
||||
for (label ia = amin; ia <= amax; ++ia)
|
||||
{
|
||||
for (label ib = bmin; ib <= bmax; ++ib)
|
||||
{
|
||||
ab_olap(ia,ib) = a_olap[ia] * b_olap[ib];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**************************************************************************************************/
|
||||
|
||||
void Foam::PDRutils::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
|
||||
)
|
||||
{
|
||||
/* This routine calculates the proportion of each (two-dimensional) grid cell
|
||||
overlapped by the circle or angled rectangle. Coordinates are labelled a and b.
|
||||
On entry:
|
||||
ac, bc coordinates of centre of circle or rectangle
|
||||
dia diameter of circle (zero for rectangle)
|
||||
theta, wa, wb parameters for rectangle
|
||||
agrid[] locations of grid lines of a-grid
|
||||
amin, amax first and last cells in a-grid overlapped by object
|
||||
(similarly for b)
|
||||
On exit:
|
||||
abolap 2-D array of (proportionate) area blockage by grid cell
|
||||
a_lblock 2-D array of (proportionate) blockage to a-direction flow
|
||||
(This will be area blockage when extruded in the third coordinate).
|
||||
a_count (2-D array)The contribution of this object to the count of obstacles blocking
|
||||
a-direction flow. This is only non-zero if the object is inside the
|
||||
lateral boundaries of the cell. It is large negative if the cell is
|
||||
totally blocked in this direction.
|
||||
(similarly for b)
|
||||
c_drag 2-D array of tensor that will give tensor drag in each cell (when multiplied
|
||||
Cd, cylinder length, and 0.5 rho*U^2) Dimension: L.
|
||||
|
||||
Note that this routine does not zero array elements outside the amin to amax, bmin to bmax area.
|
||||
*/
|
||||
|
||||
scalar count, a_lblk, b_lblk, perim, dummy;
|
||||
|
||||
symmTensor2D vdrag(Zero);
|
||||
|
||||
// Prevent stepping outside of the array when the obstacle is on the
|
||||
// upper boundary
|
||||
|
||||
// Upper limit of inclusive range is nCells-1
|
||||
amin = Foam::max(0, amin);
|
||||
bmin = Foam::max(0, bmin);
|
||||
amax = Foam::min(amax, agrid.nCells()-1);
|
||||
bmax = Foam::min(bmax, bgrid.nCells()-1);
|
||||
|
||||
for (label ia = amin; ia <= amax; ++ia)
|
||||
{
|
||||
// Cell-centred grid
|
||||
const scalar a1 = agrid[ia];
|
||||
const scalar a2 = agrid[ia+1];
|
||||
|
||||
// Left-shifted staggered face grid (-1 addressing is OK)
|
||||
const scalar af1 = agrid.C(ia-1);
|
||||
const scalar af2 = agrid.C(ia);
|
||||
|
||||
for (label ib = bmin; ib <= bmax; ++ib)
|
||||
{
|
||||
// Cell-centred grid
|
||||
const scalar b1 = bgrid[ib];
|
||||
const scalar b2 = bgrid[ib+1];
|
||||
|
||||
// Left-shifted staggered face grid (-1 addressing is OK)
|
||||
const scalar bf1 = bgrid.C(ib-1);
|
||||
const scalar bf2 = bgrid.C(ib);
|
||||
|
||||
// Do the centred cell
|
||||
if ( dia > 0.0 )
|
||||
{
|
||||
ab_olap(ia,ib) = inters_cy
|
||||
(
|
||||
ac, bc, 0.5*dia, a1, a2, b1, b2, &perim,
|
||||
&dummy, &dummy, &b_lblk, &a_lblk
|
||||
);
|
||||
/* The last two arguments of the above call appear to be reversed, but the inters_cy routine returns
|
||||
the amount of overlap in the a and b direcvtions, which are the blockage to the b and a directions. */
|
||||
|
||||
/* abolap * cell area is area of cylinder in this cell. Divide by PI%D^2/4 to get proportion of cylinder in cell
|
||||
For whole cylinger c_drag should be = D, so multiply by D. */
|
||||
|
||||
c_drag(ia,ib).xx() = c_drag(ia,ib).yy() = 4.0 * ab_olap(ia,ib) * (a2 - a1) * (b2 - b1) / dia / mathematical::pi;
|
||||
c_drag(ia,ib).xy() = Zero;
|
||||
c_count(ia,ib) = perim / (mathematical::pi * dia);
|
||||
|
||||
//******?????
|
||||
scalar area = (a2 - a1) * (b2 - b1);
|
||||
scalar rat = dia * dia / area - 1.5;
|
||||
if (rat > 0.0)
|
||||
{
|
||||
scalar da = ac - 0.5 * (a1 + a2);
|
||||
scalar db = bc - 0.5 * (b1 + b2);
|
||||
scalar dc = std::hypot(da, db);
|
||||
scalar rat1 = min(max((dc / sqrt(area) - 0.3) * 1.4, 0), 1);
|
||||
scalar drg0 = c_drag(ia,ib).xx();
|
||||
scalar drg1 = c_drag(ia,ib).yy();
|
||||
scalar drg = std::hypot(drg0, drg1);
|
||||
c_drag(ia,ib).xx() = drg * ( 1.0 - rat1 ) + drg * da*da/dc/dc * rat1;
|
||||
c_drag(ia,ib).yy() = drg * ( 1.0 - rat1 ) + drg * db*db/dc/dc * rat1;
|
||||
c_drag(ia,ib).xy() = drg * da*db/dc/dc *rat1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
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 );
|
||||
c_count(ia,ib) = perim / ( wa + wb ) * 0.5;
|
||||
}
|
||||
ac_lblock(ia,ib) = a_lblk;
|
||||
bc_lblock(ia,ib) = b_lblk;
|
||||
ab_perim(ia,ib) = perim;
|
||||
|
||||
// Do the a-shifted cell
|
||||
if ( dia > 0.0 ) // I.e. a cylinder, not a d.b.
|
||||
{
|
||||
if (ac >= af1 && ac < af2)
|
||||
{
|
||||
// Only want to block one layer of faces
|
||||
a_lblock(ia,ib) = l_blockage
|
||||
(
|
||||
ac, bc, 0.5*dia,
|
||||
af1, af2, b1, b2, &count, &dummy, &dummy
|
||||
);
|
||||
}
|
||||
inters_cy
|
||||
(
|
||||
ac, bc, 0.5*dia,
|
||||
af1, af2, b1, b2,
|
||||
&perim, &count, &dummy, &dummy, &dummy
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
inters_db
|
||||
(
|
||||
ac, bc, theta, wa, wb, af1, af2, b1, b2,
|
||||
&count, vdrag, &dummy, &a_lblk, &b_lblk, &dummy, &dummy
|
||||
);
|
||||
a_lblock(ia,ib) = a_lblk;
|
||||
}
|
||||
|
||||
// Do the b-shifted cell
|
||||
if ( dia > 0.0 )
|
||||
{
|
||||
if (bc >= bf1 && bc < bf2)
|
||||
{
|
||||
// Only want to block one layer of faces
|
||||
b_lblock(ia,ib) = l_blockage
|
||||
(
|
||||
bc, ac, 0.5*dia, bf1, bf2, a1, a2,
|
||||
&count, &(vdrag.yy()), &dummy
|
||||
);
|
||||
}
|
||||
|
||||
inters_cy
|
||||
(
|
||||
ac, bc, 0.5*dia,
|
||||
a1, a2, bf1, bf2,
|
||||
&perim, &dummy, &count, &dummy, &dummy
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
inters_db
|
||||
(
|
||||
ac, bc, theta, wa, wb,
|
||||
a1, a2, bf1, bf2,
|
||||
&count, vdrag, &dummy, &a_lblk, &b_lblk, &dummy, &dummy
|
||||
);
|
||||
b_lblock(ia,ib) = b_lblk;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // End circle_overlap
|
||||
|
||||
|
||||
/**************************************************************************************************/
|
||||
|
||||
scalar block_overlap
|
||||
(
|
||||
DynamicList<PDRobstacle>& blocks,
|
||||
const labelRange& range,
|
||||
const scalar multiplier
|
||||
)
|
||||
{
|
||||
// Size information
|
||||
const label nBlock = range.size();
|
||||
|
||||
// The return value
|
||||
scalar totVolume = 0;
|
||||
|
||||
if (nBlock < 2) return 0;
|
||||
|
||||
|
||||
// Sort blocks by their x-position (with sortBias)
|
||||
labelList blkOrder;
|
||||
sortedOrder(blocks[range], blkOrder);
|
||||
|
||||
DynamicList<PDRobstacle> newBlocks;
|
||||
|
||||
// Work through the sorted blocks
|
||||
for (label i1 = 0; i1 < nBlock-1; ++i1)
|
||||
{
|
||||
const PDRobstacle& blk1 = blocks[range[blkOrder[i1]]];
|
||||
|
||||
// Upper coordinates
|
||||
const vector max1 = blk1.pt + blk1.span;
|
||||
|
||||
// For second block start with the next one on the list, and
|
||||
// stop when we find the first one whose biased x-position
|
||||
// is beyond the end of the block1
|
||||
|
||||
for (label i2 = i1 + 1; i2 < nBlock; ++i2)
|
||||
{
|
||||
const PDRobstacle& blk2 = blocks[range[blkOrder[i2]]];
|
||||
|
||||
// Upper coordinates
|
||||
const vector max2 = blk2.pt + blk2.span;
|
||||
|
||||
if (max1.x() <= blk2.x())
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
if
|
||||
(
|
||||
max1.y() <= blk2.y()
|
||||
|| max1.z() <= blk2.z()
|
||||
|| max2.y() <= blk1.y()
|
||||
|| max2.z() <= blk1.z()
|
||||
|| (blk1.vbkge * blk2.vbkge <= 0)
|
||||
)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
|
||||
{
|
||||
PDRobstacle over;
|
||||
|
||||
over.pt = max(blk1.pt, blk2.pt);
|
||||
over.span = min(max1, max2) - over.pt;
|
||||
|
||||
assert(cmptProduct(over.span) > 0.0);
|
||||
|
||||
// This routine should only have been called for all +ve o r all -ve obstacles
|
||||
assert(blk1.vbkge * blk2.vbkge > 0);
|
||||
/* At the first level of intersection, we create an obstacle of blockage -1 (if both objects solid)
|
||||
to cancel out the double counting. (multiplier is 1).
|
||||
?? COMBLK does a (sign corrected) multiply; is this corrrect for porous obstacles?
|
||||
Depends on how blockages were summed in the first place. In fact this -ve obstacle
|
||||
concept only works if the blockages are summed??*/
|
||||
over.vbkge = - COMBLK( blk1.vbkge, blk2.vbkge ) * multiplier;
|
||||
over.xbkge = - COMBLK( blk1.xbkge, blk2.xbkge ) * multiplier;
|
||||
over.ybkge = - COMBLK( blk1.ybkge, blk2.ybkge ) * multiplier;
|
||||
over.zbkge = - COMBLK( blk1.zbkge, blk2.zbkge ) * multiplier;
|
||||
over.typeId = 81 + int(15 * multiplier); // Not subsequently used
|
||||
|
||||
if (obsHasMinSize(over.span, pars))
|
||||
{
|
||||
// Obstacle satisfies some minimum size checks
|
||||
totVolume -= over.volume();
|
||||
|
||||
newBlocks.append(over);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
blocks.append(std::move(newBlocks));
|
||||
|
||||
return totVolume;
|
||||
}
|
||||
|
||||
|
||||
/**************************************************************************************************/
|
||||
|
||||
using namespace Foam::PDRutils;
|
||||
|
||||
scalar block_cylinder_overlap
|
||||
(
|
||||
DynamicList<PDRobstacle>& blocks,
|
||||
const labelRange& range,
|
||||
const UList<PDRobstacle>& cylinders
|
||||
)
|
||||
{
|
||||
// Size information
|
||||
const label nBlock = range.size();
|
||||
const label nCyl = cylinders.size();
|
||||
|
||||
// The return value
|
||||
scalar totVolume = 0;
|
||||
|
||||
if (!nBlock || !nCyl) return 0;
|
||||
|
||||
scalar area, a_lblk, b_lblk, dummy, a_centre, b_centre;
|
||||
symmTensor2D dum2;
|
||||
|
||||
|
||||
// Sort blocks and cylinders by their x-position (with sortBias)
|
||||
labelList blkOrder;
|
||||
sortedOrder(blocks[range], blkOrder);
|
||||
|
||||
labelList cylOrder;
|
||||
sortedOrder(cylinders, cylOrder);
|
||||
|
||||
DynamicList<PDRobstacle> newBlocks;
|
||||
|
||||
// Work through the sorted blocks
|
||||
for (label i1 = 0; i1 < nBlock; i1++)
|
||||
{
|
||||
const PDRobstacle& blk1 = blocks[range[blkOrder[i1]]];
|
||||
|
||||
// Upper coordinates
|
||||
const vector max1 = blk1.pt + blk1.span;
|
||||
|
||||
// Cyls whose end is before start of this block no longer
|
||||
// need to be considered
|
||||
|
||||
label i2 = 0;
|
||||
while (i2 < nCyl-1 && cylinders[cylOrder[i2]] < blk1)
|
||||
{
|
||||
++i2;
|
||||
}
|
||||
|
||||
for (/*nil*/; i2 < nCyl; ++i2)
|
||||
{
|
||||
const PDRobstacle& cyl2 = cylinders[cylOrder[i2]];
|
||||
|
||||
// Calculate overlap in axis direction; if zero continue.
|
||||
// Calculate 2-d overlap and c 0f g; if area zero continue.
|
||||
|
||||
PDRobstacle over;
|
||||
|
||||
|
||||
switch (cyl2.orient)
|
||||
{
|
||||
case vector::Z:
|
||||
{
|
||||
const scalar zm2 = cyl2.z() + cyl2.len();
|
||||
if (blk1.z() > zm2 || cyl2.z() > max1.z()) continue;
|
||||
|
||||
if ( cyl2.dia() == 0.0 )
|
||||
{
|
||||
area = inters_db
|
||||
(
|
||||
cyl2.x(), cyl2.y(), cyl2.theta(), cyl2.wa, cyl2.wb,
|
||||
blk1.x(), max1.x(),
|
||||
blk1.y(), max1.y(),
|
||||
&dummy, dum2, &dummy, &a_lblk, &b_lblk,
|
||||
&a_centre, &b_centre
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
area = inters_cy
|
||||
(
|
||||
cyl2.x(), cyl2.y(), 0.5*cyl2.dia(),
|
||||
blk1.x(), max1.x(),
|
||||
blk1.y(), max1.y(),
|
||||
&dummy, &dummy, &dummy, &dummy, &dummy
|
||||
);
|
||||
b_lblk = l_blockage
|
||||
(
|
||||
cyl2.x(), cyl2.y(), 0.5*cyl2.dia(),
|
||||
blk1.x(), max1.x(),
|
||||
blk1.y(), max1.y(),
|
||||
&dummy, &dummy, &b_centre
|
||||
);
|
||||
a_lblk = l_blockage
|
||||
(
|
||||
cyl2.y(), cyl2.x(), 0.5*cyl2.dia(),
|
||||
blk1.y(), max1.y(),
|
||||
blk1.x(), max1.x(),
|
||||
&dummy, &dummy, &a_centre
|
||||
);
|
||||
}
|
||||
if (equal(area, 0)) continue;
|
||||
assert(a_lblk >0.0);
|
||||
assert(b_lblk >0.0);
|
||||
|
||||
// The intersection between a circle and a rectangle can be an odd shape.
|
||||
// We have its area. a_lblk and b_lblk are dimensions of enclosing rectangle
|
||||
// and a_centre and b_centre its centre. We scale this rectangle down to
|
||||
// the correct areacorrect area, as a rectangular approximation to the intersection.
|
||||
const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
|
||||
|
||||
a_lblk *= blk1.span.x() * ratio;
|
||||
b_lblk *= blk1.span.y() * ratio;
|
||||
assert(b_lblk >0.0);
|
||||
assert(a_lblk >0.0);
|
||||
|
||||
over.x() = a_centre - 0.5 * a_lblk;
|
||||
over.y() = b_centre - 0.5 * b_lblk;
|
||||
over.z() = max(blk1.z(), cyl2.z());
|
||||
|
||||
over.span.x() = a_lblk;
|
||||
over.span.y() = b_lblk;
|
||||
over.span.z() = min(max1.z(), cyl2.z() + cyl2.len()) - over.z();
|
||||
assert(over.x() > -200.0);
|
||||
assert(over.x() < 2000.0);
|
||||
}
|
||||
break;
|
||||
|
||||
case vector::Y:
|
||||
{
|
||||
const scalar ym2 = cyl2.y() + cyl2.len();
|
||||
if (blk1.y() > ym2 || cyl2.y() > max1.y()) continue;
|
||||
|
||||
if ( cyl2.dia() == 0.0 )
|
||||
{
|
||||
area = inters_db
|
||||
(
|
||||
cyl2.z(), cyl2.x(), cyl2.theta(), cyl2.wa, cyl2.wb,
|
||||
blk1.z(), max1.z(),
|
||||
blk1.x(), max1.x(),
|
||||
&dummy, dum2, &dummy, &a_lblk, &b_lblk,
|
||||
&a_centre, &b_centre
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
area = inters_cy
|
||||
(
|
||||
cyl2.z(), cyl2.x(), 0.5*cyl2.dia(),
|
||||
blk1.z(), max1.z(),
|
||||
blk1.x(), max1.x(),
|
||||
&dummy, &dummy, &dummy, &dummy, &dummy
|
||||
);
|
||||
|
||||
b_lblk = l_blockage
|
||||
(
|
||||
cyl2.z(), cyl2.x(), 0.5*cyl2.dia(),
|
||||
blk1.z(), max1.z(),
|
||||
blk1.x(), max1.x(),
|
||||
&dummy, &dummy, &b_centre
|
||||
);
|
||||
|
||||
a_lblk = l_blockage
|
||||
(
|
||||
cyl2.x(), cyl2.z(), 0.5*cyl2.dia(),
|
||||
blk1.x(), max1.x(),
|
||||
blk1.z(), max1.z(),
|
||||
&dummy, &dummy, &a_centre
|
||||
);
|
||||
}
|
||||
|
||||
if (equal(area, 0)) continue;
|
||||
assert(a_lblk >0.0);
|
||||
assert(b_lblk >0.0);
|
||||
|
||||
// a_lblk and b_lblk are dimensions of enclosing rectangle.
|
||||
// Need to scale to correct area
|
||||
const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
|
||||
a_lblk *= blk1.span.z() * ratio;
|
||||
b_lblk *= blk1.span.x() * ratio;
|
||||
|
||||
over.z() = a_centre - a_lblk * 0.5;
|
||||
over.x() = b_centre - b_lblk * 0.5;
|
||||
over.y() = max(blk1.y(), cyl2.y());
|
||||
|
||||
over.span.z() = a_lblk;
|
||||
over.span.x() = b_lblk;
|
||||
over.span.y() = min(max1.y(), cyl2.y() + cyl2.len()) - over.y();
|
||||
}
|
||||
break;
|
||||
|
||||
case vector::X:
|
||||
{
|
||||
const scalar xm2 = cyl2.x() + cyl2.len();
|
||||
if (blk1.x() > xm2 || cyl2.x() > max1.x()) continue;
|
||||
|
||||
if ( cyl2.dia() == 0.0 )
|
||||
{
|
||||
area = inters_db
|
||||
(
|
||||
cyl2.y(), cyl2.z(), cyl2.theta(), cyl2.wa, cyl2.wb,
|
||||
blk1.y(), max1.y(),
|
||||
blk1.z(), max1.z(),
|
||||
&dummy, dum2, &dummy, &a_lblk, &b_lblk,
|
||||
&a_centre, &b_centre
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
area = inters_cy
|
||||
(
|
||||
cyl2.y(), cyl2.z(), 0.5*cyl2.dia(),
|
||||
blk1.y(), max1.y(),
|
||||
blk1.z(), max1.z(),
|
||||
&dummy, &dummy, &dummy, &dummy, &dummy
|
||||
);
|
||||
|
||||
b_lblk = l_blockage
|
||||
(
|
||||
cyl2.y(), cyl2.z(), 0.5*cyl2.dia(),
|
||||
blk1.y(), max1.y(),
|
||||
blk1.z(), max1.z(),
|
||||
&dummy, &dummy, &b_centre
|
||||
);
|
||||
|
||||
a_lblk = l_blockage
|
||||
(
|
||||
cyl2.z(), cyl2.y(), 0.5*cyl2.dia(),
|
||||
blk1.z(), max1.z(),
|
||||
blk1.y(), max1.y(),
|
||||
&dummy, &dummy, &a_centre
|
||||
);
|
||||
|
||||
}
|
||||
|
||||
if (equal(area, 0)) continue;
|
||||
assert(a_lblk >0.0);
|
||||
assert(b_lblk >0.0);
|
||||
|
||||
// a_lblk and b_lblk are dimensions of enclosing rectangle.
|
||||
// Need to scale to correct area
|
||||
const scalar ratio = std::sqrt( area / a_lblk / b_lblk );
|
||||
assert(ratio >-10000.0);
|
||||
assert(ratio <10000.0);
|
||||
a_lblk *= blk1.span.y() * ratio;
|
||||
b_lblk *= blk1.span.z() * ratio;
|
||||
|
||||
over.y() = a_centre - a_lblk * 0.5;
|
||||
over.z() = b_centre - b_lblk * 0.5;
|
||||
over.x() = max(blk1.x(), cyl2.x());
|
||||
|
||||
over.span.y() = a_lblk;
|
||||
over.span.z() = b_lblk;
|
||||
over.span.x() = min(max1.x(), cyl2.x() + cyl2.len()) - over.x();
|
||||
}
|
||||
break;
|
||||
}
|
||||
over.vbkge = over.xbkge = over.ybkge = over.zbkge = -1.0;
|
||||
over.typeId = PDRobstacle::IGNORE;
|
||||
|
||||
assert(cmptProduct(over.span) > 0.0);
|
||||
assert(b_lblk >0.0);
|
||||
assert(a_lblk >0.0);
|
||||
assert(over.x() > -10000.0);
|
||||
|
||||
if (obsHasMinSize(over.span, pars))
|
||||
{
|
||||
// Obstacle satisfies some minimum size checks
|
||||
totVolume -= over.volume();
|
||||
|
||||
newBlocks.append(over);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
blocks.append(std::move(newBlocks));
|
||||
|
||||
return totVolume;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user