/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- 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 . \*---------------------------------------------------------------------------*/ #include "PDRblock.H" #include "ListOps.H" #include "emptyPolyPatch.H" // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // namespace Foam { //- Calculate geometric expansion factor from expansion ratio inline scalar calcGexp(const scalar expRatio, const label nDiv) { return nDiv > 1 ? pow(expRatio, 1.0/(nDiv - 1)) : 0.0; } } // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // bool Foam::PDRblock::checkMonotonic ( const direction cmpt, const UList& pts ) { const label len = pts.size(); if (!len) { return false; } const scalar& minVal = pts[0]; for (label i=1; i < len; ++i) { if (pts[i] <= minVal) { FatalErrorInFunction << "Points in " << vector::componentNames[cmpt] << " direction do not increase monotonically" << nl << flatOutput(pts) << nl << nl << exit(FatalError); } } return true; } const Foam::PDRblock& Foam::PDRblock::null() { return NullObjectRef(); } // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // void Foam::PDRblock::adjustSizes() { // Adjust i-j-k addressing sizes().x() = grid_.x().nCells(); sizes().y() = grid_.y().nCells(); sizes().z() = grid_.z().nCells(); if (sizes().x() <= 0 || sizes().y() <= 0 || sizes().z() <= 0) { // Sanity check. Silently disallow bad sizing ijkMesh::clear(); grid_.x().clear(); grid_.y().clear(); grid_.z().clear(); bounds_ = boundBox::invertedBox; minEdgeLen_ = Zero; return; } // Adjust boundBox bounds_.min().x() = grid_.x().first(); bounds_.min().y() = grid_.y().first(); bounds_.min().z() = grid_.z().first(); bounds_.max().x() = grid_.x().last(); bounds_.max().y() = grid_.y().last(); bounds_.max().z() = grid_.z().last(); // Min edge length minEdgeLen_ = GREAT; for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt) { const label nEdge = grid_[cmpt].nCells(); for (label edgei=0; edgei < nEdge; ++edgei) { const scalar len = grid_[cmpt].width(edgei); minEdgeLen_ = min(minEdgeLen_, len); } } } void Foam::PDRblock::readGridControl ( const direction cmpt, const dictionary& dict, const scalar scaleFactor ) { //- The begin/end nodes for each segment scalarList knots; // The number of division per segment labelList divisions; // The expansion ratio per segment scalarList expansion; if (verbose_) { Info<< "Reading grid control for " << vector::componentNames[cmpt] << " direction" << nl; } // Points // ~~~~~~ dict.readEntry("points", knots); const label nSegments = (knots.size()-1); if (nSegments < 1) { FatalIOErrorInFunction(dict) << "Must define at least two control points:" << " in " << vector::componentNames[cmpt] << " direction" << nl << exit(FatalIOError); } // Apply scaling if (scaleFactor > 0) { for (scalar& pt : knots) { pt *= scaleFactor; } } // Do points increase monotonically? checkMonotonic(cmpt, knots); if (verbose_) { Info<< " points : " << flatOutput(knots) << nl; } // Divisions // ~~~~~~~~~ dict.readEntry("nCells", divisions); label nTotalDivs = 0; for (const label ndiv : divisions) { nTotalDivs += ndiv; if (ndiv < 1) { FatalIOErrorInFunction(dict) << "Negative or zero divisions:" << " in " << vector::componentNames[cmpt] << " direction" << nl << exit(FatalIOError); } } if (divisions.size() != nSegments) { FatalIOErrorInFunction(dict) << "Mismatch in number of divisions and number of points:" << " in " << vector::componentNames[cmpt] << " direction" << nl << exit(FatalIOError); } else if (!nTotalDivs) { FatalIOErrorInFunction(dict) << "No grid divisions at all:" << " in " << vector::componentNames[cmpt] << " direction" << nl << exit(FatalIOError); } if (verbose_) { Info<< " nCells : " << flatOutput(divisions) << nl; } // Expansion ratios // ~~~~~~~~~~~~~~~~ dict.readIfPresent("ratios", expansion); // Correct expansion ratios - negative is the same as inverse. for (scalar& expRatio : expansion) { if (expRatio < 0) { expRatio = 1.0/(-expRatio); } } if (expansion.size() && expansion.size() != nSegments) { FatalIOErrorInFunction(dict) << "Mismatch in number of expansion ratios and divisions:" << " in " << vector::componentNames[cmpt] << " direction" << nl << exit(FatalIOError); } if (expansion.empty()) { expansion.resize(nSegments, scalar(1)); if (verbose_) { Info<< "Warning: 'ratios' not specified, using constant spacing" << nl; } } else { if (verbose_) { Info<< " ratios : " << flatOutput(expansion) << nl; } } // Define the grid points scalarList& gridPoint = grid_[cmpt]; gridPoint.resize(nTotalDivs+1); label start = 0; for (label segmenti=0; segmenti < nSegments; ++segmenti) { const label nDiv = divisions[segmenti]; const scalar expRatio = expansion[segmenti]; SubList subPoint(gridPoint, nDiv+1, start); start += nDiv; subPoint[0] = knots[segmenti]; subPoint[nDiv] = knots[segmenti+1]; const scalar dist = (subPoint.last() - subPoint.first()); if (equal(expRatio, 1.0)) { for (label i=1; i < nDiv; ++i) { subPoint[i] = (subPoint[0] + (dist * i)/nDiv); } } else { // Geometric expansion factor from the expansion ratio const scalar expFact = calcGexp(expRatio, nDiv); for (label i=1; i < nDiv; ++i) { subPoint[i] = ( subPoint[0] + dist * (1.0 - pow(expFact, i))/(1.0 - pow(expFact, nDiv)) ); } } } } void Foam::PDRblock::readBoundary(const dictionary& dict) { if (verbose_) { Info<< "Reading boundary entries" << nl; } PtrList patchEntries; if (dict.found("boundary")) { PtrList newEntries(dict.lookup("boundary")); patchEntries.transfer(newEntries); } patches_.clear(); patches_.resize(patchEntries.size()); // Hex cell has 6 sides const labelRange validFaceId(0, 6); // Track which sides have been handled labelHashSet handled; forAll(patchEntries, patchi) { if (!patchEntries.set(patchi)) { continue; } const entry& e = patchEntries[patchi]; if (!e.isDict()) { FatalIOErrorInFunction(dict) << "Entry " << e << " in boundary section is not a valid dictionary." << exit(FatalIOError); } const dictionary& patchDict = e.dict(); const word& patchName = e.keyword(); const word patchType(patchDict.get("type")); labelList faceLabels(patchDict.get("faces")); labelHashSet patchFaces(faceLabels); if (faceLabels.size() != patchFaces.size()) { WarningInFunction << "Patch: " << patchName << ": Ignoring duplicate face ids" << nl; } label nErased = patchFaces.filterKeys(validFaceId); if (nErased) { WarningInFunction << "Patch: " << patchName << ": Ignoring " << nErased << " faces with invalid identifiers" << nl; } nErased = patchFaces.erase(handled); if (nErased) { WarningInFunction << "Patch: " << patchName << ": Ignoring " << nErased << " faces, which were already handled" << nl; } // Mark these faces as having been handled handled += patchFaces; // Save information for later access during mesh creation. patches_.set(patchi, new boundaryEntry()); boundaryEntry& bentry = patches_[patchi]; bentry.name_ = patchName; bentry.type_ = patchType; bentry.size_ = 0; bentry.faces_ = patchFaces.sortedToc(); // Count patch faces for (const label shapeFacei : bentry.faces_) { bentry.size_ += nBoundaryFaces(shapeFacei); } } // Deal with unhandled faces labelHashSet missed(identity(6)); missed.erase(handled); if (missed.size()) { patches_.append(new boundaryEntry()); boundaryEntry& bentry = patches_.last(); bentry.name_ = "defaultFaces"; bentry.type_ = emptyPolyPatch::typeName; bentry.size_ = 0; bentry.faces_ = missed.sortedToc(); // Count patch faces for (const label shapeFacei : bentry.faces_) { bentry.size_ += nBoundaryFaces(shapeFacei); } // Use name/type from "defaultPatch" entry if available // - be generous with error handling const dictionary* defaultEntry = dict.findDict("defaultPatch"); if (defaultEntry) { defaultEntry->readIfPresent("name", bentry.name_); defaultEntry->readIfPresent("type", bentry.type_); } } if (verbose_) { label patchi = 0; for (const boundaryEntry& bentry : patches_) { Info<< " patch: " << patchi << " (size: " << bentry.size_ << " type: " << bentry.type_ << ") name: " << bentry.name_ << " faces: " << flatOutput(bentry.faces_) << nl; ++patchi; } } } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::PDRblock::PDRblock ( const UList& xgrid, const UList& ygrid, const UList& zgrid ) : PDRblock() { // Default boundaries with patchi == shapeFacei patches_.resize(6); for (label patchi=0; patchi < 6; ++patchi) { patches_.set(patchi, new boundaryEntry()); boundaryEntry& bentry = patches_[patchi]; bentry.name_ = "patch" + Foam::name(patchi); bentry.type_ = "patch"; bentry.size_ = 0; bentry.faces_ = labelList(one(), patchi); } reset(xgrid, ygrid, zgrid); } Foam::PDRblock::PDRblock(const dictionary& dict, bool verbose) : PDRblock() { verbose_ = verbose; read(dict); } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // bool Foam::PDRblock::read(const dictionary& dict) { // Mark no scaling with invalid value const scalar scaleFactor = dict.lookupOrDefault("scale", -1); readGridControl(0, dict.subDict("x"), scaleFactor); readGridControl(1, dict.subDict("y"), scaleFactor); readGridControl(2, dict.subDict("z"), scaleFactor); adjustSizes(); readBoundary(dict); return true; } void Foam::PDRblock::reset ( const UList& xgrid, const UList& ygrid, const UList& zgrid ) { static_cast(grid_.x()) = xgrid; static_cast(grid_.y()) = ygrid; static_cast(grid_.z()) = zgrid; #ifdef FULLDEBUG for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt) { checkMonotonic(cmpt, grid_[cmpt]); } #endif adjustSizes(); // Adjust boundaries for (boundaryEntry& bentry : patches_) { bentry.size_ = 0; // Count patch faces for (const label shapeFacei : bentry.faces_) { bentry.size_ += nBoundaryFaces(shapeFacei); } } } bool Foam::PDRblock::findCell(const point& pt, labelVector& pos) const { // Out-of-bounds is handled explicitly, for efficiency and consistency, // but principally to ensure that findLower() returns a valid // result when the point is to the right of the bounds. // Since findLower returns the lower index, it corresponds to the // cell in which the point is found if (!bounds_.contains(pt)) { return false; } for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt) { // Binary search pos[cmpt] = findLower(grid_[cmpt], pt[cmpt]); } return true; } bool Foam::PDRblock::gridIndex ( const point& pt, labelVector& pos, const scalar relTol ) const { const scalar tol = relTol * minEdgeLen_; for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt) { // Linear search pos[cmpt] = grid_[cmpt].findIndex(pt[cmpt], tol); if (pos[cmpt] < 0) return false; } return true; } Foam::labelVector Foam::PDRblock::findCell(const point& pt) const { labelVector pos; if (findCell(pt, pos)) { return pos; } return labelVector(-1,-1,-1); } Foam::labelVector Foam::PDRblock::gridIndex ( const point& pt, const scalar relTol ) const { labelVector pos; if (gridIndex(pt, pos, relTol)) { return pos; } return labelVector(-1,-1,-1); } // ************************************************************************* //