songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
1339c3357b
openfoam/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacle.C
Mark Olesen 1339c3357b ENH: general boundBox/treeBoundBox improvements 
- null() static method
  * as const reference to the invertedBox with the appropriate casting.

- boundBox inflate(random)
  * refactored from treeBoundBox::extend, but allows in-place modification

- boundBox::hexFaces() instead of boundBox::faces
  * rarely used, but avoids confusion with treeBoundBox::faces
    and reuses hexCell face definitions without code duplication

- boundBox::hexCorners() for corner points corresponding to a hexCell.
  Can also be accessed from a treeBoundBox without ambiguity with
  points(), which could be hex corners (boundBox) or octant corners
  (treeBoundBox)

- boundBox::add with pairs of points
  * convenient (for example) when adding edges or a 'box' that has
    been extracted from a primitive mesh shape.

- declare boundBox nPoints(), nFaces(), nEdges() as per hexCell

ENH: return invertedBox instead of FatalError for empty trees

- similar to #2612

ENH: cellShape(HEX, ...) + boundBox hexCorners for block meshes

STYLE: cellModel::ref(...) instead of de-reference cellModel::ptr(...)
2022-11-24 12:21:01 +00:00

736 lines
18 KiB
C
Raw Blame History

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2019-2022 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 "PDRobstacle.H"
#include "boundBox.H"
#include "meshedSurf.H"
#include "axisAngleRotation.H"
#include "coordinateSystem.H"
#include "foamVtkSurfaceWriter.H"
#include "unitConversion.H"
#include "addToMemberFunctionSelectionTable.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineMemberFunctionSelectionTable(PDRobstacle, read, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::PDRobstacle::PDRobstacle()
:
groupId(0),
typeId(0),
orient(vector::X),
sortBias(0),
pt(Zero),
span(Zero),
wa(0),
wb(0),
vbkge(0),
xbkge(0),
ybkge(0),
zbkge(0),
blowoff_type(0),
identifier()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::PDRobstacle::clear()
{
groupId = 0;
typeId = 0;
orient = vector::X;
sortBias = 0;
pt = Zero;
span = Zero;
wa = 0;
wb = 0;
vbkge = 0;
xbkge = 0;
ybkge = 0;
zbkge = 0;
blowoff_type = 0;
identifier.clear();
}
void Foam::PDRobstacle::readProperties(const dictionary& dict)
{
PDRobstacle::clear();
// Read as word, which handles quoted or unquoted entries
word obsName;
if (dict.readIfPresent("name", obsName))
{
identifier = std::move(obsName);
}
}
void Foam::PDRobstacle::scale(const scalar factor)
{
if (factor <= 0)
{
return;
}
sortBias *= factor;
switch (typeId)
{
case PDRobstacle::CYLINDER:
{
pt *= factor;
dia() *= factor;
len() *= factor;
break;
}
case PDRobstacle::DIAG_BEAM:
{
pt *= factor;
len() *= factor;
wa *= factor;
wb *= factor;
break;
}
case PDRobstacle::CUBOID_1:
case PDRobstacle::LOUVRE_BLOWOFF:
case PDRobstacle::CUBOID:
case PDRobstacle::WALL_BEAM:
case PDRobstacle::GRATING:
case PDRobstacle::RECT_PATCH:
{
pt *= factor;
span *= factor;
if (typeId == PDRobstacle::GRATING)
{
slat_width *= factor;
}
break;
}
}
}
Foam::scalar Foam::PDRobstacle::volume() const
{
scalar vol = 0;
switch (typeId)
{
case PDRobstacle::CYLINDER:
vol = 0.25 * mathematical::pi * sqr(dia()) * len();
break;
case PDRobstacle::DIAG_BEAM:
vol = wa * wb * len();
break;
case PDRobstacle::CUBOID_1:
case PDRobstacle::LOUVRE_BLOWOFF:
case PDRobstacle::CUBOID:
case PDRobstacle::WALL_BEAM:
case PDRobstacle::GRATING:
vol = cmptProduct(span) * vbkge;
break;
}
return vol;
}
bool Foam::PDRobstacle::tooSmall(const scalar minWidth) const
{
if (minWidth <= 0)
{
return false;
}
switch (typeId)
{
case PDRobstacle::CYLINDER:
{
// The effective half-width
if ((0.25 * dia() * sqrt(mathematical::pi)) <= minWidth)
{
return true;
}
break;
}
case PDRobstacle::DIAG_BEAM:
{
if
(
(len() <= minWidth && wa <= minWidth)
|| (len() <= minWidth && wb <= minWidth)
|| (wa <= minWidth && wb <= minWidth)
)
{
return true;
}
break;
}
case PDRobstacle::CUBOID_1:
case PDRobstacle::LOUVRE_BLOWOFF:
case PDRobstacle::CUBOID:
case PDRobstacle::WALL_BEAM:
case PDRobstacle::GRATING:
case PDRobstacle::RECT_PATCH:
{
if
(
(span.x() <= minWidth && span.y() <= minWidth)
|| (span.y() <= minWidth && span.z() <= minWidth)
|| (span.z() <= minWidth && span.x() <= minWidth)
)
{
return true;
}
break;
}
}
return false;
}
Foam::volumeType Foam::PDRobstacle::trim(const boundBox& bb)
{
volumeType::type vt = volumeType::UNKNOWN;
if (!bb.valid() || !typeId)
{
return vt;
}
switch (typeId)
{
case PDRobstacle::CYLINDER:
{
const scalar rad = 0.5*dia();
direction e1 = vector::X;
direction e2 = vector::Y;
direction e3 = vector::Z;
if (orient == vector::X)
{
e1 = vector::Y;
e2 = vector::Z;
e3 = vector::X;
}
else if (orient == vector::Y)
{
e1 = vector::Z;
e2 = vector::X;
e3 = vector::Y;
}
else
{
orient = vector::Z; // extra safety?
}
if
(
(pt[e1] + rad <= bb.min()[e1])
|| (pt[e2] + rad <= bb.min()[e2])
|| (pt[e3] + len() <= bb.min()[e3])
|| (pt[e1] - rad >= bb.max()[e1])
|| (pt[e2] - rad >= bb.max()[e2])
|| (pt[e3] >= bb.max()[e3])
)
{
// No overlap
return volumeType::OUTSIDE;
}
vt = volumeType::INSIDE;
// Trim obstacle length as required
if (pt[e3] < bb.min()[e3])
{
vt = volumeType::MIXED;
len() -= bb.min()[e3] - pt[e3];
pt[e3] = bb.min()[e3];
}
if (pt[e3] + len() > bb.max()[e3])
{
vt = volumeType::MIXED;
len() = bb.max()[e3] - pt[e3];
}
// Cannot trim diameter very well, so just mark as protruding
if
(
(pt[e1] - rad < bb.min()[e1]) || (pt[e1] + rad > bb.max()[e1])
|| (pt[e2] - rad < bb.min()[e2]) || (pt[e2] + rad > bb.max()[e2])
)
{
vt = volumeType::MIXED;
}
break;
}
case PDRobstacle::DIAG_BEAM:
{
// Not implemented
break;
}
case PDRobstacle::CUBOID_1:
case PDRobstacle::LOUVRE_BLOWOFF:
case PDRobstacle::CUBOID:
case PDRobstacle::WALL_BEAM:
case PDRobstacle::GRATING:
case PDRobstacle::RECT_PATCH:
{
for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
{
if
(
((pt[cmpt] + span[cmpt]) < bb.min()[cmpt])
|| (pt[cmpt] > bb.max()[cmpt])
)
{
// No overlap
return volumeType::OUTSIDE;
}
}
vt = volumeType::INSIDE;
// Trim obstacle as required
for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
{
if (pt[cmpt] < bb.min()[cmpt])
{
vt = volumeType::MIXED;
if (span[cmpt] > 0)
{
span[cmpt] -= bb.min()[cmpt] - pt[cmpt];
}
pt[cmpt] = bb.min()[cmpt];
}
if (pt[cmpt] + span[cmpt] > bb.max()[cmpt])
{
vt = volumeType::MIXED;
span[cmpt] -= bb.max()[cmpt] - pt[cmpt];
}
}
break;
}
}
return vt;
}
Foam::meshedSurface Foam::PDRobstacle::surface() const
{
meshedSurface surf;
const PDRobstacle& obs = *this;
switch (obs.typeId)
{
case PDRobstacle::CUBOID_1 :
case PDRobstacle::CUBOID :
{
boundBox box(obs.pt, obs.pt + obs.span);
pointField pts(box.points());
faceList fcs(boundBox::hexFaces());
surf.transfer(pts, fcs);
break;
}
case PDRobstacle::DIAG_BEAM :
{
boundBox box(Zero);
switch (orient)
{
case vector::X:
{
box.min() = vector(0, -0.5*obs.wa, -0.5*obs.wb);
box.max() = vector(obs.len(), 0.5*obs.wa, 0.5*obs.wb);
break;
}
case vector::Y:
{
box.min() = vector(-0.5*obs.wb, 0, -0.5*obs.wa);
box.max() = vector(0.5*obs.wb, obs.len(), 0.5*obs.wa);
break;
}
case vector::Z:
{
box.min() = vector(-0.5*obs.wa, -0.5*obs.wb, 0);
box.max() = vector(0.5*obs.wa, 0.5*obs.wb, obs.len());
break;
}
}
coordinateSystem cs
(
obs.pt,
coordinateRotations::axisAngle
(
vector::components(obs.orient),
obs.theta(),
false
)
);
pointField pts0(box.points());
faceList fcs(boundBox::hexFaces());
pointField pts(cs.globalPosition(pts0));
surf.transfer(pts, fcs);
break;
}
case PDRobstacle::CYLINDER :
{
// Tessellation 12 looks fairly reasonable
constexpr int nDiv = 12;
point org(obs.pt);
direction e1 = vector::X;
direction e2 = vector::Y;
direction e3 = vector::Z;
if (obs.orient == vector::X)
{
e1 = vector::Y;
e2 = vector::Z;
e3 = vector::X;
}
else if (obs.orient == vector::Y)
{
e1 = vector::Z;
e2 = vector::X;
e3 = vector::Y;
}
pointField pts(2*nDiv, org);
faceList fcs(2 + nDiv);
// Origin for back
org[e3] += obs.len();
SubList<point>(pts, nDiv, nDiv) = org;
const scalar radius = 0.5*obs.dia();
for (label i=0; i < nDiv; ++i)
{
const scalar angle = (i * mathematical::twoPi) / nDiv;
const scalar s = radius * sin(angle);
const scalar c = radius * cos(angle);
pts[i][e1] += s;
pts[i][e2] += c;
pts[nDiv+i][e1] += s;
pts[nDiv+i][e2] += c;
}
// Side-faces
for (label facei=0; facei < nDiv; ++facei)
{
face& f = fcs[facei];
f.resize(4);
f[0] = facei;
f[3] = (facei + 1) % nDiv;
f[1] = f[0] + nDiv;
f[2] = f[3] + nDiv;
}
{
// Front face
face& f1 = fcs[nDiv];
f1.resize(nDiv);
f1[0] = 0;
for (label pti=1; pti < nDiv; ++pti)
{
f1[pti] = nDiv-pti;
}
// Back face
labelList& f2 = fcs[nDiv+1];
f2 = identity(nDiv, nDiv);
}
surf.transfer(pts, fcs);
break;
}
case PDRobstacle::RECT_PATCH :
{
pointField pts(4, obs.span);
pts[0] = Zero;
switch (obs.inlet_dirn)
{
case -1:
case 1:
{
for (auto& p : pts)
{
p.x() = 0;
}
pts[1].z() = 0;
pts[3].y() = 0;
break;
}
case -2:
case 2:
{
for (auto& p : pts)
{
p.y() = 0;
}
pts[1].x() = 0;
pts[3].z() = 0;
break;
}
default:
{
for (auto& p : pts)
{
p.z() = 0;
}
pts[1].y() = 0;
pts[3].x() = 0;
break;
}
}
// pts += obs.pt;
faceList fcs(one{}, face(identity(4)));
surf.transfer(pts, fcs);
break;
}
default:
break;
// LOUVRE_BLOWOFF = 5,
// WALL_BEAM = 7,
// GRATING = 8,
// CIRC_PATCH = 12,
// MESH_PLANE = 46,
}
return surf;
}
Foam::label Foam::PDRobstacle::addPieces
(
vtk::surfaceWriter& surfWriter,
const UList<PDRobstacle>& list,
label pieceId
)
{
for (const PDRobstacle& obs : list)
{
meshedSurface surf(obs.surface());
if (!surf.empty())
{
surfWriter.piece(surf.points(), surf.surfFaces());
surfWriter.writeGeometry();
surfWriter.beginCellData(2);
surfWriter.writeUniform("group", label(obs.groupId));
surfWriter.writeUniform("type", label(obs.typeId));
surfWriter.writeUniform("obstacle", pieceId);
++pieceId;
}
}
return pieceId;
}
void Foam::PDRobstacle::generateVtk
(
const fileName& outputDir,
const UList<PDRobstacle>& obslist,
const UList<PDRobstacle>& cyllist
)
{
label pieceId = 0;
vtk::surfaceWriter surfWriter
(
pointField::null(),
faceList::null(),
// vtk::formatType::INLINE_ASCII,
(outputDir / "Obstacles"),
false // serial only
);
pieceId = addPieces(surfWriter, obslist, pieceId);
pieceId = addPieces(surfWriter, cyllist, pieceId);
Info<< "Wrote " << pieceId << " obstacles (VTK) to "
<< outputDir/"Obstacles" << nl;
}
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const InfoProxy<PDRobstacle>& iproxy
)
{
const PDRobstacle& obs = iproxy.t_;
switch (obs.typeId)
{
case PDRobstacle::CUBOID_1 :
case PDRobstacle::CUBOID :
os << "box { point " << obs.pt
<< "; size " << obs.span
<< "; }";
break;
case PDRobstacle::CYLINDER :
os << "cyl { point " << obs.pt
<< "; length " << obs.len() << "; diameter " << obs.dia()
<< "; direction " << vector::componentNames[obs.orient]
<< "; }";
break;
case PDRobstacle::DIAG_BEAM :
os << "diag { point " << obs.pt
<< "; length " << obs.len()
<< "; width (" << obs.wa << ' ' << obs.wb << ')'
<< "; angle " << radToDeg(obs.theta())
<< "; direction " << vector::componentNames[obs.orient]
<< "; }";
break;
case PDRobstacle::WALL_BEAM :
os << "wallbeam { point " << obs.pt
<< " size " << obs.span
<< "; }";
break;
case PDRobstacle::GRATING :
os << "grate { point " << obs.pt
<< "; size " << obs.span
<< "; slats " << obs.slat_width
<< "; }";
break;
case PDRobstacle::LOUVER_BLOWOFF :
os << "louver { point " << obs.pt
<< "; size " << obs.span
<< "; pressure " << paToBar(obs.blowoff_press)
<< "; }";
break;
case PDRobstacle::RECT_PATCH :
os << "patch { " << obs.pt
<< "; size " << obs.span
<< "; name " << obs.identifier
<< "; }";
break;
case PDRobstacle::OLD_INLET :
case PDRobstacle::OLD_BLOWOFF :
case PDRobstacle::IGNITION :
os << "/* ignored: " << obs.typeId << " */";
break;
default:
os << "/* unknown: " << obs.typeId << " */";
break;
}
return os;
}
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
bool Foam::operator<(const PDRobstacle& a, const PDRobstacle& b)
{
return (a.pt.x() + a.sortBias) < (b.pt.x() + b.sortBias);
}
// ************************************************************************* //
Reference in New Issue View Git Blame Copy Permalink