37 inline scalar
calcGexp(
const scalar expRatio,
const label nDiv)
39 return nDiv > 1 ?
pow(expRatio, 1.0/(nDiv - 1)) : 0.0;
46 bool Foam::PDRblock::checkMonotonic
49 const UList<scalar>& pts
52 const label len = pts.size();
59 const scalar& minVal = pts[0];
61 for (
label i=1; i < len; ++i)
67 <<
" direction do not increase monotonically" <<
nl
79 return NullObjectRef<PDRblock>();
85 void Foam::PDRblock::adjustSizes()
88 sizes().x() = grid_.x().nCells();
89 sizes().y() = grid_.y().nCells();
90 sizes().z() = grid_.z().nCells();
92 if (sizes().
x() <= 0 || sizes().
y() <= 0 || sizes().z() <= 0)
107 bounds_.min().x() = grid_.x().first();
108 bounds_.min().y() = grid_.y().first();
109 bounds_.min().z() = grid_.z().first();
111 bounds_.max().x() = grid_.x().last();
112 bounds_.max().y() = grid_.y().last();
113 bounds_.max().z() = grid_.z().last();
120 const label nEdge = grid_[cmpt].nCells();
122 for (
label edgei=0; edgei < nEdge; ++edgei)
124 const scalar len = grid_[cmpt].width(edgei);
125 minEdgeLen_ =
min(minEdgeLen_, len);
131 void Foam::PDRblock::readGridControl
134 const dictionary&
dict,
135 const scalar scaleFactor
149 Info<<
"Reading grid control for "
158 const label nSegments = (knots.size()-1);
163 <<
"Must define at least two control points:"
165 <<
" direction" <<
nl
172 for (scalar& pt : knots)
179 checkMonotonic(cmpt, knots);
192 label nTotalDivs = 0;
200 <<
"Negative or zero divisions:"
202 <<
" direction" <<
nl
207 if (divisions.size() != nSegments)
210 <<
"Mismatch in number of divisions and number of points:"
212 <<
" direction" <<
nl
215 else if (!nTotalDivs)
218 <<
"No grid divisions at all:"
220 <<
" direction" <<
nl
237 for (scalar& expRatio : expansion)
241 expRatio = 1.0/(-expRatio);
245 if (expansion.size() && expansion.size() != nSegments)
248 <<
"Mismatch in number of expansion ratios and divisions:"
250 <<
" direction" <<
nl
254 if (expansion.empty())
256 expansion.resize(nSegments, scalar(1));
260 Info<<
"Warning: 'ratios' not specified, using constant spacing"
277 gridPoint.
resize(nTotalDivs+1);
281 for (
label segmenti=0; segmenti < nSegments; ++segmenti)
283 const label nDiv = divisions[segmenti];
284 const scalar expRatio = expansion[segmenti];
286 SubList<scalar> subPoint(gridPoint, nDiv+1,
start);
289 subPoint[0] = knots[segmenti];
290 subPoint[nDiv] = knots[segmenti+1];
292 const scalar dist = (subPoint.last() - subPoint.first());
294 if (
equal(expRatio, 1))
296 for (
label i=1; i < nDiv; ++i)
298 subPoint[i] = (subPoint[0] + (dist * i)/nDiv);
304 const scalar expFact =
calcGexp(expRatio, nDiv);
306 for (
label i=1; i < nDiv; ++i)
311 + dist * (1.0 -
pow(expFact, i))/(1.0 -
pow(expFact, nDiv))
319 void Foam::PDRblock::readBoundary(
const dictionary&
dict)
323 Info<<
"Reading boundary entries" <<
nl;
326 PtrList<entry> patchEntries;
330 PtrList<entry> newEntries(
dict.
lookup(
"boundary"));
331 patchEntries.transfer(newEntries);
335 patches_.resize(patchEntries.size());
338 const labelRange validFaceId(0, 6);
343 forAll(patchEntries, patchi)
345 if (!patchEntries.set(patchi))
350 const entry&
e = patchEntries[patchi];
356 <<
" in boundary section is not a valid dictionary."
360 const dictionary& patchDict =
e.dict();
362 const word& patchName =
e.keyword();
364 const word patchType(patchDict.get<word>(
"type"));
370 if (faceLabels.size() != patchFaces.size())
373 <<
"Patch: " << patchName
374 <<
": Ignoring duplicate face ids" <<
nl;
377 label nErased = patchFaces.filterKeys(validFaceId);
381 <<
"Patch: " << patchName <<
": Ignoring " << nErased
382 <<
" faces with invalid identifiers" <<
nl;
385 nErased = patchFaces.erase(handled);
389 <<
"Patch: " << patchName <<
": Ignoring " << nErased
390 <<
" faces, which were already handled" <<
nl;
394 handled += patchFaces;
399 patches_.set(patchi,
new boundaryEntry());
401 boundaryEntry& bentry = patches_[patchi];
403 bentry.name_ = patchName;
404 bentry.type_ = patchType;
406 bentry.faces_ = patchFaces.sortedToc();
409 for (
const label shapeFacei : bentry.faces_)
411 bentry.size_ += nBoundaryFaces(shapeFacei);
419 missed.erase(handled);
423 patches_.append(
new boundaryEntry());
425 boundaryEntry& bentry = patches_.last();
427 bentry.name_ =
"defaultFaces";
428 bentry.type_ = emptyPolyPatch::typeName;
430 bentry.faces_ = missed.sortedToc();
433 for (
const label shapeFacei : bentry.faces_)
435 bentry.size_ += nBoundaryFaces(shapeFacei);
442 const dictionary* defaultEntry =
dict.
findDict(
"defaultPatch");
446 defaultEntry->readIfPresent(
"type", bentry.type_);
453 for (
const boundaryEntry& bentry : patches_)
455 Info<<
" patch: " << patchi
456 <<
" (size: " << bentry.size_
457 <<
" type: " << bentry.type_
458 <<
") name: " << bentry.name_
480 for (
label patchi=0; patchi < 6; ++patchi)
482 patches_.set(patchi,
new boundaryEntry());
484 boundaryEntry& bentry = patches_[patchi];
487 bentry.type_ =
"patch";
492 reset(xgrid, ygrid, zgrid);
513 readGridControl(0,
dict.
subDict(
"x"), scaleFactor);
514 readGridControl(1,
dict.
subDict(
"y"), scaleFactor);
515 readGridControl(2,
dict.
subDict(
"z"), scaleFactor);
539 checkMonotonic(cmpt, grid_[cmpt]);
546 for (boundaryEntry& bentry : patches_)
551 for (
const label shapeFacei : bentry.faces_)
553 bentry.size_ += nBoundaryFaces(shapeFacei);
568 if (!bounds_.contains(pt))
583 bool Foam::PDRblock::gridIndex
590 const scalar tol = relTol * minEdgeLen_;
595 pos[cmpt] = grid_[cmpt].findIndex(pt[cmpt], tol);
597 if (
pos[cmpt] < 0)
return false;
608 if (findCell(pt,
pos))
625 if (gridIndex(pt,
pos, relTol))