PDRobstacle.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) 2019-2021 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "PDRobstacle.H"
29#include "boundBox.H"
30#include "meshedSurf.H"
31#include "axisAngleRotation.H"
32#include "coordinateSystem.H"
34#include "unitConversion.H"
36
37using namespace Foam::constant;
38
39// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40
41namespace Foam
42{
44}
45
46
47// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48
50:
51 groupId(0),
52 typeId(0),
53 orient(vector::X),
54 sortBias(0),
55 pt(Zero),
56 span(Zero),
57 wa(0),
58 wb(0),
59 vbkge(0),
60 xbkge(0),
61 ybkge(0),
62 zbkge(0),
63 blowoff_type(0),
64 identifier()
65{}
66
67
68// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
69
71{
72 groupId = 0;
73 typeId = 0;
75 sortBias = 0;
76 pt = Zero;
77 span = Zero;
78 wa = 0;
79 wb = 0;
80 vbkge = 0;
81 xbkge = 0;
82 ybkge = 0;
83 zbkge = 0;
84 blowoff_type = 0;
85 identifier.clear();
86}
87
88
90{
92
93 // Read as word, which handles quoted or unquoted entries
94 word obsName;
95
96 if (dict.readIfPresent("name", obsName))
97 {
98 identifier = std::move(obsName);
99 }
100}
101
102
103void Foam::PDRobstacle::scale(const scalar factor)
104{
105 if (factor <= 0)
106 {
107 return;
108 }
109
110 sortBias *= factor;
111
112 switch (typeId)
113 {
115 {
116 pt *= factor;
117
118 dia() *= factor;
119 len() *= factor;
120 break;
121 }
122
124 {
125 pt *= factor;
126
127 len() *= factor;
128 wa *= factor;
129 wb *= factor;
130 break;
131 }
132
139 {
140 pt *= factor;
141 span *= factor;
142
143 if (typeId == PDRobstacle::GRATING)
144 {
145 slat_width *= factor;
146 }
147 break;
148 }
149 }
150}
151
152
153Foam::scalar Foam::PDRobstacle::volume() const
154{
155 scalar vol = 0;
156
157 switch (typeId)
158 {
160 vol = 0.25 * mathematical::pi * sqr(dia()) * len();
161 break;
162
164 vol = wa * wb * len();
165 break;
166
172 vol = cmptProduct(span) * vbkge;
173 break;
174 }
175
176 return vol;
177}
178
179
180bool Foam::PDRobstacle::tooSmall(const scalar minWidth) const
181{
182 if (minWidth <= 0)
183 {
184 return false;
185 }
186
187 switch (typeId)
188 {
190 {
191 // The effective half-width
192 if ((0.25 * dia() * sqrt(mathematical::pi)) <= minWidth)
193 {
194 return true;
195 }
196 break;
197 }
198
200 {
201 if
202 (
203 (len() <= minWidth && wa <= minWidth)
204 || (len() <= minWidth && wb <= minWidth)
205 || (wa <= minWidth && wb <= minWidth)
206 )
207 {
208 return true;
209 }
210 break;
211 }
212
219 {
220 if
221 (
222 (span.x() <= minWidth && span.y() <= minWidth)
223 || (span.y() <= minWidth && span.z() <= minWidth)
224 || (span.z() <= minWidth && span.x() <= minWidth)
225 )
226 {
227 return true;
228 }
229
230 break;
231 }
232 }
233
234 return false;
235}
236
237
239{
241
242 if (!bb.valid() || !typeId)
243 {
244 return vt;
245 }
246
247 switch (typeId)
248 {
250 {
251 const scalar rad = 0.5*dia();
252
253 direction e1 = vector::X;
254 direction e2 = vector::Y;
255 direction e3 = vector::Z;
256
257 if (orient == vector::X)
258 {
259 e1 = vector::Y;
260 e2 = vector::Z;
261 e3 = vector::X;
262 }
263 else if (orient == vector::Y)
264 {
265 e1 = vector::Z;
266 e2 = vector::X;
267 e3 = vector::Y;
268 }
269 else
270 {
271 orient = vector::Z; // extra safety?
272 }
273
274 if
275 (
276 (pt[e1] + rad <= bb.min()[e1])
277 || (pt[e2] + rad <= bb.min()[e2])
278 || (pt[e3] + len() <= bb.min()[e3])
279 || (pt[e1] - rad >= bb.max()[e1])
280 || (pt[e2] - rad >= bb.max()[e2])
281 || (pt[e3] >= bb.max()[e3])
282 )
283 {
284 // No overlap
285 return volumeType::OUTSIDE;
286 }
287
289
290 // Trim obstacle length as required
291 if (pt[e3] < bb.min()[e3])
292 {
294 len() -= bb.min()[e3] - pt[e3];
295 pt[e3] = bb.min()[e3];
296 }
297
298 if (pt[e3] + len() > bb.max()[e3])
299 {
301 len() = bb.max()[e3] - pt[e3];
302 }
303
304 // Cannot trim diameter very well, so just mark as protruding
305 if
306 (
307 (pt[e1] - rad < bb.min()[e1]) || (pt[e1] + rad > bb.max()[e1])
308 || (pt[e2] - rad < bb.min()[e2]) || (pt[e2] + rad > bb.max()[e2])
309 )
310 {
312 }
313
314 break;
315 }
316
317
319 {
320 // Not implemented
321 break;
322 }
323
324
331 {
332 for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
333 {
334 if
335 (
336 ((pt[cmpt] + span[cmpt]) < bb.min()[cmpt])
337 || (pt[cmpt] > bb.max()[cmpt])
338 )
339 {
340 // No overlap
341 return volumeType::OUTSIDE;
342 }
343 }
344
345
347
348 // Trim obstacle as required
349
350 for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
351 {
352 if (pt[cmpt] < bb.min()[cmpt])
353 {
355 if (span[cmpt] > 0)
356 {
357 span[cmpt] -= bb.min()[cmpt] - pt[cmpt];
358 }
359 pt[cmpt] = bb.min()[cmpt];
360 }
361
362
363 if (pt[cmpt] + span[cmpt] > bb.max()[cmpt])
364 {
366 span[cmpt] -= bb.max()[cmpt] - pt[cmpt];
367 }
368 }
369
370 break;
371 }
372 }
373
374 return vt;
375}
376
377
379{
380 meshedSurface surf;
381
382 const PDRobstacle& obs = *this;
383
384 switch (obs.typeId)
385 {
388 {
389 boundBox box(obs.pt, obs.pt + obs.span);
390
391 pointField pts(box.points());
393
394 surf.transfer(pts, fcs);
395
396 break;
397 }
398
400 {
401 boundBox box(Zero);
402
403 switch (orient)
404 {
405 case vector::X:
406 {
407 box.min() = vector(0, -0.5*obs.wa, -0.5*obs.wb);
408 box.max() = vector(obs.len(), 0.5*obs.wa, 0.5*obs.wb);
409 break;
410 }
411
412 case vector::Y:
413 {
414 box.min() = vector(-0.5*obs.wb, 0, -0.5*obs.wa);
415 box.max() = vector(0.5*obs.wb, obs.len(), 0.5*obs.wa);
416 break;
417 }
418
419 case vector::Z:
420 {
421 box.min() = vector(-0.5*obs.wa, -0.5*obs.wb, 0);
422 box.max() = vector(0.5*obs.wa, 0.5*obs.wb, obs.len());
423 break;
424 }
425 }
426
428 (
429 obs.pt,
431 (
433 obs.theta(),
434 false
435 )
436 );
437
438 pointField pts0(box.points());
440
441 pointField pts(cs.globalPosition(pts0));
442
443 surf.transfer(pts, fcs);
444
445 break;
446 }
447
449 {
450 // Tessellation 12 looks fairly reasonable
451
452 constexpr int nDiv = 12;
453
454 point org(obs.pt);
455
456 direction e1 = vector::X;
457 direction e2 = vector::Y;
458 direction e3 = vector::Z;
459
460 if (obs.orient == vector::X)
461 {
462 e1 = vector::Y;
463 e2 = vector::Z;
464 e3 = vector::X;
465 }
466 else if (obs.orient == vector::Y)
467 {
468 e1 = vector::Z;
469 e2 = vector::X;
470 e3 = vector::Y;
471 }
472
473 pointField pts(2*nDiv, org);
474 faceList fcs(2 + nDiv);
475
476 // Origin for back
477 org[e3] += obs.len();
478 SubList<point>(pts, nDiv, nDiv) = org;
479
480 const scalar radius = 0.5*obs.dia();
481
482 for (label i=0; i < nDiv; ++i)
483 {
484 const scalar angle = (i * mathematical::twoPi) / nDiv;
485 const scalar s = radius * sin(angle);
486 const scalar c = radius * cos(angle);
487
488 pts[i][e1] += s;
489 pts[i][e2] += c;
490
491 pts[nDiv+i][e1] += s;
492 pts[nDiv+i][e2] += c;
493 }
494
495 // Side-faces
496 for (label facei=0; facei < nDiv; ++facei)
497 {
498 face& f = fcs[facei];
499 f.resize(4);
500
501 f[0] = facei;
502 f[3] = (facei + 1) % nDiv;
503 f[1] = f[0] + nDiv;
504 f[2] = f[3] + nDiv;
505 }
506
507 {
508 // Front face
509 face& f1 = fcs[nDiv];
510 f1.resize(nDiv);
511
512 f1[0] = 0;
513 for (label pti=1; pti < nDiv; ++pti)
514 {
515 f1[pti] = nDiv-pti;
516 }
517
518 // Back face
519 labelList& f2 = fcs[nDiv+1];
520 f2 = identity(nDiv, nDiv);
521 }
522
523 surf.transfer(pts, fcs);
524
525 break;
526 }
527
529 {
530 pointField pts(4, obs.span);
531 pts[0] = Zero;
532
533 switch (obs.inlet_dirn)
534 {
535 case -1:
536 case 1:
537 {
538 for (auto& p : pts)
539 {
540 p.x() = 0;
541 }
542
543 pts[1].z() = 0;
544 pts[3].y() = 0;
545 break;
546 }
547 case -2:
548 case 2:
549 {
550 for (auto& p : pts)
551 {
552 p.y() = 0;
553 }
554
555 pts[1].x() = 0;
556 pts[3].z() = 0;
557 break;
558 }
559 default:
560 {
561 for (auto& p : pts)
562 {
563 p.z() = 0;
564 }
565
566 pts[1].y() = 0;
567 pts[3].x() = 0;
568 break;
569 }
570 }
571
572 // pts += obs.pt;
573
574 faceList fcs(one{}, face(identity(4)));
575
576 surf.transfer(pts, fcs);
577
578 break;
579 }
580
581 default:
582 break;
583
584// LOUVRE_BLOWOFF = 5,
585// WALL_BEAM = 7,
586// GRATING = 8,
587// CIRC_PATCH = 12,
588// MESH_PLANE = 46,
589 }
590
591 return surf;
592}
593
594
596(
597 vtk::surfaceWriter& surfWriter,
598 const UList<PDRobstacle>& list,
599 label pieceId
600)
601{
602 for (const PDRobstacle& obs : list)
603 {
604 meshedSurface surf(obs.surface());
605
606 if (!surf.empty())
607 {
608 surfWriter.piece(surf.points(), surf.surfFaces());
609
610 surfWriter.writeGeometry();
611 surfWriter.beginCellData(2);
612 surfWriter.writeUniform("group", label(obs.groupId));
613 surfWriter.writeUniform("type", label(obs.typeId));
614 surfWriter.writeUniform("obstacle", pieceId);
615 ++pieceId;
616 }
617 }
618
619 return pieceId;
620}
621
622
624(
625 const fileName& outputDir,
626 const UList<PDRobstacle>& obslist,
627 const UList<PDRobstacle>& cyllist
628)
629{
630 label pieceId = 0;
631
632 vtk::surfaceWriter surfWriter
633 (
636 // vtk::formatType::INLINE_ASCII,
637 (outputDir / "Obstacles"),
638 false // serial only
639 );
640
641 pieceId = addPieces(surfWriter, obslist, pieceId);
642 pieceId = addPieces(surfWriter, cyllist, pieceId);
643
644 Info<< "Wrote " << pieceId << " obstacles (VTK) to "
645 << outputDir/"Obstacles" << nl;
646}
647
648
649// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
650
651Foam::Ostream& Foam::operator<<
652(
653 Ostream& os,
654 const InfoProxy<PDRobstacle>& iproxy
655)
656{
657 const PDRobstacle& obs = iproxy.t_;
658
659 switch (obs.typeId)
660 {
663 os << "box { point " << obs.pt
664 << "; size " << obs.span
665 << "; }";
666 break;
667
669 os << "cyl { point " << obs.pt
670 << "; length " << obs.len() << "; diameter " << obs.dia()
671 << "; direction " << vector::componentNames[obs.orient]
672 << "; }";
673 break;
674
676 os << "diag { point " << obs.pt
677 << "; length " << obs.len()
678 << "; width (" << obs.wa << ' ' << obs.wb << ')'
679 << "; angle " << radToDeg(obs.theta())
680 << "; direction " << vector::componentNames[obs.orient]
681 << "; }";
682 break;
683
685 os << "wallbeam { point " << obs.pt
686 << " size " << obs.span
687 << "; }";
688 break;
689
691 os << "grate { point " << obs.pt
692 << "; size " << obs.span
693 << "; slats " << obs.slat_width
694 << "; }";
695 break;
696
698 os << "louver { point " << obs.pt
699 << "; size " << obs.span
700 << "; pressure " << paToBar(obs.blowoff_press)
701 << "; }";
702 break;
703
705 os << "patch { " << obs.pt
706 << "; size " << obs.span
707 << "; name " << obs.identifier
708 << "; }";
709 break;
710
714 os << "/* ignored: " << obs.typeId << " */";
715 break;
716
717 default:
718 os << "/* unknown: " << obs.typeId << " */";
719 break;
720
721 }
722
723 return os;
724}
725
726
727// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
728
729bool Foam::operator<(const PDRobstacle& a, const PDRobstacle& b)
730{
731 return (a.pt.x() + a.sortBias) < (b.pt.x() + b.sortBias);
732}
733
734
735// ************************************************************************* //
Macros for easy insertion into member function selection tables.
radiation::radiationModel & rad
A helper class for outputting values to Ostream.
Definition: InfoProxy.H:52
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
void transfer(pointField &pointLst, List< Face > &faceLst)
Transfer the components.
const List< Face > & surfFaces() const
Return const access to the faces.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
Obstacle definitions for PDR.
Definition: PDRobstacle.H:75
scalar volume() const
Volume of the obstacle.
static label addPieces(vtk::surfaceWriter &surfWriter, const UList< PDRobstacle > &list, label pieceId=0)
Add pieces to vtp output.
scalar sortBias
Bias for position sorting.
Definition: PDRobstacle.H:119
point pt
The obstacle location.
Definition: PDRobstacle.H:123
direction orient
The x/y/z orientation (0,1,2)
Definition: PDRobstacle.H:116
scalar theta() const
Definition: PDRobstacle.H:131
scalar len() const
Definition: PDRobstacle.H:132
static void generateVtk(const fileName &outputDir, const UList< PDRobstacle > &obslist, const UList< PDRobstacle > &cyllist)
Generate multi-piece VTK (vtp) file of obstacles.
bool tooSmall(const scalar minWidth) const
True if the obstacle is considered to be too small.
meshedSurface surface() const
Surface (points, faces) representation.
label groupId
The group-id.
Definition: PDRobstacle.H:110
void readProperties(const dictionary &dict)
Read the 'name' identifier if present.
void clear()
Reset to a zero obstacle.
scalar dia() const
Definition: PDRobstacle.H:130
vector span
The obstacle dimensions (for boxes)
Definition: PDRobstacle.H:126
@ OLD_BLOWOFF
ignored (old)
Definition: PDRobstacle.H:90
@ IGNITION
ignored (old)
Definition: PDRobstacle.H:94
@ OLD_INLET
ignored (old)
Definition: PDRobstacle.H:89
int typeId
The obstacle type-id.
Definition: PDRobstacle.H:113
PDRobstacle()
Construct zero-initialized.
scalar scale
Overall scale factor.
Definition: PDRparams.H:140
const Field< point_type > & points() const noexcept
Return reference to global points.
A List obtained as a section of another List.
Definition: SubList.H:70
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
components
Component labeling enumeration.
Definition: Vector.H:81
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
bool valid() const
Bounding box is non-inverted.
Definition: boundBoxI.H:76
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
static const faceList faces
Faces to point addressing, as per a 'hex' cell.
Definition: boundBox.H:89
A coordinateRotation specified by a rotation axis and a rotation angle about that axis.
Base class for coordinate system specification, the default coordinate system type is cartesian .
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
void trim()
Inplace trim leading and trailing whitespace.
Definition: exprString.C:96
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
static const char *const componentNames[]
Definition: bool.H:104
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition: volumeType.H:61
type
Volume classification types.
Definition: volumeType.H:66
@ OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
@ MIXED
A location that is partly inside and outside.
Definition: volumeType.H:70
@ UNKNOWN
Unknown state.
Definition: volumeType.H:67
@ INSIDE
A location inside the volume.
Definition: volumeType.H:68
virtual bool beginCellData(label nFields=0)
Begin CellData output section for specified number of fields.
Write faces/points (optionally with fields) as a vtp file or a legacy vtk file.
void writeUniform(const word &fieldName, const Type &val)
Write a uniform field of Cell (Poly) or Point values.
void piece(const pointField &points, const faceList &faces)
Reset point/face references to begin a new piece.
virtual bool writeGeometry()
Write patch topology.
A class for handling words, derived from Foam::string.
Definition: word.H:68
volScalarField & p
OBJstream os(runTime.globalPath()/outputName)
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define defineMemberFunctionSelectionTable(baseType, funcName, argNames)
Define run-time selection table.
const dimensionedScalar c
Speed of light in a vacuum.
Different types of constants.
Namespace for OpenFOAM.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
bool read(const char *buf, int32_t &val)
Same as readInt32.
Definition: int32.H:108
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sin(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Cmpt cmptProduct(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:598
dimensionedScalar sqrt(const dimensionedScalar &ds)
constexpr scalar radToDeg() noexcept
Multiplication factor for radians to degrees conversion.
uint8_t direction
Definition: direction.H:56
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
bool operator<(const IOstreamOption::versionNumber &a, const IOstreamOption::versionNumber &b) noexcept
Version A older than B.
constexpr scalar paToBar(const scalar pa) noexcept
Conversion from Pa to bar.
Vector< scalar > vector
Definition: vector.H:61
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
dictionary dict
volScalarField & b
Definition: createFields.H:27
Unit conversion functions.