From 25adff2503ceedde854075ee5bfa9215b5f20f36 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Tue, 30 Jan 2018 10:15:19 +0100 Subject: [PATCH] ENH: improved handling of bounded sampled planes (issue #714) - now warn about the following: * the bounding box does not overlap wih the global mesh * plane does not intersect the (valid) bounding box * plane does not intersect the global mesh - add bounding to the "plane" variant of a sampled plane. --- .../sampledCuttingPlane/sampledCuttingPlane.C | 48 +++++++-- .../sampledCuttingPlane/sampledCuttingPlane.H | 2 +- .../sampledPlane/sampledPlane.C | 101 +++++++++++++++--- .../sampledPlane/sampledPlane.H | 7 +- .../sampledPlane/sampledPlaneTemplates.C | 17 +-- .../surface/cuttingPlane/cuttingPlane.C | 82 +++++++------- .../surface/cuttingPlane/cuttingPlane.H | 11 +- 7 files changed, 185 insertions(+), 83 deletions(-) diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C index ef1217b7ed..1d5eaaf38a 100644 --- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C +++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C @@ -114,7 +114,7 @@ void Foam::sampledCuttingPlane::createGeometry() false ), mesh, - dimensionedScalar("zero", dimLength, 0) + dimensionedScalar(dimLength) ) ); volScalarField& cellDistance = cellDistancePtr_(); @@ -212,7 +212,7 @@ void Foam::sampledCuttingPlane::createGeometry() false ), pointMesh::New(mesh), - dimensionedScalar("zero", dimLength, 0) + dimensionedScalar(dimLength) ); pDist.primitiveFieldRef() = pointDistance_; @@ -220,7 +220,6 @@ void Foam::sampledCuttingPlane::createGeometry() pDist.write(); } - //- Direct from cell field and point field. isoSurfPtr_.reset ( @@ -244,6 +243,43 @@ void Foam::sampledCuttingPlane::createGeometry() //) ); + // Verify specified bounding box + if (!bounds_.empty()) + { + // Bounding box does not overlap with (global) mesh! + if (!bounds_.overlaps(mesh.bounds())) + { + WarningInFunction + << nl + << name() << " : " + << "Bounds " << bounds_ + << " do not overlap the mesh bounding box " << mesh.bounds() + << nl << endl; + } + + // Plane does not intersect the bounding box + if (!bounds_.intersects(plane_)) + { + WarningInFunction + << nl + << name() << " : " + << "Plane "<< plane_ << " does not intersect the bounds " + << bounds_ + << nl << endl; + } + } + + // Plane does not intersect the (global) mesh! + if (!mesh.bounds().intersects(plane_)) + { + WarningInFunction + << nl + << name() << " : " + << "Plane "<< plane_ << " does not intersect the mesh bounds " + << mesh.bounds() + << nl << endl; + } + if (debug) { print(Pout); @@ -297,12 +333,6 @@ Foam::sampledCuttingPlane::sampledCuttingPlane } -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::sampledCuttingPlane::~sampledCuttingPlane() -{} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // bool Foam::sampledCuttingPlane::needsUpdate() const diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H index 05e9a2247e..756ee7ba9d 100644 --- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H +++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H @@ -133,7 +133,7 @@ public: //- Destructor - virtual ~sampledCuttingPlane(); + virtual ~sampledCuttingPlane() = default; // Member Functions diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C index 333bb94527..ed7babbb00 100644 --- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C +++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C @@ -52,12 +52,13 @@ Foam::sampledPlane::sampledPlane sampledSurface(name, mesh), cuttingPlane(planeDesc), zoneKey_(zoneKey), + bounds_(), triangulate_(triangulate), needsUpdate_(true) { if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1) { - Info<< "cellZone " << zoneKey_ + Info<< "cellZone(s) " << zoneKey_ << " not found - using entire mesh" << endl; } } @@ -72,7 +73,8 @@ Foam::sampledPlane::sampledPlane : sampledSurface(name, mesh, dict), cuttingPlane(plane(dict)), - zoneKey_(keyType::null), + zoneKey_(dict.lookupOrDefault("zone", keyType::null)), + bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)), triangulate_(dict.lookupOrDefault("triangulate", true)), needsUpdate_(true) { @@ -82,29 +84,21 @@ Foam::sampledPlane::sampledPlane { coordinateSystem cs(mesh, dict.subDict("coordinateSystem")); - point base = cs.globalPosition(planeDesc().refPoint()); - vector norm = cs.globalVector(planeDesc().normal()); + const point base = cs.globalPosition(planeDesc().refPoint()); + const vector norm = cs.globalVector(planeDesc().normal()); // Assign the plane description static_cast(*this) = plane(base, norm); } - dict.readIfPresent("zone", zoneKey_); - if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1) { - Info<< "cellZone " << zoneKey_ + Info<< "cellZone(s) " << zoneKey_ << " not found - using entire mesh" << endl; } } -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::sampledPlane::~sampledPlane() -{} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // bool Foam::sampledPlane::needsUpdate() const @@ -137,9 +131,88 @@ bool Foam::sampledPlane::update() sampledSurface::clearGeom(); + const plane& pln = static_cast(*this); + + // Verify specified bounding box + if (!bounds_.empty()) + { + // Bounding box does not overlap with (global) mesh! + if (!bounds_.overlaps(mesh().bounds())) + { + WarningInFunction + << nl + << name() << " : " + << "Bounds " << bounds_ + << " do not overlap the mesh bounding box " << mesh().bounds() + << nl << endl; + } + + // Plane does not intersect the bounding box + if (!bounds_.intersects(pln)) + { + WarningInFunction + << nl + << name() << " : " + << "Plane "<< pln << " does not intersect the bounds " + << bounds_ + << nl << endl; + } + } + + // Plane does not intersect the (global) mesh! + if (!mesh().bounds().intersects(pln)) + { + WarningInFunction + << nl + << name() << " : " + << "Plane "<< pln << " does not intersect the mesh bounds " + << mesh().bounds() + << nl << endl; + } + labelList selectedCells = mesh().cellZones().findMatching(zoneKey_).used(); - if (returnReduce(selectedCells.empty(), andOp())) + bool fullMesh = returnReduce(selectedCells.empty(), andOp()); + + if (!bounds_.empty()) + { + const auto& cellCentres = static_cast(mesh()).C(); + + if (fullMesh) + { + const label len = mesh().nCells(); + + selectedCells.setSize(len); + + label count = 0; + for (label celli=0; celli < len; ++celli) + { + if (bounds_.contains(cellCentres[celli])) + { + selectedCells[count++] = celli; + } + } + + selectedCells.setSize(count); + } + else + { + label count = 0; + for (const label celli : selectedCells) + { + if (bounds_.contains(cellCentres[celli])) + { + selectedCells[count++] = celli; + } + } + + selectedCells.setSize(count); + } + + fullMesh = false; + } + + if (fullMesh) { reCut(mesh(), triangulate_); } diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.H b/src/sampling/sampledSurface/sampledPlane/sampledPlane.H index 8c8867070a..8a64dfe896 100644 --- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.H +++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.H @@ -59,7 +59,10 @@ class sampledPlane // Private data //- If restricted to zones, name of this zone or a regular expression - keyType zoneKey_; + const keyType zoneKey_; + + //- Optional bounding box to trim triangles against + const boundBox bounds_; //- Triangulated faces or keep faces as is const bool triangulate_; @@ -110,7 +113,7 @@ public: //- Destructor - virtual ~sampledPlane(); + virtual ~sampledPlane() = default; // Member Functions diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C b/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C index f3ff590c1b..f081e27d74 100644 --- a/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C +++ b/src/sampling/sampledSurface/sampledPlane/sampledPlaneTemplates.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -45,8 +45,10 @@ Foam::sampledPlane::interpolateField const interpolation& interpolator ) const { - // One value per point - tmp> tvalues(new Field(points().size())); + // One value per point. + // Initialize with Zero to handle missed/degenerate faces + + tmp> tvalues(new Field(points().size(), Zero)); Field& values = tvalues.ref(); boolList pointDone(points().size(), false); @@ -54,19 +56,18 @@ Foam::sampledPlane::interpolateField forAll(faces(), cutFacei) { const face& f = faces()[cutFacei]; + const label celli = meshCells()[cutFacei]; - forAll(f, faceVertI) + for (const label pointi : f) { - label pointi = f[faceVertI]; - if (!pointDone[pointi]) { + pointDone[pointi] = true; values[pointi] = interpolator.interpolate ( points()[pointi], - meshCells()[cutFacei] + celli ); - pointDone[pointi] = true; } } } diff --git a/src/sampling/surface/cuttingPlane/cuttingPlane.C b/src/sampling/surface/cuttingPlane/cuttingPlane.C index be5bf0dce6..d48ab71761 100644 --- a/src/sampling/surface/cuttingPlane/cuttingPlane.C +++ b/src/sampling/surface/cuttingPlane/cuttingPlane.C @@ -33,8 +33,8 @@ License // Set values for what is close to zero and what is considered to // be positive (and not just rounding noise) //! \cond localScope -const Foam::scalar zeroish = Foam::SMALL; -const Foam::scalar positive = Foam::SMALL * 1E3; +static const Foam::scalar zeroish = Foam::SMALL; +static const Foam::scalar positive = Foam::SMALL * 1E3; //! \endcond // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -49,31 +49,25 @@ void Foam::cuttingPlane::calcCutCells const labelListList& cellEdges = mesh.cellEdges(); const edgeList& edges = mesh.edges(); - label listSize = cellEdges.size(); - if (notNull(cellIdLabels)) - { - listSize = cellIdLabels.size(); - } + const label len = + (notNull(cellIdLabels) ? cellIdLabels.size() : cellEdges.size()); - meshCells_.setSize(listSize); + meshCells_.setSize(len); label cutcelli(0); // Find the cut cells by detecting any cell that uses points with // opposing dotProducts. - for (label listI = 0; listI < listSize; ++listI) + for (label listi = 0; listi < len; ++listi) { - label celli = listI; - - if (notNull(cellIdLabels)) - { - celli = cellIdLabels[listI]; - } + const label celli = + (notNull(cellIdLabels) ? cellIdLabels[listi] : listi); const labelList& cEdges = cellEdges[celli]; label nCutEdges = 0; - forAll(cEdges, i) + + for (const label edgei : cEdges) { - const edge& e = edges[cEdges[i]]; + const edge& e = edges[edgei]; if ( @@ -81,9 +75,7 @@ void Foam::cuttingPlane::calcCutCells || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive) ) { - nCutEdges++; - - if (nCutEdges > 2) + if (++nCutEdges > 2) { meshCells_[cutcelli++] = celli; break; @@ -129,7 +121,7 @@ void Foam::cuttingPlane::intersectEdges const point& p0 = points[e[0]]; const point& p1 = points[e[1]]; - scalar alpha = lineIntersect(linePointRef(p0, p1)); + const scalar alpha = lineIntersect(linePointRef(p0, p1)); if (alpha < zeroish) { @@ -186,10 +178,8 @@ bool Foam::cuttingPlane::walkCell // If so should e.g. decompose the cells on both faces and redo // the calculation. - forAll(fEdges, i) + for (const label edge2I : fEdges) { - label edge2I = fEdges[i]; - if (edge2I != edgeI && edgePoint[edge2I] != -1) { nextEdgeI = edge2I; @@ -212,7 +202,7 @@ bool Foam::cuttingPlane::walkCell edgeI = nextEdgeI; - nIter++; + ++nIter; if (nIter > 1000) { @@ -232,16 +222,14 @@ bool Foam::cuttingPlane::walkCell { return true; } - else - { - WarningInFunction - << "Did not find closed walk along surface of cell " << celli - << " starting from edge " << startEdgeI << nl - << "Collected cutPoints so far:" << faceVerts - << endl; - return false; - } + WarningInFunction + << "Did not find closed walk along surface of cell " << celli + << " starting from edge " << startEdgeI << nl + << "Collected cutPoints so far:" << faceVerts + << endl; + + return false; } @@ -261,19 +249,15 @@ void Foam::cuttingPlane::walkCellCuts // scratch space for calculating the face vertices DynamicList