diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C index 08a2d84f9f..69f971030e 100644 --- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C @@ -108,7 +108,7 @@ void Foam::sampledIsoSurface::getIsoFields() const fvm ) ); - volFieldPtr_ = storedVolFieldPtr_.operator->(); + volFieldPtr_ = storedVolFieldPtr_.get(); } else { @@ -132,7 +132,7 @@ void Foam::sampledIsoSurface::getIsoFields() const // (volPointInterpolation::interpolate with cache=false deletes any // registered one or if mesh.changing()) - if (!subMeshPtr_.valid()) + if (subMeshPtr_.empty()) { const word pointFldName = "volPointInterpolate_" @@ -200,8 +200,8 @@ void Foam::sampledIsoSurface::getIsoFields() const } - // If averaging redo the volField. Can only be done now since needs the - // point field. + // If averaging redo the volField. + // Can only be done now since needs the point field. if (average_) { storedVolFieldPtr_.reset @@ -261,7 +261,7 @@ void Foam::sampledIsoSurface::getIsoFields() const } - // Pointfield on submesh + // The point field on subMesh const word pointFldName = "volPointInterpolate_" @@ -360,30 +360,30 @@ bool Foam::sampledIsoSurface::updateGeometry() const return false; } - // Get any subMesh - if (zoneID_.index() != -1 && !subMeshPtr_.valid()) + // Get sub-mesh if any + if + ( + (-1 != mesh().cellZones().findIndex(zoneNames_)) + && subMeshPtr_.empty() + ) { const polyBoundaryMesh& patches = mesh().boundaryMesh(); // Patch to put exposed internal faces into const label exposedPatchi = patches.findPatchID(exposedPatchName_); - if (debug) - { - Info<< "Allocating subset of size " - << mesh().cellZones()[zoneID_.index()].size() - << " with exposed faces into patch " - << patches[exposedPatchi].name() << endl; - } + DebugInfo + << "Allocating subset of size " + << mesh().cellZones().selection(zoneNames_).count() + << " with exposed faces into patch " + << patches[exposedPatchi].name() << endl; - // Remove old mesh if any - subMeshPtr_.clear(); subMeshPtr_.reset ( new fvMeshSubset ( fvm, - mesh().cellZones()[zoneID_.index()], + mesh().cellZones().selection(zoneNames_), exposedPatchi ) ); @@ -474,8 +474,8 @@ Foam::sampledIsoSurface::sampledIsoSurface mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)), regularise_(dict.lookupOrDefault("regularise", true)), average_(dict.lookupOrDefault("average", false)), - zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()), - exposedPatchName_(word::null), + zoneNames_(), + exposedPatchName_(), surfPtr_(nullptr), prevTimeIndex_(-1), storedVolFieldPtr_(nullptr), @@ -491,27 +491,30 @@ Foam::sampledIsoSurface::sampledIsoSurface << " span across cells." << exit(FatalIOError); } - if (zoneID_.index() != -1) + + if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone")) + { + zoneNames_.resize(1); + dict.readEntry("zone", zoneNames_.first()); + } + + if (-1 != mesh.cellZones().findIndex(zoneNames_)) { dict.readEntry("exposedPatchName", exposedPatchName_); if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1) { - FatalIOErrorInFunction - ( - dict - ) << "Cannot find patch " << exposedPatchName_ + FatalIOErrorInFunction(dict) + << "Cannot find patch " << exposedPatchName_ << " in which to put exposed faces." << endl << "Valid patches are " << mesh.boundaryMesh().names() << exit(FatalIOError); } - if (debug && zoneID_.index() != -1) - { - Info<< "Restricting to cellZone " << zoneID_.name() - << " with exposed internal faces into patch " - << exposedPatchName_ << endl; - } + DebugInfo + << "Restricting to cellZone(s) " << flatOutput(zoneNames_) + << " with exposed internal faces into patch " + << exposedPatchName_ << endl; } } diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H index 9ecc52b9cc..beb4503242 100644 --- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H @@ -54,7 +54,8 @@ Usage regularise | point snapping | yes | average | cell values from averaged point values | no | false bounds | limit with bounding box | no | - zone | limit to specified zone | no | + zone | limit to cell zone (name or regex) | no | + zones | limit to cell zones (names, regexs) | no | exposedPatchName | name for zone subset | partly | \endtable @@ -68,7 +69,6 @@ SourceFiles #include "isoSurface.H" #include "sampledSurface.H" -#include "ZoneIDs.H" #include "fvMeshSubset.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -104,8 +104,8 @@ class sampledIsoSurface //- Whether to recalculate cell values as average of point values const bool average_; - //- Zone name/index (if restricted to zones) - mutable cellZoneID zoneID_; + //- The zone or zones for the iso-surface + wordRes zoneNames_; //- For zones: patch to put exposed faces into mutable word exposedPatchName_; diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C index cd869ea898..f51b7f1ccf 100644 --- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C @@ -201,7 +201,6 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)), regularise_(dict.lookupOrDefault("regularise", true)), average_(dict.lookupOrDefault("average", true)), - zoneKey_(keyType::null), prevTimeIndex_(-1), meshCells_() {} diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H index a197649e5d..93e2e4a3a9 100644 --- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H @@ -56,6 +56,9 @@ Usage bounds | limit with bounding box | no | \endtable +Note + Does not currently support cell zones. + SourceFiles sampledIsoSurfaceCell.C @@ -102,9 +105,6 @@ class sampledIsoSurfaceCell //- Whether to recalculate cell values as average of point values const bool average_; - //- If restricted to zones, name of this zone or a regular expression - keyType zoneKey_; - // Recreated for every isoSurface diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C index 6dc7894ac6..661dc7b31f 100644 --- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C +++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C @@ -46,6 +46,51 @@ namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +void Foam::sampledCuttingPlane::checkBoundsIntersection +( + const plane& pln, + const boundBox& meshBb +) const +{ + // Verify specified bounding box + if (!bounds_.empty()) + { + // Bounding box does not overlap with (global) mesh! + if (!bounds_.overlaps(meshBb)) + { + WarningInFunction + << nl + << name() << " : " + << "Bounds " << bounds_ + << " do not overlap the mesh bounding box " << meshBb + << 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 (!meshBb.intersects(pln)) + { + WarningInFunction + << nl + << name() << " : " + << "Plane "<< pln << " does not intersect the mesh bounds " + << meshBb + << nl << endl; + } +} + + void Foam::sampledCuttingPlane::createGeometry() { if (debug) @@ -62,29 +107,51 @@ void Foam::sampledCuttingPlane::createGeometry() // Clear derived data clearGeom(); - // Get any subMesh - if (zoneID_.index() != -1 && !subMeshPtr_.valid()) + const fvMesh& fvm = static_cast(this->mesh()); + + // Get sub-mesh if any + if + ( + (-1 != mesh().cellZones().findIndex(zoneNames_)) + && subMeshPtr_.empty() + ) { const polyBoundaryMesh& patches = mesh().boundaryMesh(); // Patch to put exposed internal faces into const label exposedPatchi = patches.findPatchID(exposedPatchName_); + bitSet cellsToSelect = mesh().cellZones().selection(zoneNames_); + DebugInfo << "Allocating subset of size " - << mesh().cellZones()[zoneID_.index()].size() + << cellsToSelect.count() << " with exposed faces into patch " << patches[exposedPatchi].name() << endl; - subMeshPtr_.reset - ( - new fvMeshSubset - ( - static_cast(mesh()), - mesh().cellZones()[zoneID_.index()], - exposedPatchi - ) - ); + + // If we will use a fvMeshSubset so can apply bounds as well to make + // the initial selection smaller. + if (!bounds_.empty() && cellsToSelect.any()) + { + const auto& cellCentres = fvm.C(); + + for (const label celli : cellsToSelect) + { + const point& cc = cellCentres[celli]; + + if (!bounds_.contains(cc)) + { + cellsToSelect.unset(celli); + } + } + + DebugInfo + << "Bounded subset of size " + << cellsToSelect.count() << endl; + } + + subMeshPtr_.reset(new fvMeshSubset(fvm, cellsToSelect, exposedPatchi)); } @@ -93,9 +160,11 @@ void Foam::sampledCuttingPlane::createGeometry() ( subMeshPtr_.valid() ? subMeshPtr_().subMesh() - : static_cast(this->mesh()) + : fvm ); + checkBoundsIntersection(plane_, mesh.bounds()); + // Distance to cell centres // ~~~~~~~~~~~~~~~~~~~~~~~~ @@ -121,7 +190,7 @@ void Foam::sampledCuttingPlane::createGeometry() // Internal field { - const pointField& cc = mesh.cellCentres(); + const auto& cc = mesh.cellCentres(); scalarField& fld = cellDistance.primitiveFieldRef(); forAll(cc, i) @@ -167,6 +236,7 @@ void Foam::sampledCuttingPlane::createGeometry() } else { + // Other side cell centres? const pointField& cc = mesh.C().boundaryField()[patchi]; fvPatchScalarField& fld = cellDistanceBf[patchi]; @@ -242,43 +312,6 @@ 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); @@ -302,32 +335,36 @@ Foam::sampledCuttingPlane::sampledCuttingPlane mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)), regularise_(dict.lookupOrDefault("regularise", true)), average_(dict.lookupOrDefault("average", false)), - zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()), - exposedPatchName_(word::null), + zoneNames_(), + exposedPatchName_(), needsUpdate_(true), subMeshPtr_(nullptr), cellDistancePtr_(nullptr), isoSurfPtr_(nullptr) { - if (zoneID_.index() != -1) + if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone")) + { + zoneNames_.resize(1); + dict.readEntry("zone", zoneNames_.first()); + } + + if (-1 != mesh.cellZones().findIndex(zoneNames_)) { dict.readEntry("exposedPatchName", exposedPatchName_); - if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1) + if (-1 == mesh.boundaryMesh().findPatchID(exposedPatchName_)) { - FatalErrorInFunction + FatalIOErrorInFunction(dict) << "Cannot find patch " << exposedPatchName_ << " in which to put exposed faces." << endl << "Valid patches are " << mesh.boundaryMesh().names() << exit(FatalError); } - if (debug && zoneID_.index() != -1) - { - Info<< "Restricting to cellZone " << zoneID_.name() - << " with exposed internal faces into patch " - << exposedPatchName_ << endl; - } + DebugInfo + << "Restricting to cellZone(s) " << flatOutput(zoneNames_) + << " with exposed internal faces into patch " + << exposedPatchName_ << endl; } } diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H index 3f67de6f5e..45b94674b0 100644 --- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H +++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H @@ -55,10 +55,14 @@ Usage mergeTol | tolerance for merging points | no | 1e-6 regularise | point snapping | no | true bounds | limit with bounding box | no | - zone | limit to specified zone | no | + zone | limit to cell zone (name or regex) | no | + zones | limit to cell zones (names, regexs) | no | exposedPatchName | name for zone subset | partly | \endtable +Note + The keyword \c zones has priority over \c zone. + SeeAlso Foam::plane @@ -74,7 +78,6 @@ SourceFiles #include "isoSurface.H" //#include "isoSurfaceCell.H" #include "plane.H" -#include "ZoneIDs.H" #include "fvMeshSubset.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -107,8 +110,8 @@ class sampledCuttingPlane //- Whether to recalculate cell values as average of point values const bool average_; - //- Zone name/index (if restricted to zones) - mutable cellZoneID zoneID_; + //- The zone or zones in which cutting is to occur + wordRes zoneNames_; //- For zones: patch to put exposed faces into mutable word exposedPatchName_; @@ -117,7 +120,7 @@ class sampledCuttingPlane mutable bool needsUpdate_; - //- Optional subsetted mesh + //- Mesh subset (optional: only used with zones) autoPtr subMeshPtr_; //- Distance to cell centres @@ -133,6 +136,13 @@ class sampledCuttingPlane // Private Member Functions + //- Check and warn if bounding box does not intersect mesh or plane + void checkBoundsIntersection + ( + const plane& pln, + const boundBox& meshBb + ) const; + //- Create iso surface void createGeometry(); diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C index 8d2b895ed9..ad9725dcf5 100644 --- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C +++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C @@ -44,6 +44,159 @@ namespace Foam ); } + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::sampledPlane::checkBoundsIntersection +( + const plane& pln, + const boundBox& meshBb +) const +{ + // Verify specified bounding box + if (!bounds_.empty()) + { + // Bounding box does not overlap with (global) mesh! + if (!bounds_.overlaps(meshBb)) + { + WarningInFunction + << nl + << name() << " : " + << "Bounds " << bounds_ + << " do not overlap the mesh bounding box " << meshBb + << 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 (!meshBb.intersects(pln)) + { + WarningInFunction + << nl + << name() << " : " + << "Plane "<< pln << " does not intersect the mesh bounds " + << meshBb + << nl << endl; + } +} + + +Foam::bitSet Foam::sampledPlane::cellSelection +( + const bool warnIntersect +) const +{ + // Zones requested and in use? + const bool hasZones = + returnReduce + ( + (-1 != mesh().cellZones().findIndex(zoneNames_)), + andOp() + ); + + + bitSet cellsToSelect; + + if (hasZones) + { + cellsToSelect = mesh().cellZones().selection(zoneNames_); + } + + + // Subset the zoned cells with the bounds_. + // For a full mesh, use the bounds to define the cell selection. + + // If there are zones cells, use them to build the effective mesh + // bound box. + // Note that for convenience we use cell centres here instead of the + // cell points, since it will only be used for checking. + + boundBox meshBb; + + if (bounds_.empty()) + { + // No bounds restriction, but will need the effective mesh boundBox + // for checking intersections + + if (hasZones && warnIntersect) + { + const auto& cellCentres = static_cast(mesh()).C(); + + for (const label celli : cellsToSelect) + { + const point& cc = cellCentres[celli]; + + meshBb.add(cc); + } + + meshBb.reduce(); + } + else + { + meshBb = mesh().bounds(); // use the regular mesh bounding box + } + } + else + { + const auto& cellCentres = static_cast(mesh()).C(); + + // Only examine cells already set + if (hasZones) + { + for (const label celli : cellsToSelect) + { + const point& cc = cellCentres[celli]; + + meshBb.add(cc); + + if (!bounds_.contains(cc)) + { + cellsToSelect.unset(celli); + } + } + + meshBb.reduce(); + } + else + { + const label len = mesh().nCells(); + + cellsToSelect.resize(len); + + for (label celli=0; celli < len; ++celli) + { + const point& cc = cellCentres[celli]; + + if (bounds_.contains(cc)) + { + cellsToSelect.set(celli); + } + } + + meshBb = mesh().bounds(); // use the regular mesh bounding box + } + } + + if (warnIntersect) + { + checkBoundsIntersection(*this, meshBb); + } + + return cellsToSelect; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::sampledPlane::sampledPlane @@ -51,21 +204,29 @@ Foam::sampledPlane::sampledPlane const word& name, const polyMesh& mesh, const plane& planeDesc, - const keyType& zoneKey, + const wordRes& zones, const bool triangulate ) : sampledSurface(name, mesh), cuttingPlane(planeDesc), - zoneKey_(zoneKey), + zoneNames_(zones), bounds_(), triangulate_(triangulate), needsUpdate_(true) { - if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1) + if (debug) { - Info<< "cellZone(s) " << zoneKey_ - << " not found - using entire mesh" << endl; + if (!zoneNames_.empty()) + { + Info<< " cellZones " << flatOutput(zoneNames_); + + if (-1 == mesh.cellZones().findIndex(zoneNames_)) + { + Info<< " not found!"; + } + Info<< endl; + } } } @@ -79,28 +240,61 @@ Foam::sampledPlane::sampledPlane : sampledSurface(name, mesh, dict), cuttingPlane(plane(dict)), - zoneKey_(dict.lookupOrDefault("zone", keyType::null)), + zoneNames_(), bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)), triangulate_(dict.lookupOrDefault("triangulate", true)), needsUpdate_(true) { + if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone")) + { + zoneNames_.resize(1); + dict.readEntry("zone", zoneNames_.first()); + } + + // Make plane relative to the coordinateSystem (Cartesian) // allow lookup from global coordinate systems if (dict.found("coordinateSystem")) { coordinateSystem cs(mesh, dict.subDict("coordinateSystem")); - const point base = cs.globalPosition(planeDesc().origin()); + const point orig = cs.globalPosition(planeDesc().origin()); const vector norm = cs.globalVector(planeDesc().normal()); + if (debug) + { + Info<< "plane " << name << " :" + << " origin:" << origin() + << " normal:" << normal() + << " defined within a local coordinateSystem" << endl; + } + // Assign the plane description - static_cast(*this) = plane(base, norm); + static_cast(*this) = plane(orig, norm); } - if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1) + + if (debug) { - Info<< "cellZone(s) " << zoneKey_ - << " not found - using entire mesh" << endl; + Info<< "plane " << name << " :" + << " origin:" << origin() + << " normal:" << normal(); + + if (!bounds_.empty()) + { + Info<< " bounds:" << bounds_; + } + + if (!zoneNames_.empty()) + { + Info<< " cellZones " << flatOutput(zoneNames_); + + if (-1 == mesh.cellZones().findIndex(zoneNames_)) + { + Info<< " not found!"; + } + } + Info<< endl; } } @@ -137,99 +331,7 @@ 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().selection(zoneKey_).sortedToc() - ); - - 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_); - } - else - { - reCut(mesh(), triangulate_, selectedCells); - } + performCut(mesh(), triangulate_, this->cellSelection(true)); if (debug) { @@ -334,11 +436,11 @@ Foam::tmp Foam::sampledPlane::interpolate void Foam::sampledPlane::print(Ostream& os) const { os << "sampledPlane: " << name() << " :" - << " base:" << plane::origin() - << " normal:" << plane::normal() - << " triangulate:" << triangulate_ - << " faces:" << faces().size() - << " points:" << points().size(); + << " origin:" << plane::origin() + << " normal:" << plane::normal() + << " triangulate:" << triangulate_ + << " faces:" << faces().size() + << " points:" << points().size(); } diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.H b/src/sampling/sampledSurface/sampledPlane/sampledPlane.H index 731c46fbc5..42657ccb8d 100644 --- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.H +++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.H @@ -54,7 +54,8 @@ Usage planeType | plane description (pointAndNormal etc) | yes | triangulate | triangulate faces | no | true bounds | limit with bounding box | no | - zone | limit to specified zone | no | + zone | limit to cell zone (name or regex) | no | + zones | limit to cell zones (names, regexs) | no | coordinateSystem | plane relative to given coordinate system | no | \endtable @@ -65,6 +66,7 @@ SeeAlso Note Does not actually cut until update() called. + The keyword \c zones has priority over \c zone. SourceFiles sampledPlane.C @@ -93,8 +95,8 @@ class sampledPlane { // Private data - //- If restricted to zones, name of this zone or a regular expression - const keyType zoneKey_; + //- The zone or zones in which cutting is to occur + wordRes zoneNames_; //- Optional bounding box to trim triangles against const boundBox bounds_; @@ -108,6 +110,20 @@ class sampledPlane // Private Member Functions + //- Check and warn if bounding box does not intersect mesh or plane + void checkBoundsIntersection + ( + const plane& pln, + const boundBox& meshBb + ) const; + + //- Define cell selection from zones and bounding box. + // Optionally check and warn if the plane does not intersect + // with the bounds of the mesh (or submesh) or if the bounding box + // does not overlap with the mesh (or submesh) + bitSet cellSelection(const bool warnIntersect=false) const; + + //- Sample volume field onto surface faces template tmp> sampleOnFaces @@ -137,7 +153,7 @@ public: const word& name, const polyMesh& mesh, const plane& planeDesc, - const keyType& zoneKey = word::null, + const wordRes& zones = wordRes(), const bool triangulate = true ); diff --git a/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C b/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C index 0bb7e2b65b..82e158111b 100644 --- a/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C +++ b/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.C @@ -150,10 +150,9 @@ Foam::sampledThresholdCellFaces::sampledThresholdCellFaces fieldName_(dict.get("field")), lowerThreshold_(dict.lookupOrDefault("lowerLimit", -VGREAT)), upperThreshold_(dict.lookupOrDefault("upperLimit", VGREAT)), - zoneKey_(keyType::null), triangulate_(dict.lookupOrDefault("triangulate", false)), prevTimeIndex_(-1), - meshCells_(0) + meshCells_() { if (!dict.found("lowerLimit") && !dict.found("upperLimit")) { diff --git a/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.H b/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.H index 262585543f..bfce213113 100644 --- a/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.H +++ b/src/sampling/sampledSurface/thresholdCellFaces/sampledThresholdCellFaces.H @@ -99,9 +99,6 @@ class sampledThresholdCellFaces //- Threshold value const scalar upperThreshold_; - //- If restricted to zones, name of this zone or a regular expression - keyType zoneKey_; - //- Triangulated faces or keep faces as is bool triangulate_; diff --git a/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C b/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C index 9b05c160c7..744c3adb63 100644 --- a/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C +++ b/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C @@ -45,6 +45,157 @@ namespace Foam } +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::surfMeshSamplePlane::checkBoundsIntersection +( + const plane& pln, + const boundBox& meshBb +) const +{ + // Verify specified bounding box + if (!bounds_.empty()) + { + // Bounding box does not overlap with (global) mesh! + if (!bounds_.overlaps(meshBb)) + { + 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 (!meshBb.intersects(pln)) + { + WarningInFunction + << nl + << name() << " : " + << "Plane "<< pln << " does not intersect the mesh bounds " + << meshBb + << nl << endl; + } +} + + +Foam::bitSet Foam::surfMeshSamplePlane::cellSelection +( + const bool warnIntersect +) const +{ + // Zones requested and in use? + const bool hasZones = + returnReduce + ( + (-1 != mesh().cellZones().findIndex(zoneNames_)), + andOp() + ); + + + bitSet cellsToSelect; + + if (hasZones) + { + cellsToSelect = mesh().cellZones().selection(zoneNames_); + } + + // Subset the zoned cells with the bounds_. + // For a full mesh, use the bounds to define the cell selection. + + // If there are zones cells, use them to build the effective mesh + // bound box. + // Note that for convenience we use cell centres here instead of the + // cell points, since it will only be used for checking. + + boundBox meshBb; + + if (bounds_.empty()) + { + // No bounds restriction, but will need the effective mesh boundBox + // for checking intersections + + if (hasZones && warnIntersect) + { + const auto& cellCentres = static_cast(mesh()).C(); + + for (const label celli : cellsToSelect) + { + const point& cc = cellCentres[celli]; + + meshBb.add(cc); + } + + meshBb.reduce(); + } + else + { + meshBb = mesh().bounds(); // use the regular mesh bounding box + } + } + else + { + const auto& cellCentres = static_cast(mesh()).C(); + + // Only examine cells already set + if (hasZones) + { + for (const label celli : cellsToSelect) + { + const point& cc = cellCentres[celli]; + + meshBb.add(cc); + + if (!bounds_.contains(cc)) + { + cellsToSelect.unset(celli); + } + } + + meshBb.reduce(); + } + else + { + const label len = mesh().nCells(); + + cellsToSelect.resize(len); + + for (label celli=0; celli < len; ++celli) + { + const point& cc = cellCentres[celli]; + + if (bounds_.contains(cc)) + { + cellsToSelect.set(celli); + } + } + + meshBb = mesh().bounds(); // use the regular mesh bounding box + } + } + + if (warnIntersect) + { + checkBoundsIntersection(*this, meshBb); + } + + return cellsToSelect; +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::surfMeshSamplePlane::surfMeshSamplePlane @@ -52,21 +203,29 @@ Foam::surfMeshSamplePlane::surfMeshSamplePlane const word& name, const polyMesh& mesh, const plane& planeDesc, - const keyType& zoneKey, + const wordRes& zones, const bool triangulate ) : surfMeshSample(name, mesh), SurfaceSource(planeDesc), - zoneKey_(zoneKey), + zoneNames_(zones), bounds_(), triangulate_(triangulate), needsUpdate_(true) { - if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1) + if (debug) { - Info<< "cellZone(s) " << zoneKey_ - << " not found - using entire mesh" << endl; + if (!zoneNames_.empty()) + { + Info<< "cellZones " << flatOutput(zoneNames_); + + if (-1 == mesh.cellZones().findIndex(zoneNames_)) + { + Info<< " not found!"; + } + Info<< endl; + } } } @@ -80,28 +239,60 @@ Foam::surfMeshSamplePlane::surfMeshSamplePlane : surfMeshSample(name, mesh, dict), SurfaceSource(plane(dict)), - zoneKey_(dict.lookupOrDefault("zone", keyType::null)), + zoneNames_(), bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)), triangulate_(dict.lookupOrDefault("triangulate", true)), needsUpdate_(true) { + if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone")) + { + zoneNames_.resize(1); + dict.readEntry("zone", zoneNames_.first()); + } + + // Make plane relative to the coordinateSystem (Cartesian) // allow lookup from global coordinate systems if (dict.found("coordinateSystem")) { coordinateSystem cs(mesh, dict.subDict("coordinateSystem")); - const point base = cs.globalPosition(planeDesc().origin()); + const point orig = cs.globalPosition(planeDesc().origin()); const vector norm = cs.globalVector(planeDesc().normal()); + if (debug) + { + Info<< "plane " << name << " :" + << " origin:" << origin() + << " normal:" << normal() + << " defined within a local coordinateSystem" << endl; + } + // Assign the plane description - static_cast(*this) = plane(base, norm); + static_cast(*this) = plane(orig, norm); } - if (debug && zoneKey_.size() && mesh.cellZones().findIndex(zoneKey_) == -1) + if (debug) { - Info<< "cellZone(s) " << zoneKey_ - << " not found - using entire mesh" << endl; + Info<< "plane " << name << " :" + << " origin:" << origin() + << " normal:" << normal(); + + if (!bounds_.empty()) + { + Info<< " bounds:" << bounds_; + } + + if (!zoneNames_.empty()) + { + Info<< " cellZones " << flatOutput(zoneNames_); + + if (-1 == mesh.cellZones().findIndex(zoneNames_)) + { + Info<< " not found!"; + } + } + Info<< endl; } } @@ -134,97 +325,7 @@ bool Foam::surfMeshSamplePlane::update() return false; } - 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().selection(zoneKey_).sortedToc() - ); - - 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_); - } - else - { - reCut(mesh(), triangulate_, selectedCells); - } + performCut(mesh(), triangulate_, this->cellSelection(true)); if (debug) { @@ -263,11 +364,11 @@ bool Foam::surfMeshSamplePlane::sample void Foam::surfMeshSamplePlane::print(Ostream& os) const { os << "surfMeshSamplePlane: " << name() << " :" - << " base:" << plane::origin() - << " normal:" << plane::normal() - << " triangulate:" << triangulate_ - << " faces:" << SurfaceSource::surfFaces().size() - << " points:" << SurfaceSource::points().size(); + << " base:" << plane::origin() + << " normal:" << plane::normal() + << " triangulate:" << triangulate_ + << " faces:" << SurfaceSource::surfFaces().size() + << " points:" << SurfaceSource::points().size(); } diff --git a/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.H b/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.H index 923b7b2c37..c1c5959b92 100644 --- a/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.H +++ b/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.H @@ -63,8 +63,8 @@ class surfMeshSamplePlane // Private data - //- If restricted to zones, name of this zone or a regular expression - const keyType zoneKey_; + //- The zone or zones in which cutting is to occur + wordRes zoneNames_; //- Optional bounding box to trim triangles against const boundBox bounds_; @@ -78,6 +78,20 @@ class surfMeshSamplePlane // Private Member Functions + //- Check and warn if bounding box does not intersect mesh or plane + void checkBoundsIntersection + ( + const plane& pln, + const boundBox& meshBb + ) const; + + //- Define cell selection from zones and bounding box. + // Optionally check and warn if the plane does not intersect + // with the bounds of the mesh (or submesh) or if the bounding box + // does not overlap with the mesh (or submesh) + bitSet cellSelection(const bool warnIntersect=false) const; + + //- Sample field on surface template tmp> sampleOnFaces @@ -109,7 +123,7 @@ public: const word& name, const polyMesh& mesh, const plane& planeDesc, - const keyType& zoneKey = word::null, + const wordRes& zones = wordRes(), const bool triangulate = true ); diff --git a/src/sampling/surface/cuttingPlane/cuttingPlane.C b/src/sampling/surface/cuttingPlane/cuttingPlane.C index 580f817aaa..cd8f314482 100644 --- a/src/sampling/surface/cuttingPlane/cuttingPlane.C +++ b/src/sampling/surface/cuttingPlane/cuttingPlane.C @@ -24,117 +24,236 @@ License \*---------------------------------------------------------------------------*/ #include "cuttingPlane.H" -#include "primitiveMesh.H" +#include "fvMesh.H" +#include "volFields.H" #include "linePointRef.H" #include "meshTools.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -// Set values for what is close to zero and what is considered to -// be positive (and not just rounding noise) -//! \cond localScope -static const Foam::scalar zeroish = Foam::SMALL; -static const Foam::scalar positive = Foam::SMALL * 1E3; -//! \endcond +int Foam::cuttingPlane::debug(Foam::debug::debugSwitch("cuttingPlane", 0)); + + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + // Check edge/plane intersection based on crossings ... trivial check. + inline bool intersectsEdge(const PackedList<2>& sides, const edge& e) + { + return (sides[e.first()] != sides[e.last()]); + } + + + // Check for face/plane intersection based on crossings + // Took (-1,0,+1) from plane::sign and packed as (0,1,2). + // Now use for left shift to obtain (1,2,4). + // + // Test accumulated value for an intersection with the plane. + inline bool intersectsFace + ( + const PackedList<2>& sides, + const labelUList& indices + ) + { + unsigned accum = 0u; + + for (const label pointi : indices) + { + accum |= (1u << sides[pointi]); + } + + // Accumulated value 3,5,6,7 indicates an intersection + return (accum == 3 || accum >= 5); + } + +} // End namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::cuttingPlane::calcCutCells +Foam::PackedList<2> Foam::cuttingPlane::classifySides ( - const primitiveMesh& mesh, - const scalarField& dotProducts, - const labelUList& cellIdLabels + const plane& pln, + const pointField& pts ) { - const labelListList& cellEdges = mesh.cellEdges(); - const edgeList& edges = mesh.edges(); + const label len = pts.size(); - const label len = - (notNull(cellIdLabels) ? cellIdLabels.size() : cellEdges.size()); + PackedList<2> output(len); - 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 < len; ++listi) + // From signed (-1,0,+1) to (0,1,2) for PackedList + for (label i=0; i < len; ++i) { - const label celli = - (notNull(cellIdLabels) ? cellIdLabels[listi] : listi); + output.set(i, unsigned(1 + pln.sign(pts[i], SMALL))); + } - const labelList& cEdges = cellEdges[celli]; - label nCutEdges = 0; + return output; +} - for (const label edgei : cEdges) + +Foam::label Foam::cuttingPlane::calcCellCuts +( + const primitiveMesh& mesh, + const PackedList<2>& sides, + bitSet& cellCuts +) +{ + const faceList& faces = mesh.faces(); + + const label nCells = mesh.nCells(); + const label nFaces = mesh.nFaces(); + const label nInternFaces = mesh.nInternalFaces(); + + // Detect cells cuts from the face cuts + + label nFaceCuts = 0; + + // 1st face cut of cell + bitSet hasCut1(nCells); + + // 2nd face cut of cell + bitSet hasCut2(nCells); + + for (label facei = 0; facei < nInternFaces; ++facei) + { + if (intersectsFace(sides, faces[facei])) { - const edge& e = edges[edgei]; + const label own = mesh.faceOwner()[facei]; + const label nei = mesh.faceNeighbour()[facei]; - if - ( - (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive) - || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive) - ) + ++nFaceCuts; + + if (!hasCut1.set(own)) { - if (++nCutEdges > 2) - { - meshCells_[cutcelli++] = celli; - break; - } + hasCut2.set(own); + } + if (!hasCut1.set(nei)) + { + hasCut2.set(nei); } } } - // Set correct list size - meshCells_.resize(cutcelli); + for (label facei = nInternFaces; facei < nFaces; ++facei) + { + if (intersectsFace(sides, faces[facei])) + { + const label own = mesh.faceOwner()[facei]; + + ++nFaceCuts; + + if (!hasCut1.set(own)) + { + hasCut2.set(own); + } + } + } + + hasCut1.clearStorage(); // Not needed now + + if (cellCuts.size()) + { + cellCuts.resize(nCells); // safety + cellCuts &= hasCut2; // restrict to cell subset + + if (debug) + { + Pout<<"detected " << cellCuts.count() << "/" << nCells + << " cells cut, subsetted from " + << hasCut2.count() << "/" << nCells << " cells." << endl; + } + } + else + { + cellCuts = std::move(hasCut2); + + if (debug) + { + Pout<<"detected " << cellCuts.count() << "/" << nCells + << " cells cut." << endl; + } + } + + + if (debug && isA(mesh)) + { + const fvMesh& fvm = dynamicCast(mesh); + + volScalarField debugField + ( + IOobject + ( + "cuttingPlane.cellCuts", + fvm.time().timeName(), + fvm.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + fvm, + dimensionedScalar(dimless, Zero) + ); + + auto& debugFld = debugField.primitiveFieldRef(); + + for (const label celli : cellCuts) + { + debugFld[celli] = 1; + } + + Pout<< "Writing cut types:" + << debugField.objectPath() << endl; + + debugField.write(); + } + + + return nFaceCuts; } void Foam::cuttingPlane::intersectEdges ( const primitiveMesh& mesh, - const scalarField& dotProducts, + const PackedList<2>& sides, + const bitSet& cellCuts, List