PDRobstacleLegacyRead.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-2020 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 "PDRobstacle.H"
31#include "volumeType.H"
32
33using namespace Foam;
34using namespace Foam::constant;
35
36#undef USE_ZERO_INSTANCE_GROUPS
37// #define USE_ZERO_INSTANCE_GROUPS
38
39
40// Counting
42(
43 const fileName& obsFileDir,
44 const wordList& obsFileNames,
46)
47{
48 // Default group with single instance and position (0,0,0)
49 groups(0).clear();
50 groups(0).append(point::zero);
51
52 string buffer;
53
54 if (!obsFileNames.empty())
55 {
56 Info<< "Counting groups in obstacle files" << nl;
57 }
58 for (const word& inputFile : obsFileNames)
59 {
60 Info<< " file: " << inputFile << nl;
61
62 fileName path = (obsFileDir / inputFile);
63
64 IFstream is(path);
65
66 while (is.good())
67 {
68 // Process each line of obstacle files
69 is.getLine(buffer);
70
71 const auto firstCh = buffer.find_first_not_of(" \t\n\v\f\r");
72
73 if
74 (
75 firstCh == std::string::npos
76 || buffer[firstCh] == '#'
77 )
78 {
79 // Empty line or comment line
80 continue;
81 }
82
83 int typeId;
84 double x, y, z; // Double (not scalar) to match sscanf spec
85
86 if
87 (
88 sscanf(buffer.c_str(), "%d %lf %lf %lf", &typeId, &x, &y, &z)<4
89 || typeId == 0
90 || typeId == PDRobstacle::MESH_PLANE
91 )
92 {
93 continue;
94 }
95
96 x *= pars.scale;
97 y *= pars.scale;
98 z *= pars.scale;
99
100 const label groupId = typeId / 100;
101 typeId %= 100;
102
103 if (typeId == PDRobstacle::OLD_INLET)
104 {
105 Info<< "Ignored old-inlet type" << nl;
106 continue;
107 }
108 else if (typeId == PDRobstacle::GRATING && pars.ignoreGratings)
109 {
110 Info<< "Ignored grating" << nl;
111 continue;
112 }
113
114 if (typeId == 0)
115 {
116 // Defining a group location
117 groups(groupId).append(x, y, z);
118 }
119 else if (PDRobstacle::isCylinder(typeId))
120 {
121 // Increment cylinder count for the group
122 groups(groupId).addCylinder();
123 }
124 else
125 {
126 // Increment obstacle count for the group
127 groups(groupId).addObstacle();
128 }
129 }
130 }
131
132
133 label nTotalObs = 0;
134 label nTotalCyl = 0;
135
136 label nMissedObs = 0;
137 label nMissedCyl = 0;
138
139 forAllConstIters(groups, iter)
140 {
141 const auto& group = iter.val();
142
143 nTotalObs += group.nTotalObstacle();
144 nTotalCyl += group.nTotalCylinder();
145
146 if (group.empty())
147 {
148 nMissedObs += group.nObstacle();
149 nMissedCyl += group.nCylinder();
150 }
151 }
152
153 for (const label groupId : groups.sortedToc())
154 {
155 const auto& group = groups[groupId];
156
157 if (groupId)
158 {
159 if (group.size())
160 {
161 Info<< "Found " << group.size()
162 << " instances of group " << groupId << " ("
163 << group.nObstacle() << " obstacles "
164 << group.nCylinder() << " cylinders)"
165 << nl;
166 }
167 }
168 else
169 {
170 // The group 0 is for ungrouped obstacles
171 Info<< "Found "
172 << group.nObstacle() << " obstacles "
173 << group.nCylinder() << " cylinders not in groups" << nl;
174 }
175 }
176
177 Info<< "Number of obstacles: "
178 << (nTotalObs + nTotalCyl) << " ("
179 << nTotalCyl << " cylinders)" << nl;
180
181 if (nMissedObs + nMissedCyl)
182 {
183 #ifdef USE_ZERO_INSTANCE_GROUPS
184
185 nTotalObs += nMissedObs;
186 nTotalCyl += nMissedCyl;
187 Info<< "Adding " << (nMissedObs + nMissedCyl)
188 << " obstacles in groups without instances to default group" << nl;
189
190 #else
191
192 Warning
193 << nl << "Found " << (nMissedObs + nMissedCyl)
194 << " obstacles in groups without instances" << nl << nl;
195
196 if (pars.debugLevel > 1)
197 {
198 for (const label groupId : groups.sortedToc())
199 {
200 const auto& group = groups[groupId];
201
202 if
203 (
204 groupId && group.empty()
205 && (group.nObstacle() || group.nCylinder())
206 )
207 {
208 Info<< " Group " << groupId << " ("
209 << group.nObstacle() << " obstacles "
210 << group.nCylinder() << " cylinders)"
211 << nl;
212 }
213 }
214 }
215 #endif
216 }
217
218 return labelPair(nTotalObs, nTotalCyl);
219}
220
221
223(
224 const fileName& obsFileDir,
225 const wordList& obsFileNames,
226 const Map<obstacleGrouping>& groups,
227 const boundBox& meshBb,
228
230 DynamicList<PDRobstacle>& cylinders
231)
232{
233 // Catch programming errors
234 if (!groups.found(0))
235 {
237 << "No default group 0 defined!" << nl
238 << exit(FatalError);
239 }
240
241 scalar totVolume = 0;
242 label nOutside = 0;
243 label nProtruding = 0;
244
245 scalar shift = pars.obs_expand;
246
247 string buffer;
248
249 if (!obsFileNames.empty())
250 {
251 Info<< "Reading obstacle files" << nl;
252 }
253
254 for (const word& inputFile : obsFileNames)
255 {
256 Info<< " file: " << inputFile << nl;
257
258 fileName path = (obsFileDir / inputFile);
259
260 IFstream is(path);
261
262 label lineNo = 0;
263 while (is.good())
264 {
265 // Process each line of obstacle files
266 ++lineNo;
267 is.getLine(buffer);
268
269 const auto firstCh = buffer.find_first_not_of(" \t\n\v\f\r");
270
271 if
272 (
273 firstCh == std::string::npos
274 || buffer[firstCh] == '#'
275 )
276 {
277 // Empty line or comment line
278 continue;
279 }
280
281 // Quick reject
282
283 int typeId; // Int (not label) to match sscanf spec
284 double x, y, z; // Double (not scalar) to match sscanf spec
285
286 if
287 (
288 sscanf(buffer.c_str(), "%d %lf %lf %lf", &typeId, &x, &y, &z) < 4
289 || typeId == 0
290 || typeId == PDRobstacle::MESH_PLANE
291 )
292 {
293 continue;
294 }
295
296 int groupId = typeId / 100;
297 typeId %= 100;
298
299 if
300 (
301 typeId == PDRobstacle::OLD_INLET
302 || (typeId == PDRobstacle::GRATING && pars.ignoreGratings)
303 )
304 {
305 // Silent - already warned during counting
306 continue;
307 }
308
309 if (typeId == 0)
310 {
311 // Group location - not an obstacle
312 continue;
313 }
314
315 if (!groups.found(groupId))
316 {
317 // Catch programming errors.
318 // - group should be there after the previous read
319 Warning
320 << "Encountered undefined group: " << groupId << nl;
321 continue;
322 }
323
324 #ifdef USE_ZERO_INSTANCE_GROUPS
325 const obstacleGrouping& group =
326 (
327 groups[groups[groupId].size() ? groupId : 0]
328 );
329 #else
330 const obstacleGrouping& group = groups[groupId];
331 #endif
332
333 // Add the obstacle to the list with different position
334 // offsets according to its group. Non-group obstacles
335 // are treated as group 0, which has a single instance
336 // with position (0,0,0) and are added only once.
337
338 PDRobstacle scanObs;
339
340 if
341 (
342 !scanObs.setFromLegacy
343 (
344 (groupId * 100) + typeId,
345 buffer,
346 lineNo,
347 inputFile
348 )
349 )
350 {
351 continue;
352 }
353
354 scanObs.scale(pars.scale);
355
356 // Ignore anything below minimum width
357 if (scanObs.tooSmall(pars.min_width))
358 {
359 continue;
360 }
361
362
363 for (const point& origin : group)
364 {
365 // A different (very small) shift for each obstacle
366 // so that faces cannot be coincident
367 shift += floatSMALL;
368 const scalar shift2 = shift * 2.0;
369
370 switch (typeId)
371 {
372 case PDRobstacle::CYLINDER:
373 {
374 // Make a copy
375 PDRobstacle obs(scanObs);
376
377 // Offset for the position of the group
378 obs.pt += origin;
379
380 // Shift the end outwards so, if exactly on
381 // cell boundary, now overlap cell.
382 // So included in Aw.
383 obs.pt -= point::uniform(shift);
384 obs.len() += shift2;
385 obs.dia() -= floatSMALL;
386
387
388 // Trim against the mesh bounds.
389 // Ignore if it doesn't overlap, or bounds are invalid
390 const volumeType vt = obs.trim(meshBb);
391
392 switch (vt)
393 {
394 case volumeType::OUTSIDE:
395 ++nOutside;
396 continue; // Can ignore the rest
397 break;
398
399 case volumeType::MIXED:
400 ++nProtruding;
401 break;
402
403 default:
404 break;
405 }
406
407 // Later for position sorting
408 switch (obs.orient)
409 {
410 case vector::X:
411 obs.sortBias = obs.len();
412 break;
413 case vector::Y:
414 obs.sortBias = 0.5*obs.dia();
415 break;
416 case vector::Z:
417 obs.sortBias = 0.5*obs.dia();
418 break;
419 }
420
421 totVolume += obs.volume();
422 cylinders.append(obs);
423
424 break;
425 }
426
427 case PDRobstacle::DIAG_BEAM:
428 {
429 // Make a copy
430 PDRobstacle obs(scanObs);
431
432 // Offset for the position of the group
433 obs.pt += origin;
434
435 // Shift the end outwards so, if exactly on
436 // cell boundary, now overlap cell.
437 // So included in Aw.
438 obs.pt -= point::uniform(shift);
439 obs.len() += shift2;
440 obs.wa += shift2;
441 obs.wb += shift2;
442
443 totVolume += obs.volume();
444 cylinders.append(obs);
445
446 break;
447 }
448
449 case PDRobstacle::CUBOID_1: // Cuboid "Type 1"
450 case PDRobstacle::LOUVRE_BLOWOFF: // Louvred wall or blow-off panel
451 case PDRobstacle::CUBOID: // Cuboid
452 case PDRobstacle::WALL_BEAM: // Beam against wall (treated here as normal cuboid)
453 case PDRobstacle::GRATING: // Grating
454 case PDRobstacle::RECT_PATCH: // Inlet, outlet or ather b.c. (rectangular)
455 {
456 // Make a copy
457 PDRobstacle obs(scanObs);
458
459 // Offset for the position of the group
460 obs.pt += origin;
461
462 if (typeId == PDRobstacle::GRATING)
463 {
464 if (obs.slat_width <= 0)
465 {
466 obs.slat_width = pars.def_grating_slat_w;
467 }
468 }
469
470 // Shift the end outwards so, if exactly on
471 // cell boundary, now overlap cell.
472 // So included in Aw.
473 obs.pt -= point::uniform(shift);
474 obs.span += point::uniform(shift2);
475
476
477 // Trim against the mesh bounds.
478 // Ignore if it doesn't overlap, or bounds are invalid
479 const volumeType vt = obs.trim(meshBb);
480
481 switch (vt)
482 {
483 case volumeType::OUTSIDE:
484 ++nOutside;
485 continue; // Can ignore the rest
486 break;
487
488 case volumeType::MIXED:
489 ++nProtruding;
490 break;
491
492 default:
493 break;
494 }
495
496 totVolume += obs.volume();
497
498 blocks.append(obs);
499
500 break;
501 }
502 }
503 }
504 }
505
506 if (nOutside || nProtruding)
507 {
508 Info<< "Warning: " << nOutside << " obstacles outside domain, "
509 << nProtruding << " obstacles partly outside domain" << nl;
510 }
511 }
512
513
514 // #ifdef FULLDEBUG
515 // Info<< blocks << nl << cylinders << nl;
516 // #endif
517
518 return totVolume;
519}
520
521
523(
524 const fileName& obsFileDir,
525 const wordList& obsFileNames,
526 const boundBox& meshBb,
528 DynamicList<PDRobstacle>& cylinders
529)
530{
531 // Still just with legacy reading
532
533 // Count the obstacles and get the group locations
535
536 const labelPair obsCounts =
537 PDRlegacy::readObstacleFiles(obsFileDir, obsFileNames, groups);
538
539 const label nObstacle = obsCounts.first();
540 const label nCylinder = obsCounts.second();
541
542 // Info<< "grouping: " << groups << endl;
543
544 if (!nObstacle && !nCylinder)
545 {
547 << "No obstacles in domain" << nl
548 << exit(FatalError);
549 }
550
551 blocks.clear();
552 blocks.reserve(4 * max(4, nObstacle));
553
554 cylinders.clear();
555 cylinders.reserve(4 * max(4, nCylinder));
556
558 (
559 obsFileDir, obsFileNames, groups,
560 meshBb,
561 blocks,
562 cylinders
563 );
564}
565
566
567// ************************************************************************* //
scalar y
Preparation of fields for PDRFoam.
#define floatSMALL
Definition: PDRsetFields.H:54
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:391
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
void reserve(const label len)
Definition: DynamicListI.H:333
T & first() noexcept
The first element of the list, position [0].
Definition: FixedListI.H:207
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
Definition: HashTable.C:137
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
void clear()
Clear all entries from table.
Definition: HashTable.C:678
Input from file stream, using an ISstream.
Definition: IFstream.H:57
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
Obstacle definitions for PDR.
Definition: PDRobstacle.H:75
bool setFromLegacy(const int groupTypeId, const string &buffer, const label lineNo=-1, const word &inputFile=word::null)
Set values from single-line, multi-column format.
bool tooSmall(const scalar minWidth) const
True if the obstacle is considered to be too small.
static scalar legacyReadFiles(const fileName &obsFileDir, const wordList &obsFileNames, const boundBox &meshBb, DynamicList< PDRobstacle > &blocks, DynamicList< PDRobstacle > &cylinders)
Read obstacle files and add to the lists.
void scale(const scalar factor)
Scale obstacle dimensions by specified scaling factor.
scalar def_grating_slat_w
Default slat thickness grating.
Definition: PDRparams.H:124
scalar scale
Overall scale factor.
Definition: PDRparams.H:140
scalar obs_expand
Definition: PDRparams.H:121
bool ignoreGratings
Definition: PDRparams.H:83
scalar min_width
Ignore obstacles with second dimension (or diameter) less than this.
Definition: PDRparams.H:113
const T & second() const noexcept
Return second element, which is also the last element.
Definition: PairI.H:120
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
Particle-size distribution model wherein random samples are drawn from the doubly-truncated uniform p...
Definition: uniform.H:164
A class for handling file names.
Definition: fileName.H:76
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition: volumeType.H:61
A class for handling words, derived from Foam::string.
Definition: word.H:68
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
labelPair readObstacleFiles(const fileName &obsFileDir, const wordList &obsFileNames, Map< obstacleGrouping > &groups)
Read obstacle files, do counting only.
constexpr const char *const group
Group name for atomic constants.
Different types of constants.
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
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::PDRparams pars
Globals for program parameters (ugly hack)
error FatalError
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
messageStream Warning
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
List< treeBoundBox > meshBb(1, treeBoundBox(boundBox(coarseMesh.points(), false)).extend(rndGen, 1e-3))
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278
A non-counting (dummy) refCount.
Definition: refCount.H:59