From b1bbd50fe4489a9a006096a4b22d0c52244e2058 Mon Sep 17 00:00:00 2001 From: Mark Olesen <Mark.Olesen@esi-group.com> Date: Mon, 25 Mar 2019 12:53:45 +0100 Subject: [PATCH] ENH: additional constructor and methods for PDRblock (#1216) - construct from positions - reset positions, minEdgeLen, find grid index ENH: add accessor for obtaining the size of a single i-j-k dimension - eg, obtain the 'i' dimension with any of these methods: ijkAddr.size(vector::X) ijkAddr.sizes().x() ijkAddr.sizes()[0] --- src/OpenFOAM/meshes/ijkMesh/IjkField.H | 6 + src/OpenFOAM/meshes/ijkMesh/IjkFieldI.H | 10 + src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H | 18 +- src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H | 13 +- src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H | 4 +- src/mesh/blockMesh/Make/files | 1 + src/mesh/blockMesh/PDRblockMesh/PDRblock.C | 245 +++++++++++++++--- src/mesh/blockMesh/PDRblockMesh/PDRblock.H | 82 +++++- src/mesh/blockMesh/PDRblockMesh/PDRblockI.H | 65 ++++- .../blockMesh/PDRblockMesh/PDRblockLocation.C | 108 ++++++++ 10 files changed, 487 insertions(+), 65 deletions(-) create mode 100644 src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C diff --git a/src/OpenFOAM/meshes/ijkMesh/IjkField.H b/src/OpenFOAM/meshes/ijkMesh/IjkField.H index 0e21e473e1b..b748b0ea5cc 100644 --- a/src/OpenFOAM/meshes/ijkMesh/IjkField.H +++ b/src/OpenFOAM/meshes/ijkMesh/IjkField.H @@ -104,6 +104,12 @@ public: //- Return i,j,k addressing sizes for modification inline labelVector& sizes(); + //- The field size + using Field<Type>::size; + + //- The addressing dimension in the given direction + inline const label& size(const vector::components cmpt) const; + // Edit diff --git a/src/OpenFOAM/meshes/ijkMesh/IjkFieldI.H b/src/OpenFOAM/meshes/ijkMesh/IjkFieldI.H index 8943bfa5dcf..09da6a5c323 100644 --- a/src/OpenFOAM/meshes/ijkMesh/IjkFieldI.H +++ b/src/OpenFOAM/meshes/ijkMesh/IjkFieldI.H @@ -157,6 +157,16 @@ inline Foam::labelVector& Foam::IjkField<Type>::sizes() } +template<class Type> +inline const Foam::label& Foam::IjkField<Type>::size +( + const vector::components cmpt +) const +{ + return ijk_.size(cmpt); +} + + template<class Type> inline void Foam::IjkField<Type>::clear() { diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H b/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H index 215acbd8ad8..93889dce082 100644 --- a/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H +++ b/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H @@ -36,6 +36,7 @@ SourceFiles #define ijkAddressing_H #include "labelVector.H" +#include "vector.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -72,17 +73,20 @@ public: // Access - //- Return the i,j,k addressing sizes + //- Addressing is considered empty if any component is zero + inline bool empty() const; + + //- The (i,j,k) addressing dimensions inline const labelVector& sizes() const; - //- Return the i,j,k addressing sizes for modification + //- Return the (i,j,k) dimensions for modification inline labelVector& sizes(); //- Return the total i*j*k size inline label size() const; - //- Addressing is considered empty if any component is zero - inline bool empty() const; + //- The addressing dimension in the given direction + inline const label& size(const vector::components cmpt) const; //- Reset to (0,0,0) sizing inline void clear(); @@ -93,13 +97,13 @@ public: //- Change the sizing parameters inline void reset(const labelVector& newSizes); - //- Linear addressing index (offset) for an i,j,k position. + //- Linear addressing index (offset) for an (i,j,k) position. inline label index(const label i, const label j, const label k) const; - //- Linear addressing index (offset) for an i,j,k position. + //- Linear addressing index (offset) for an (i,j,k) position. inline label index(const labelVector& ijk) const; - //- The i,j,k indexing from linear addressing. + //- The (i,j,k) indexing from linear addressing. inline labelVector index(const label idx) const; diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H b/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H index adb11da7935..9affedece16 100644 --- a/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H +++ b/src/OpenFOAM/meshes/ijkMesh/ijkAddressingI.H @@ -58,6 +58,12 @@ inline Foam::ijkAddressing::ijkAddressing // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +inline bool Foam::ijkAddressing::empty() const +{ + return (!sizes_.x() || !sizes_.y() || !sizes_.z()); +} + + inline const Foam::labelVector& Foam::ijkAddressing::sizes() const { return sizes_; @@ -77,9 +83,12 @@ inline Foam::label Foam::ijkAddressing::size() const } -inline bool Foam::ijkAddressing::empty() const +inline const Foam::label& Foam::ijkAddressing::size +( + const vector::components cmpt +) const { - return (!sizes_.x() || !sizes_.y() || !sizes_.z()); + return sizes_[cmpt]; } diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H b/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H index da73e24be32..7a5c2af7314 100644 --- a/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H +++ b/src/OpenFOAM/meshes/ijkMesh/ijkMeshI.H @@ -143,13 +143,13 @@ inline Foam::label Foam::ijkMesh::nBoundaryFaces return n.y()*n.z(); break; - // Face 2, 3 == y-min, y-max + // Face 2,3 == y-min, y-max case 2: case 3: return n.z()*n.x(); break; - // Face 4, 5 == z-min, z-max + // Face 4,5 == z-min, z-max case 4: case 5: return n.x()*n.y(); diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files index 84473bfc0f2..8193f56e1a1 100644 --- a/src/mesh/blockMesh/Make/files +++ b/src/mesh/blockMesh/Make/files @@ -41,5 +41,6 @@ blockMeshTools/blockMeshTools.C PDRblockMesh/PDRblock.C PDRblockMesh/PDRblockCreate.C +PDRblockMesh/PDRblockLocation.C LIB = $(FOAM_LIBBIN)/libblockMesh diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C index cca8eb39a29..b591fa42e94 100644 --- a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C @@ -39,8 +39,87 @@ namespace Foam } +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +bool Foam::PDRblock::checkMonotonic +( + const direction cmpt, + const UList<scalar>& 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; +} + + // * * * * * * * * * * * * * 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, @@ -89,22 +168,7 @@ void Foam::PDRblock::readGridControl } // Do points increase monotonically? - { - const scalar& minVal = knots.first(); - - for (label pointi = 1; pointi < knots.size(); ++pointi) - { - if (knots[pointi] <= minVal) - { - FatalIOErrorInFunction(dict) - << "Points are not monotonically increasing:" - << " in " << vector::componentNames[cmpt] - << " direction" << nl - << exit(FatalIOError); - } - } - } - + checkMonotonic(cmpt, knots); if (verbose_) { @@ -394,6 +458,33 @@ void Foam::PDRblock::readBoundary(const dictionary& dict) // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +Foam::PDRblock::PDRblock +( + const UList<scalar>& xgrid, + const UList<scalar>& ygrid, + const UList<scalar>& 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() @@ -411,50 +502,124 @@ bool Foam::PDRblock::read(const dictionary& dict) // Mark no scaling with invalid value const scalar scaleFactor = dict.lookupOrDefault<scalar>("scale", -1); - // Grid controls readGridControl(0, dict.subDict("x"), scaleFactor); readGridControl(1, dict.subDict("y"), scaleFactor); readGridControl(2, dict.subDict("z"), scaleFactor); - // Adjust i-j-k addressing - sizes().x() = grid_.x().nCells(); - sizes().y() = grid_.y().nCells(); - sizes().z() = grid_.z().nCells(); + adjustSizes(); - // Adjust boundBox - bounds_.min().x() = grid_.x().first(); - bounds_.min().y() = grid_.y().first(); - bounds_.min().z() = grid_.z().first(); + readBoundary(dict); - bounds_.max().x() = grid_.x().last(); - bounds_.max().y() = grid_.y().last(); - bounds_.max().z() = grid_.z().last(); + return true; +} - // Boundaries - readBoundary(dict); +void Foam::PDRblock::reset +( + const UList<scalar>& xgrid, + const UList<scalar>& ygrid, + const UList<scalar>& zgrid +) +{ + static_cast<scalarList&>(grid_.x()) = xgrid; + static_cast<scalarList&>(grid_.y()) = ygrid; + static_cast<scalarList&>(grid_.z()) = zgrid; - return true; + #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); + } + } } -Foam::labelVector Foam::PDRblock::findCell(const point& pt) const +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. - labelVector pos(-1,-1,-1); - if (bounds_.contains(pt)) + // Since findLower returns the lower index, it corresponds to the + // cell in which the point is found + + if (!bounds_.contains(pt)) { - for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt) - { - // Binary search - pos[cmpt] = findLower(grid_[cmpt], pt[cmpt]); - } + 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 pos; + return labelVector(-1,-1,-1); } diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H index 8ee0df830af..91e16b66e28 100644 --- a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H @@ -44,10 +44,10 @@ Description Grid coordinate controls \table - Property| Description | Required | Default - points | Locations defining the mesh segment | yes | - nCells | Divisions per mesh segment | yes | - ratios | Expansion ratio (end/start) per segment | no | uniform + Property| Description | Required | Default + points | Locations defining the mesh segment | yes | + nCells | Divisions per mesh segment | yes | + ratios | Expansion ratio (end/start) per segment | no | uniform \endtable The expansion ratios are defined as per blockMesh and represent the ratio @@ -98,14 +98,20 @@ public: : public scalarList { + //- The locations are valid if they contain 2 or more points + inline bool valid() const; + //- The number of cells in this direction. inline label nCells() const; //- The number of points in this direction. inline label nPoints() const; - //- The number of divisions (nCells) in this direction. - inline label size() const; + //- True if the location is within the range + inline bool contains(const scalar p) const; + + //- Mid-point location, zero for an empty list. + inline scalar centre() const; //- Check that element index is within valid range. inline void checkIndex(const label i) const; @@ -114,7 +120,22 @@ public: inline scalar width(const label i) const; //- Cell centre at element position. + // Treats -1 and nCells positions like a halo cell. inline scalar C(const label i) const; + + //- Find the cell index enclosing this location + // \return -1 for out-of-bounds + label findCell(const scalar p) const; + + //- Find the grid index, within the given tolerance + // Return -1 for out-of-bounds and -2 for a point that is + // within bounds, but not aligned with a grid point. + label findIndex(const scalar p, const scalar tol) const; + + //- If out of range, return the respective min/max limits, + //- otherwise return the value itself. + // If the range is invalid, always return the value. + inline const scalar& clip(const scalar& val) const; }; @@ -153,9 +174,18 @@ private: //- The boundary patch information PtrList<boundaryEntry> patches_; + //- The min edge length + scalar minEdgeLen_; + // Private Member Functions + //- Check that points increase monotonically + static bool checkMonotonic(const direction cmpt, const UList<scalar>& pts); + + //- Adjust sizing for updated grid points + void adjustSizes(); + //- Read and define grid points in given direction void readGridControl ( @@ -167,7 +197,6 @@ private: //- Read "boundary" information void readBoundary(const dictionary& dict); - //- Populate point field for the block void createPoints(pointField& pts) const; @@ -181,7 +210,6 @@ private: labelList::iterator& neiIter ) const; - //- Add boundary faces for the shape face to lists // Lists must be properly sized! // \return the number of faces added @@ -192,6 +220,19 @@ private: labelList::iterator& ownIter ) const; + //- Obtain i,j,k index for cell enclosing this location + // \return false for out-of-bounds + bool findCell(const point& pt, labelVector& pos) const; + + //- Obtain i,j,k grid index for point location + // \return false for out-of-bounds and off-grid + bool gridIndex + ( + const point& pt, + labelVector& pos, + const scalar tol + ) const; + public: @@ -200,6 +241,14 @@ public: //- Construct zero-size inline PDRblock(); + //- Construct from components + PDRblock + ( + const UList<scalar>& xgrid, + const UList<scalar>& ygrid, + const UList<scalar>& zgrid + ); + //- Construct from dictionary explicit PDRblock(const dictionary& dict, bool verbose=false); @@ -209,6 +258,14 @@ public: //- Read dictionary bool read(const dictionary& dict); + //- Reset grid locations and mesh i-j-k sizing + void reset + ( + const UList<scalar>& xgrid, + const UList<scalar>& ygrid, + const UList<scalar>& zgrid + ); + // Access @@ -221,6 +278,9 @@ public: //- The mesh bounding box inline const boundBox& bounds() const; + //- The min edge length + inline const scalar& minEdgeLen() const; + //- Cell size in x-direction at i position. inline scalar dx(const label i) const; @@ -272,6 +332,12 @@ public: // The value (-1,-1,-1) is returned for out-of-bounds (not found). labelVector findCell(const point& pt) const; + //- Obtain i,j,k grid index for point location within specified + // relative tolerance of the minEdgeLen. + // The value (-1,-1,-1) is returned for out-of-bounds (not found). + // and off-grid + labelVector gridIndex(const point& pt, const scalar relTol=0.01) const; + // Mesh generation diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H index 64becbd4b86..a624ee68a65 100644 --- a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H @@ -31,12 +31,19 @@ inline Foam::PDRblock::PDRblock() verbose_(false), grid_(), bounds_(), - patches_() + patches_(), + minEdgeLen_(Zero) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +inline bool Foam::PDRblock::location::valid() const +{ + return (scalarList::size() > 1); +} + + inline Foam::label Foam::PDRblock::location::nCells() const { return (scalarList::size()-1); @@ -49,19 +56,25 @@ inline Foam::label Foam::PDRblock::location::nPoints() const } -inline Foam::label Foam::PDRblock::location::size() const +inline bool Foam::PDRblock::location::contains(const scalar p) const { - return (scalarList::size()-1); + return (scalarList::size() > 1 && first() <= p && p <= last()); +} + + +inline Foam::scalar Foam::PDRblock::location::centre() const +{ + return scalarList::empty() ? 0 : (0.5*first() + 0.5*last()); } inline void Foam::PDRblock::location::checkIndex(const label i) const { - if (i < 0 || i >= size()) + if (i < 0 || i >= nCells()) { FatalErrorInFunction << "The index " << i - << " is out of range [0," << size() << ']' << nl + << " is out of range [0," << nCells() << ']' << nl << abort(FatalError); } } @@ -79,11 +92,45 @@ inline Foam::scalar Foam::PDRblock::location::width(const label i) const inline Foam::scalar Foam::PDRblock::location::C(const label i) const { + if (i == -1) + { + #ifdef FULLDEBUG + checkIndex(0); + #endif + + // "Halo" centre [-1] == x0 - 1/2 (x1 - x0) + return first() - 0.5*(operator[](1) - first()); + } + else if (i > 1 && i == scalarList::size()-1) + { + // "Halo" centre [nCells] == xN + 1/2 (xN - xN1) + return last() + 0.5*(last() - operator[](scalarList::size()-2)); + } + #ifdef FULLDEBUG checkIndex(i); #endif - return (operator[](i+1) + operator[](i)); + return 0.5*(operator[](i+1) + operator[](i)); +} + + +inline const Foam::scalar& +Foam::PDRblock::location::clip(const scalar& val) const +{ + if (scalarList::size()) + { + if (val < first()) + { + return first(); + } + else if (last() < val) + { + return last(); + } + } + + return val; // Pass-through } @@ -94,6 +141,12 @@ Foam::PDRblock::grid() const } +inline const Foam::scalar& Foam::PDRblock::minEdgeLen() const +{ + return minEdgeLen_; +} + + inline const Foam::boundBox& Foam::PDRblock::bounds() const { return bounds_; diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C new file mode 100644 index 00000000000..ca92808d837 --- /dev/null +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C @@ -0,0 +1,108 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "PDRblock.H" +#include "ListOps.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::label Foam::PDRblock::location::findCell(const scalar p) const +{ + if (scalarList::empty()) + { + return -1; + } + else if (equal(p, first())) + { + return 0; + } + else if (equal(p, last())) + { + return nCells()-1; + } + else if (p < first() || p > last()) + { + // The point is out-of-bounds + return -1; + } + + // Binary search, finds lower index and thus corresponds to the + // cell in which the point is found + return findLower(*this, p); +} + + +Foam::label Foam::PDRblock::location::findIndex +( + const scalar p, + const scalar tol +) const +{ + if (scalarList::empty()) + { + return -1; + } + else if (Foam::mag(p - first()) <= tol) + { + return 0; + } + else if (Foam::mag(p - last()) <= tol) + { + return scalarList::size()-1; + } + else if (p < first() || p > last()) + { + // The point is out-of-bounds + return -1; + } + + // Linear search + label i = 0; + scalar delta = GREAT; + + for (const scalar& val : *this) + { + const scalar diff = mag(p - val); + + if (diff <= tol) + { + return i; + } + else if (delta < diff) + { + // Moving further away + break; + } + + delta = diff; + ++i; + } + + // This point is within bounds, but not near a grid-point + return -2; +} + + +// ************************************************************************* // -- GitLab