From f8d08a805ba9ac59e4a64e84ea0657dfc8956a6c Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Wed, 21 Oct 2020 22:10:41 +0200 Subject: [PATCH] ENH: additional handling for out-of-range sampling (#1891) - when sampling onto a meshed surface, the sampling surface may be outside of the mesh region, or simply too far away to be considered reasonable. Can now specify a max search distance and default values for samples that are too distant. If a default value is not specified, uses Type(Zero). Eg, maxDistance 0.005; defaultValue { "p.*" 1e5; T 273.15; U (-100 -100 -100); } --- .../sampledMeshedSurface.C | 144 +++++++++++------- .../sampledMeshedSurface.H | 67 +++++--- .../sampledMeshedSurfaceTemplates.C | 93 +++++++---- .../sampledSurface/sampledSurface.H | 8 +- .../sampledSurface/sampledSurfaceTemplates.C | 16 +- .../rhoSimpleFoam/squareBend/system/sampling | 38 +++++ 6 files changed, 250 insertions(+), 116 deletions(-) diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C index 9830373ee5..2eed7d3d43 100644 --- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C +++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.C @@ -170,7 +170,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) // Does approximation by looking at the face centres only const pointField& fc = surface_.faceCentres(); - List nearest(fc.size(), nearInfo(GREAT, labelMax)); + List nearest(fc.size(), nearInfo(Foam::sqr(GREAT), labelMax)); if (sampleSource_ == samplingSource::cells) { @@ -181,13 +181,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) forAll(fc, facei) { const point& pt = fc[facei]; + auto& near = nearest[facei]; - pointIndexHit info = cellTree.findNearest(pt, sqr(GREAT)); + pointIndexHit info = cellTree.findNearest(pt, Foam::sqr(GREAT)); if (info.hit()) { - nearest[facei].first() = magSqr(info.hitPoint()-pt); - nearest[facei].second() = globalCells.toGlobal(info.index()); + near.first() = magSqr(info.hitPoint()-pt); + near.second() = globalCells.toGlobal(info.index()); } } } @@ -200,14 +201,15 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) forAll(fc, facei) { const point& pt = fc[facei]; + auto& near = nearest[facei]; if (cellTree.bb().contains(pt)) { const label index = cellTree.findInside(pt); if (index != -1) { - nearest[facei].first() = 0; - nearest[facei].second() = globalCells.toGlobal(index); + near.first() = 0; + near.second() = globalCells.toGlobal(index); } } } @@ -222,13 +224,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) forAll(fc, facei) { const point& pt = fc[facei]; + auto& near = nearest[facei]; pointIndexHit info = bndTree.findNearest(pt, sqr(GREAT)); if (info.hit()) { - nearest[facei].first() = magSqr(info.hitPoint()-pt); - nearest[facei].second() = + near.first() = magSqr(info.hitPoint()-pt); + near.second() = globalCells.toGlobal ( bndTree.shapes().faceLabels()[info.index()] @@ -250,16 +253,26 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) forAll(nearest, facei) { - const label index = nearest[facei].second(); + const auto& near = nearest[facei]; + + const label index = near.second(); if (index == labelMax) { - // Not found on any processor. How to map? + // Not found on any processor, or simply too far. + // How to map? } else if (globalCells.isLocal(index)) { - cellOrFaceLabels[facei] = globalCells.toLocal(index); facesToSubset.set(facei); + + // Mark as special if too far away + cellOrFaceLabels[facei] = + ( + (near.first() < maxDistanceSqr_) + ? globalCells.toLocal(index) + : -1 + ); } } @@ -295,7 +308,7 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) { // With point interpolation - samplePoints_.resize(pointMap.size()); + samplePoints_ = points(); sampleElements_.resize(pointMap.size(), -1); // Store any face per point (without using pointFaces()) @@ -314,32 +327,25 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) if (sampleSource_ == samplingSource::cells) { - // samplePoints_ : per surface point a location inside the cell - // sampleElements_ : per surface point the cell + // sampleElements_ : per surface point, the cell + // samplePoints_ : per surface point, a location inside the cell - forAll(points(), pointi) + forAll(samplePoints_, pointi) { - const point& pt = points()[pointi]; + // Use _copy_ of point + const point pt = samplePoints_[pointi]; const label celli = cellOrFaceLabels[pointToFace[pointi]]; sampleElements_[pointi] = celli; - // Check if point inside cell if ( - mesh().pointInCell - ( - pt, - sampleElements_[pointi], - meshSearcher.decompMode() - ) + celli >= 0 + && !mesh().pointInCell(pt, celli, meshSearcher.decompMode()) ) { - samplePoints_[pointi] = pt; - } - else - { + // Point not inside cell // Find nearest point on faces of cell scalar minDistSqr = VGREAT; @@ -348,7 +354,12 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) { const face& f = mesh().faces()[facei]; - pointHit info = f.nearestPoint(pt, mesh().points()); + pointHit info = + f.nearestPoint + ( + pt, + mesh().points() + ); if (info.distance() < minDistSqr) { @@ -361,50 +372,52 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) } else if (sampleSource_ == samplingSource::insideCells) { - // samplePoints_ : per surface point a location inside the cell // sampleElements_ : per surface point the cell + // samplePoints_ : per surface point a location inside the cell - forAll(points(), pointi) + forAll(samplePoints_, pointi) { - const point& pt = points()[pointi]; - const label celli = cellOrFaceLabels[pointToFace[pointi]]; sampleElements_[pointi] = celli; - samplePoints_[pointi] = pt; } } else // samplingSource::boundaryFaces { - // samplePoints_ : per surface point a location on the boundary - // sampleElements_ : per surface point the boundary face containing + // sampleElements_ : per surface point, the boundary face containing // the location + // samplePoints_ : per surface point, a location on the boundary - forAll(points(), pointi) + forAll(samplePoints_, pointi) { - const point& pt = points()[pointi]; + const point& pt = samplePoints_[pointi]; const label facei = cellOrFaceLabels[pointToFace[pointi]]; sampleElements_[pointi] = facei; - samplePoints_[pointi] = mesh().faces()[facei].nearestPoint - ( - pt, - mesh().points() - ).rawPoint(); + + if (facei >= 0) + { + samplePoints_[pointi] = + mesh().faces()[facei].nearestPoint + ( + pt, + mesh().points() + ).rawPoint(); + } } } } else { // if sampleSource_ == cells: - // sampleElements_ : per surface triangle the cell + // sampleElements_ : per surface face, the cell // samplePoints_ : n/a // if sampleSource_ == insideCells: - // sampleElements_ : -1 or per surface triangle the cell + // sampleElements_ : -1 or per surface face, the cell // samplePoints_ : n/a - // else: - // sampleElements_ : per surface triangle the boundary face + // if sampleSource_ == boundaryFaces: + // sampleElements_ : per surface face, the boundary face // samplePoints_ : n/a sampleElements_.transfer(cellOrFaceLabels); @@ -432,14 +445,14 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) forAll(samplePoints_, pointi) { meshTools::writeOBJ(str, points()[pointi]); - vertI++; + ++vertI; meshTools::writeOBJ(str, samplePoints_[pointi]); - vertI++; + ++vertI; label elemi = sampleElements_[pointi]; meshTools::writeOBJ(str, centres[elemi]); - vertI++; + ++vertI; str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << nl; } } @@ -449,11 +462,11 @@ bool Foam::sampledMeshedSurface::update(const meshSearch& meshSearcher) forAll(sampleElements_, triI) { meshTools::writeOBJ(str, faceCentres()[triI]); - vertI++; + ++vertI; label elemi = sampleElements_[triI]; meshTools::writeOBJ(str, centres[elemi]); - vertI++; + ++vertI; str << "l " << vertI-1 << ' ' << vertI << nl; } } @@ -475,7 +488,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface ) : sampledSurface(name, mesh), - MeshStorage(), + meshedSurface(), surfaceName_(surfaceName), surface_ ( @@ -487,7 +500,9 @@ Foam::sampledMeshedSurface::sampledMeshedSurface keepIds_(true), zoneIds_(), sampleElements_(), - samplePoints_() + samplePoints_(), + maxDistanceSqr_(Foam::sqr(GREAT)), + defaultValues_() {} @@ -499,7 +514,7 @@ Foam::sampledMeshedSurface::sampledMeshedSurface ) : sampledSurface(name, mesh, dict), - MeshStorage(), + meshedSurface(), surfaceName_ ( meshedSurface::findFile @@ -518,8 +533,25 @@ Foam::sampledMeshedSurface::sampledMeshedSurface keepIds_(dict.getOrDefault("keepIds", true)), zoneIds_(), sampleElements_(), - samplePoints_() + samplePoints_(), + maxDistanceSqr_(Foam::sqr(GREAT)), + defaultValues_(dict.subOrEmptyDict("defaultValue")) { + if (dict.readIfPresent("maxDistance", maxDistanceSqr_)) + { + // Info<< "Limit to maxDistance " << maxDistanceSqr_ << nl; + // if (defaultValues_.empty()) + // { + // Info<< "defaultValues = zero" << nl; + // } + // else + // { + // defaultValues_.writeEntry(Info); + // } + + maxDistanceSqr_ = Foam::sqr(maxDistanceSqr_); + } + wordRes includePatches; dict.readIfPresent("patches", includePatches); includePatches.uniq(); @@ -597,7 +629,7 @@ bool Foam::sampledMeshedSurface::expire() } sampledSurface::clearGeom(); - MeshStorage::clear(); + meshedSurface::clear(); zoneIds_.clear(); sampleElements_.clear(); diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H index bfc4ca4c69..5b61d0ef18 100644 --- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H +++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurface.H @@ -35,33 +35,35 @@ Description - 6 different modes: - source=cells, interpolate=false: - finds per triangle centre the \em nearest cell centre and uses its value + finds per surface face centre the \em nearest cell centre + and uses its value - source=cells, interpolate=true: - finds per triangle centre the \em nearest cell centre. + finds per surface face centre the \em nearest cell centre. Per surface point checks if this nearest cell is the one containing point; otherwise projects the point onto the nearest point on the boundary of the cell (to make sure interpolateCellPoint gets a valid location) - source=insideCells, interpolate=false: - finds per triangle centre the cell containing it and uses its value. - Trims triangles outside mesh. + finds per surface face centre the cell containing it + and uses its value. + Trims surface faces outside of the mesh. - source=insideCells, interpolate=true: Per surface point interpolate cell containing it. - source=boundaryFaces, interpolate=false: - finds per triangle centre the nearest point on the boundary + finds per surface face centre the \em nearest point on the boundary (uncoupled faces only) and uses the value (or 0 if the nearest is on an empty boundary) - source=boundaryFaces, interpolate=true: - finds per triangle centre the nearest point on the boundary + finds per surface face centre the \em nearest point on the boundary (uncoupled faces only). Per surface point projects the point onto this boundary face (to make sure interpolateCellPoint gets a valid location) - - since it finds a nearest per triangle each triangle is guaranteed - to be on one processor only. So after stitching (by sampledSurfaces) - the original surface should be complete. + - since it finds the nearest per surface face, each surface face + is guaranteed to be on one processor only. + So after stitching the original surface should be complete. This is often embedded as part of a sampled surfaces function object. @@ -75,6 +77,16 @@ Usage type meshedSurface; surface something.obj; source cells; + + //- Max sampling distance + maxDistance 0.005; + + //- Fallback for missed sampling in 'cells' mode + defaultValue + { + "p.*" 1e5; + T 273.15; + } } } \endverbatim @@ -90,8 +102,11 @@ Usage file | Alternative file name | no | fileType | The surface format | no | (extension) scale | Surface scaling factor | no | 0 + maxDistance | Max search distance | no | GREAT + defaultValue | Value beyond max distance (dictionary) | no | empty \endtable + SourceFiles sampledMeshedSurface.C sampledMeshedSurfaceTemplates.C @@ -102,8 +117,7 @@ SourceFiles #define sampledMeshedSurface_H #include "sampledSurface.H" -#include "MeshedSurface.H" -#include "MeshedSurfacesFwd.H" +#include "MeshedSurfaces.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -123,18 +137,19 @@ class sampledMeshedSurface public meshedSurface { public: + //- Types of sampling regions - enum class samplingSource + enum class samplingSource : unsigned char { cells, //!< Use nearest cell value - insideCells, //!< Face within a cell - boundaryFaces //!< Use boundary face values + insideCells, //!< Surface face within a cell, or trim + boundaryFaces //!< Use nearest boundary face values }; private: - //- Private typedefs for convenience - typedef meshedSurface MeshStorage; + //- The concrete storage type for faces and points + typedef meshedSurface Mesh; // Private Data @@ -165,6 +180,12 @@ private: //- Local points to sample per point pointField samplePoints_; + //- Max search distance squared + scalar maxDistanceSqr_; + + //- Out-of-range fallback values (when beyond the max search distance) + dictionary defaultValues_; + // Private Member Functions @@ -239,13 +260,13 @@ public: //- Points of surface virtual const pointField& points() const { - return MeshStorage::points(); + return Mesh::points(); } //- Faces of surface virtual const faceList& faces() const { - return MeshStorage::surfFaces(); + return Mesh::surfFaces(); } //- Per-face zone/region information @@ -257,31 +278,31 @@ public: //- Per-face identifier (eg, element Id) virtual const labelList& faceIds() const { - return MeshStorage::faceIds(); + return Mesh::faceIds(); } //- Face area vectors virtual const vectorField& Sf() const { - return MeshStorage::Sf(); + return Mesh::Sf(); } //- Face area magnitudes virtual const scalarField& magSf() const { - return MeshStorage::magSf(); + return Mesh::magSf(); } //- Face centres virtual const vectorField& Cf() const { - return MeshStorage::Cf(); + return Mesh::Cf(); } //- If element ids/order of the original surface are kept virtual bool hasFaceIds() const { - return keepIds_ && !MeshStorage::faceIds().empty(); + return keepIds_ && !Mesh::faceIds().empty(); } //- Sampling boundary values instead of cell values diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C index dd8cd7c259..6833dd2bda 100644 --- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C +++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C @@ -37,8 +37,14 @@ Foam::sampledMeshedSurface::sampleOnFaces const interpolation& sampler ) const { + const Type deflt + ( + defaultValues_.getOrDefault(sampler.psi().name(), Zero) + ); + const labelList& elements = sampleElements_; + if (!onBoundary()) { // Sample cells @@ -48,7 +54,8 @@ Foam::sampledMeshedSurface::sampleOnFaces sampler, elements, faces(), - points() + points(), + deflt ); } @@ -60,33 +67,33 @@ Foam::sampledMeshedSurface::sampleOnFaces auto tvalues = tmp>::New(elements.size()); auto& values = tvalues.ref(); - const polyBoundaryMesh& pbm = mesh().boundaryMesh(); - const label nBnd = mesh().nBoundaryFaces(); - // Create flat boundary field - Field bVals(nBnd, Zero); + const polyBoundaryMesh& pbm = mesh().boundaryMesh(); + + Field bVals(mesh().nBoundaryFaces(), Zero); const auto& bField = sampler.psi().boundaryField(); forAll(bField, patchi) { - const label bFacei = (pbm[patchi].start() - mesh().nInternalFaces()); - - SubList - ( - bVals, - bField[patchi].size(), - bFacei - ) = bField[patchi]; + SubList(bVals, pbm[patchi].range()) = bField[patchi]; } - // Sample in flat boundary field + // Sample within the flat boundary field forAll(elements, i) { const label bFacei = (elements[i] - mesh().nInternalFaces()); - values[i] = bVals[bFacei]; + + if (bFacei < 0) + { + values[i] = deflt; + } + else + { + values[i] = bVals[bFacei]; + } } return tvalues; @@ -100,34 +107,60 @@ Foam::sampledMeshedSurface::sampleOnPoints const interpolation& interpolator ) const { - // One value per vertex - auto tvalues = tmp>::New(sampleElements_.size()); + const Type deflt + ( + defaultValues_.getOrDefault(interpolator.psi().name(), Zero) + ); + + const labelList& elements = sampleElements_; + + + // One value per vertex. + // - sampleElements_ and samplePoints_ have the identical size + auto tvalues = tmp>::New(elements.size()); auto& values = tvalues.ref(); if (!onBoundary()) { // Sample cells - forAll(sampleElements_, pointi) + forAll(elements, i) { - values[pointi] = interpolator.interpolate - ( - samplePoints_[pointi], - sampleElements_[pointi] - ); + const label celli = elements[i]; + + if (celli < 0) + { + values[i] = deflt; + } + else + { + values[i] = interpolator.interpolate + ( + samplePoints_[i], + celli + ); + } } } - else + + + // + // Sample boundary faces + // + + forAll(elements, i) { - // Sample boundary faces + const label facei = elements[i]; - forAll(samplePoints_, pointi) + if (facei < 0) { - const label facei = sampleElements_[pointi]; - - values[pointi] = interpolator.interpolate + values[i] = deflt; + } + else + { + values[i] = interpolator.interpolate ( - samplePoints_[pointi], + samplePoints_[i], mesh().faceOwner()[facei], facei ); diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H index b21df8349e..d135ab0693 100644 --- a/src/sampling/sampledSurface/sampledSurface/sampledSurface.H +++ b/src/sampling/sampledSurface/sampledSurface/sampledSurface.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2016-2019 OpenCFD Ltd. + Copyright (C) 2016-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -132,14 +132,16 @@ protected: // Protected Member Functions - //- General loop for sampling elements to faces + //- General loop for sampling elements to faces. + // The defaultValue is used for invalid (negative) elements template static tmp> sampleOnFaces ( const interpolation& sampler, const labelUList& elements, const faceList& fcs, - const pointField& pts + const pointField& pts, + const Type& defaultValue = Type(Zero) ); diff --git a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C index 90a5e7f376..2b494d5027 100644 --- a/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C +++ b/src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2018-2019 OpenCFD Ltd. + Copyright (C) 2018-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -37,7 +37,8 @@ Foam::sampledSurface::sampleOnFaces const interpolation& sampler, const labelUList& elements, const faceList& fcs, - const pointField& pts + const pointField& pts, + const Type& defaultValue ) { const label len = elements.size(); @@ -57,9 +58,16 @@ Foam::sampledSurface::sampleOnFaces for (label i=0; i < len; ++i) { const label celli = elements[i]; - const point pt = fcs[i].centre(pts); + if (celli < 0) + { + values[i] = defaultValue; + } + else + { + const point pt = fcs[i].centre(pts); - values[i] = sampler.interpolate(pt, celli); + values[i] = sampler.interpolate(pt, celli); + } } return tvalues; diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling index b5f9317229..529a5a6c06 100644 --- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling +++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/sampling @@ -121,6 +121,44 @@ sampled } +// Sample and write for meshed surfaces +sampleMesh +{ + type surfaces; + libs (sampling); + log true; + + writeControl writeTime; + writeInterval 1; + + surfaceFormat vtk; + + fields (p rho U T); + + surfaces + { + // Oversized sampling - for general testing + meshed_sample + { + type meshedSurface; + surface oversized_sample.obj; + source cells; + maxDistance 0.0025; + defaultValue + { + T 273; + } + } + + meshed_interpolate + { + $meshed_sample; + interpolate true; + } + } +} + + // * * * * * * * * * * * * * * * Calculations * * * * * * * * * * * * * * * // massflow