34 #define Log if (verbose_) Info
43 Foam::PDRblock::expansionNames_
45 { expansionType::EXPAND_UNIFORM,
"uniform" },
46 { expansionType::EXPAND_RATIO,
"ratio" },
47 { expansionType::EXPAND_RELATIVE,
"relative"}
56 inline scalar
calcGexp(
const scalar expRatio,
const label nDiv)
58 return nDiv > 1 ?
pow(expRatio, 1.0/(nDiv - 1)) : 0.0;
64 const scalar expRatio,
68 return nDiv > 1 ?
pow(expRatio, (nDiv - 1)) : 1.0;
76 bool Foam::PDRblock::checkMonotonic
79 const UList<scalar>& pts
82 const label len = pts.size();
89 const scalar& minVal = pts[0];
91 for (label i=1; i < len; ++i)
97 <<
" direction do not increase monotonically" <<
nl
109 return NullObjectRef<PDRblock>();
115 void Foam::PDRblock::adjustSizes()
118 sizes().x() = grid_.x().nCells();
119 sizes().y() = grid_.y().nCells();
120 sizes().z() = grid_.z().nCells();
122 if (sizes().
x() <= 0 || sizes().
y() <= 0 || sizes().z() <= 0)
137 bounds_.min().x() = grid_.x().first();
138 bounds_.min().y() = grid_.y().first();
139 bounds_.min().z() = grid_.z().first();
141 bounds_.max().x() = grid_.x().last();
142 bounds_.max().y() = grid_.y().last();
143 bounds_.max().z() = grid_.z().last();
150 const label nEdge = grid_[cmpt].nCells();
152 for (label edgei=0; edgei < nEdge; ++edgei)
154 const scalar len = grid_[cmpt].width(edgei);
155 minEdgeLen_ =
min(minEdgeLen_, len);
161 void Foam::PDRblock::readGridControl
164 const dictionary&
dict,
165 const scalar scaleFactor,
166 expansionType expandType
178 Log <<
"Reading grid control for "
186 const label nSegments = (knots.size()-1);
191 <<
"Must define at least two control points:"
193 <<
" direction" <<
nl
200 for (scalar& pt : knots)
207 checkMonotonic(cmpt, knots);
217 label nTotalDivs = 0;
218 for (
const label
ndiv : divisions)
225 <<
"Negative or zero divisions:"
227 <<
" direction" <<
nl
232 if (divisions.size() != nSegments)
235 <<
"Mismatch in number of divisions and number of points:"
237 <<
" direction" <<
nl
240 else if (!nTotalDivs)
243 <<
"No grid divisions at all:"
245 <<
" direction" <<
nl
258 for (scalar& expRatio : expansion)
262 expRatio = 1.0/(-expRatio);
266 if (expansion.size() && divisions.size() != expansion.size())
269 <<
"Mismatch in number of expansion ratios and divisions:"
271 <<
" direction" <<
nl
275 if (expansion.empty())
277 expansion.resize(nSegments, scalar(1));
279 if (expandType != expansionType::EXPAND_UNIFORM)
281 expandType = expansionType::EXPAND_UNIFORM;
283 Log <<
"Warning: no 'ratios', use uniform spacing" <<
nl;
290 case expansionType::EXPAND_UNIFORM:
292 expansion = scalar(1);
296 case expansionType::EXPAND_RATIO:
302 case expansionType::EXPAND_RELATIVE:
306 auto divIter = divisions.cbegin();
308 for (scalar& expRatio : expansion)
324 gridPoint.
resize(nTotalDivs+1);
328 for (label segmenti=0; segmenti < nSegments; ++segmenti)
330 const label nDiv = divisions[segmenti];
331 const scalar expRatio = expansion[segmenti];
333 SubList<scalar> subPoint(gridPoint, nDiv+1, start);
336 subPoint[0] = knots[segmenti];
337 subPoint[nDiv] = knots[segmenti+1];
339 const scalar dist = (subPoint.last() - subPoint.first());
343 expandType == expansionType::EXPAND_UNIFORM
344 ||
equal(expRatio, 1)
347 for (label i=1; i < nDiv; ++i)
349 subPoint[i] = (subPoint[0] + (dist * i)/nDiv);
355 const scalar expFact =
calcGexp(expRatio, nDiv);
357 for (label i=1; i < nDiv; ++i)
362 + dist * (1.0 -
pow(expFact, i))/(1.0 -
pow(expFact, nDiv))
370 void Foam::PDRblock::readBoundary(
const dictionary&
dict)
372 Log <<
"Reading boundary entries" <<
nl;
374 PtrList<entry> patchEntries;
378 PtrList<entry> newEntries(
dict.
lookup(
"boundary"));
379 patchEntries.transfer(newEntries);
383 patches_.resize(patchEntries.size());
386 const labelRange validFaceId(0, 6);
391 forAll(patchEntries, patchi)
393 if (!patchEntries.set(patchi))
398 const entry&
e = patchEntries[patchi];
404 <<
" in boundary section is not a valid dictionary."
408 const dictionary& patchDict =
e.dict();
410 const word& patchName =
e.keyword();
412 const word patchType(patchDict.get<word>(
"type"));
418 if (faceLabels.size() != patchFaces.size())
421 <<
"Patch: " << patchName
422 <<
": Ignoring duplicate face ids" <<
nl;
425 label nErased = patchFaces.filterKeys(validFaceId);
429 <<
"Patch: " << patchName <<
": Ignoring " << nErased
430 <<
" faces with invalid identifiers" <<
nl;
433 nErased = patchFaces.erase(handled);
437 <<
"Patch: " << patchName <<
": Ignoring " << nErased
438 <<
" faces, which were already handled" <<
nl;
442 handled += patchFaces;
447 patches_.set(patchi,
new boundaryEntry());
449 boundaryEntry& bentry = patches_[patchi];
451 bentry.name_ = patchName;
452 bentry.type_ = patchType;
454 bentry.faces_ = patchFaces.sortedToc();
457 for (
const label shapeFacei : bentry.faces_)
459 bentry.size_ += nBoundaryFaces(shapeFacei);
467 missed.erase(handled);
471 patches_.append(
new boundaryEntry());
473 boundaryEntry& bentry = patches_.last();
475 bentry.name_ =
"defaultFaces";
476 bentry.type_ = emptyPolyPatch::typeName;
478 bentry.faces_ = missed.sortedToc();
481 for (
const label shapeFacei : bentry.faces_)
483 bentry.size_ += nBoundaryFaces(shapeFacei);
490 const dictionary* defaultEntry =
dict.
findDict(
"defaultPatch");
494 defaultEntry->readIfPresent(
"type", bentry.type_);
501 for (
const boundaryEntry& bentry : patches_)
503 Info<<
" patch: " << patchi
504 <<
" (size: " << bentry.size_
505 <<
" type: " << bentry.type_
506 <<
") name: " << bentry.name_
528 for (label patchi=0; patchi < 6; ++patchi)
530 patches_.set(patchi,
new boundaryEntry());
532 boundaryEntry& bentry = patches_[patchi];
535 bentry.type_ =
"patch";
540 reset(xgrid, ygrid, zgrid);
562 expansionNames_.getOrDefault
566 expansionType::EXPAND_RATIO
570 readGridControl(0,
dict.
subDict(
"x"), scaleFactor, expandType);
571 readGridControl(1,
dict.
subDict(
"y"), scaleFactor, expandType);
572 readGridControl(2,
dict.
subDict(
"z"), scaleFactor, expandType);
596 checkMonotonic(cmpt, grid_[cmpt]);
603 for (boundaryEntry& bentry : patches_)
608 for (
const label shapeFacei : bentry.faces_)
610 bentry.size_ += nBoundaryFaces(shapeFacei);
625 if (!bounds_.contains(pt))
640 bool Foam::PDRblock::gridIndex
647 const scalar tol = relTol * minEdgeLen_;
652 pos[cmpt] = grid_[cmpt].findIndex(pt[cmpt], tol);
654 if (
pos[cmpt] < 0)
return false;
665 if (findCell(pt,
pos))
682 if (gridIndex(pt,
pos, relTol))