Commit ccde68d4 authored by Mark Olesen's avatar Mark Olesen
Browse files

ENH: cellZones support for isoSurface cell/topo sampling variants (#1678)

- better alignment of sampling Cell/Point/Topo inputs

- make exposedPatchName optional for isoSurface, cuttingPlane. This
  was a holdover requirement from an older version of fvMeshSubset
parent 9da52157
......@@ -60,14 +60,21 @@ Foam::sampledDistanceSurface::sampleOnPoints
// Assume volPointInterpolation for the point field!
const auto& volFld = interpolator.psi();
auto tpointFld =
volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
tmp<GeometricField<Type, fvPatchField, volMesh>> tvolFld(volFld);
tmp<GeometricField<Type, pointPatchField, pointMesh>> tpointFld;
return distanceSurface::interpolate
// Interpolated point field
tpointFld.reset
(
(average_ ? pointAverage(tpointFld())() : volFld),
tpointFld()
volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
);
if (average_)
{
tvolFld.reset(pointAverage(tpointFld()));
}
return distanceSurface::interpolate(tvolFld(), tpointFld());
}
......
......@@ -327,16 +327,14 @@ bool Foam::sampledIsoSurface::updateGeometry() const
&& (-1 != mesh().cellZones().findIndex(zoneNames_))
)
{
const polyBoundaryMesh& patches = mesh().boundaryMesh();
// Patch to put exposed internal faces into
const label exposedPatchi = patches.findPatchID(exposedPatchName_);
const label exposedPatchi =
mesh().boundaryMesh().findPatchID(exposedPatchName_);
DebugInfo
<< "Allocating subset of size "
<< mesh().cellZones().selection(zoneNames_).count()
<< " with exposed faces into patch "
<< patches[exposedPatchi].name() << endl;
<< exposedPatchi << endl;
subMeshPtr_.reset
(
......@@ -359,36 +357,25 @@ bool Foam::sampledIsoSurface::updateGeometry() const
// Clear derived data
clearGeom();
refPtr<volScalarField> tvolFld(*volFieldPtr_);
refPtr<pointScalarField> tpointFld(*pointFieldPtr_);
if (subMeshPtr_)
{
const volScalarField& vfld = *volSubFieldPtr_;
isoSurfacePtr_.reset
(
new isoSurfacePoint
(
vfld,
*pointSubFieldPtr_,
isoVal_,
isoParams_
)
);
tvolFld.cref(*volSubFieldPtr_);
tpointFld.cref(*pointSubFieldPtr_);
}
else
{
const volScalarField& vfld = *volFieldPtr_;
isoSurfacePtr_.reset
isoSurfacePtr_.reset
(
new isoSurfacePoint
(
new isoSurfacePoint
(
vfld,
*pointFieldPtr_,
isoVal_,
isoParams_
)
);
}
tvolFld(),
tpointFld(),
isoVal_,
isoParams_
)
);
if (debug)
......@@ -455,21 +442,12 @@ Foam::sampledIsoSurface::sampledIsoSurface
if (-1 != mesh.cellZones().findIndex(zoneNames_))
{
dict.readEntry("exposedPatchName", exposedPatchName_);
if (mesh.boundaryMesh().findPatchID(exposedPatchName_) == -1)
{
FatalIOErrorInFunction(dict)
<< "Cannot find patch " << exposedPatchName_
<< " in which to put exposed faces." << endl
<< "Valid patches are " << mesh.boundaryMesh().names()
<< exit(FatalIOError);
}
dict.readIfPresent("exposedPatchName", exposedPatchName_);
DebugInfo
<< "Restricting to cellZone(s) " << flatOutput(zoneNames_)
<< " with exposed internal faces into patch "
<< exposedPatchName_ << endl;
<< mesh.boundaryMesh().findPatchID(exposedPatchName_) << endl;
}
}
......
......@@ -57,7 +57,7 @@ Usage
bounds | limit with bounding box | no |
zone | limit to cell zone (name or regex) | no |
zones | limit to cell zones (names, regexs) | no |
exposedPatchName | name for zone subset | partly |
exposedPatchName | name for zone subset | optional |
regularise | point snapping (bool or enum) | no | true
mergeTol | tolerance for merging points | no | 1e-6
\endtable
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2018 OpenCFD Ltd.
Copyright (C) 2016-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -27,12 +27,13 @@ License
\*---------------------------------------------------------------------------*/
#include "sampledIsoSurfaceCell.H"
#include "isoSurfaceCell.H"
#include "dictionary.H"
#include "fvMesh.H"
#include "volFields.H"
#include "volPointInterpolation.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
#include "isoSurfaceCell.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -65,8 +66,28 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
// Clear derived data
sampledSurface::clearGeom();
// Use field from database, or try to read it in
// Handle cell zones as inverse (blocked) selection
if (!ignoreCellsPtr_)
{
ignoreCellsPtr_.reset(new bitSet);
if (-1 != mesh().cellZones().findIndex(zoneNames_))
{
bitSet select(mesh().cellZones().selection(zoneNames_));
if (select.any() && !select.all())
{
// From selection to blocking
select.flip();
*ignoreCellsPtr_ = std::move(select);
}
}
}
// Use field from database, or try to read it in
const auto* cellFldPtr = fvm.findObject<volScalarField>(isoField_);
if (debug)
......@@ -111,7 +132,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
// Non-averaged? Use reference
// Field reference (assuming non-averaged)
tmp<scalarField> tcellValues(cellFld.primitiveField());
if (average_)
......@@ -139,25 +160,25 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
}
}
meshedSurface& mySurface = const_cast<sampledIsoSurfaceCell&>(*this);
{
isoSurfaceCell surf
(
fvm,
tcellValues(), // A primitiveField
tpointFld().primitiveField(),
isoVal_,
isoParams_,
*ignoreCellsPtr_
);
isoSurfaceCell surf
(
fvm,
tcellValues(),
tpointFld().primitiveField(),
isoVal_,
isoParams_
);
mySurface.transfer(static_cast<meshedSurface&>(surf));
meshCells_.transfer(surf.meshCells());
mySurface.transfer(static_cast<meshedSurface&>(surf));
meshCells_.transfer(surf.meshCells());
}
if (debug)
{
Pout<< "isoSurfaceCell::updateGeometry() : constructed iso:"
<< nl
Pout<< "isoSurfaceCell::updateGeometry() : constructed iso:" << nl
<< " isoField : " << isoField_ << nl
<< " isoValue : " << isoVal_ << nl
<< " average : " << Switch(average_) << nl
......@@ -188,10 +209,24 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
isoVal_(dict.get<scalar>("isoValue")),
isoParams_(dict),
average_(dict.getOrDefault("average", true)),
zoneNames_(),
prevTimeIndex_(-1),
meshCells_()
meshCells_(),
ignoreCellsPtr_(nullptr)
{
isoParams_.algorithm(isoSurfaceParams::ALGO_CELL); // Force
if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
{
zoneNames_.resize(1);
dict.readEntry("zone", zoneNames_.first());
}
if (-1 != mesh.cellZones().findIndex(zoneNames_))
{
DebugInfo
<< "Restricting to cellZone(s) " << flatOutput(zoneNames_) << endl;
}
}
......@@ -216,6 +251,8 @@ bool Foam::sampledIsoSurfaceCell::expire()
// Clear derived data
sampledSurface::clearGeom();
ignoreCellsPtr_.reset(nullptr);
// Already marked as expired
if (prevTimeIndex_ == -1)
{
......
......@@ -55,13 +55,12 @@ Usage
isoValue | value of iso-surface | yes |
average | cell values from averaged point values | no | false
bounds | limit with bounding box | no |
zone | limit to cell zone (name or regex) | no |
zones | limit to cell zones (names, regexs) | no |
regularise | point snapping | yes |
mergeTol | tolerance for merging points | no | 1e-6
\endtable
Note
Does not currently support cell zones.
SourceFiles
sampledIsoSurfaceCell.C
sampledIsoSurfaceCellTemplates.C
......@@ -105,19 +104,28 @@ class sampledIsoSurfaceCell
//- Parameters (filtering etc) for iso-surface
isoSurfaceParams isoParams_;
//- Whether to recalculate cell values as average of point values
//- Recalculate cell values as average of point values
bool average_;
//- The zone or zones for the iso-surface
wordRes zoneNames_;
// Recreated for every isoSurface
//- Time at last call, also track if surface needs an update
mutable label prevTimeIndex_;
//- For every triangle the original cell in mesh
//- For every face the original cell in mesh (direct storage)
mutable labelList meshCells_;
// Mesh subsetting
//- Cached ignore cells for sub-mesh (zoned)
mutable autoPtr<bitSet> ignoreCellsPtr_;
// Private Member Functions
//- Create iso surface (if time has changed)
......
......@@ -64,30 +64,27 @@ Foam::sampledIsoSurface::sampleOnPoints
// Assume volPointInterpolation for the point field!
const auto& volFld = interpolator.psi();
tmp<GeometricField<Type, fvPatchField, volMesh>> tvolFld(volFld);
tmp<GeometricField<Type, pointPatchField, pointMesh>> tpointFld;
if (subMeshPtr_)
{
auto tvolSubFld = subMeshPtr_->interpolate(volFld);
const auto& volSubFld = tvolSubFld();
auto tpointFld =
volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
return surface().interpolate
(
(average_ ? pointAverage(tpointFld())() : volSubFld),
tpointFld()
);
// Replace with subset
tvolFld.reset(subMeshPtr_->interpolate(volFld));
}
auto tpointFld =
volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
return surface().interpolate
// Interpolated point field
tpointFld.reset
(
(average_ ? pointAverage(tpointFld())() : volFld),
tpointFld()
volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
);
if (average_)
{
tvolFld.reset(pointAverage(tpointFld()));
}
return surface().interpolate(tvolFld(), tpointFld());
}
......
......@@ -65,8 +65,28 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
// Clear derived data
sampledSurface::clearGeom();
// Use field from database, or try to read it in
// Handle cell zones as inverse (blocked) selection
if (!ignoreCellsPtr_)
{
ignoreCellsPtr_.reset(new bitSet);
if (-1 != mesh().cellZones().findIndex(zoneNames_))
{
bitSet select(mesh().cellZones().selection(zoneNames_));
if (select.any() && !select.all())
{
// From selection to blocking
select.flip();
*ignoreCellsPtr_ = std::move(select);
}
}
}
// Use field from database, or try to read it in
const auto* cellFldPtr = fvm.findObject<volScalarField>(isoField_);
if (debug)
......@@ -111,19 +131,50 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
Mesh& mySurface = const_cast<sampledIsoSurfaceTopo&>(*this);
// Field reference (assuming non-averaged)
tmp<scalarField> tcellValues(cellFld.primitiveField());
isoSurfaceTopo surf
(
fvm,
cellFld.primitiveField(),
tpointFld().primitiveField(),
isoVal_,
isoParams_
);
if (average_)
{
// From point field and interpolated cell.
tcellValues = tmp<scalarField>::New(fvm.nCells(), Zero);
auto& cellAvg = tcellValues.ref();
labelField nPointCells(fvm.nCells(), Zero);
for (label pointi = 0; pointi < fvm.nPoints(); ++pointi)
{
const scalar& val = tpointFld().primitiveField()[pointi];
const labelList& pCells = fvm.pointCells(pointi);
for (const label celli : pCells)
{
cellAvg[celli] += val;
++nPointCells[celli];
}
}
forAll(cellAvg, celli)
{
cellAvg[celli] /= nPointCells[celli];
}
}
mySurface.transfer(static_cast<meshedSurface&>(surf));
meshCells_ = std::move(surf.meshCells());
meshedSurface& mySurface = const_cast<sampledIsoSurfaceTopo&>(*this);
{
isoSurfaceTopo surf
(
fvm,
cellFld.primitiveField(),
tpointFld().primitiveField(),
isoVal_,
isoParams_,
*ignoreCellsPtr_
);
mySurface.transfer(static_cast<meshedSurface&>(surf));
meshCells_.transfer(surf.meshCells());
}
// triangulate uses remapFaces()
// - this is somewhat less efficient since it recopies the faces
......@@ -141,6 +192,7 @@ bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
Pout<< "isoSurfaceTopo::updateGeometry() : constructed iso:" << nl
<< " isoField : " << isoField_ << nl
<< " isoValue : " << isoVal_ << nl
<< " average : " << Switch(average_) << nl
<< " filter : "
<< isoSurfaceParams::filterNames[isoParams_.filter()] << nl
<< " triangulate : " << Switch(triangulate_) << nl
......@@ -168,9 +220,12 @@ Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
isoField_(dict.get<word>("isoField")),
isoVal_(dict.get<scalar>("isoValue")),
isoParams_(dict),
average_(dict.getOrDefault("average", false)),
triangulate_(dict.getOrDefault("triangulate", false)),
zoneNames_(),
prevTimeIndex_(-1),
meshCells_()
meshCells_(),
ignoreCellsPtr_(nullptr)
{
isoParams_.algorithm(isoSurfaceParams::ALGO_TOPO); // Force
......@@ -184,6 +239,18 @@ Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
<< "Cannot triangulate without a regularise filter" << nl
<< exit(FatalIOError);
}
if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
{
zoneNames_.resize(1);
dict.readEntry("zone", zoneNames_.first());
}
if (-1 != mesh.cellZones().findIndex(zoneNames_))
{
DebugInfo
<< "Restricting to cellZone(s) " << flatOutput(zoneNames_) << endl;
}
}
......@@ -208,6 +275,8 @@ bool Foam::sampledIsoSurfaceTopo::expire()
// Clear derived data
sampledSurface::clearGeom();
ignoreCellsPtr_.reset(nullptr);
// Already marked as expired
if (prevTimeIndex_ == -1)
{
......
......@@ -53,13 +53,14 @@ Usage
type | isoSurfaceTopo | yes |
isoField | field name for obtaining iso-surface | yes |
isoValue | value of iso-surface | yes |
average | cell values from averaged point values | no | false
bounds | limit with bounding box | no |
zone | limit to cell zone (name or regex) | no |
zones | limit to cell zones (names, regexs) | no |
regularise | filter faces (bool or enum) | no | true
triangulate | triangulate faces (requires regularise) | no | false
\endtable
Note
Does not currently support cell zones.
SourceFiles
sampledIsoSurfaceTopo.C
sampledIsoSurfaceTopoTemplates.C
......@@ -103,19 +104,32 @@ class sampledIsoSurfaceTopo
//- Parameters (filtering etc) for iso-surface
isoSurfaceParams isoParams_;
//- Recalculate cell values as average of point values?
bool average_;
//- Whether to triangulate (after filtering)
bool triangulate_;
//- The zone or zones for the iso-surface
wordRes zoneNames_;
// Recreated for every isoSurface
//- Time at last call, also track it surface needs an update
mutable label prevTimeIndex_;
//- For every triangle/face the original cell in mesh
//- For every face the original cell in mesh (direct storage)
mutable labelList meshCells_;
// Mesh subsetting
//- Cached ignore cells for sub-mesh (zoned)
mutable autoPtr<bitSet> ignoreCellsPtr_;
// Private Member Functions
//- Create iso surface (if time has changed)
......
......@@ -124,18 +124,16 @@ void Foam::sampledCuttingPlane::createGeometry()
&& (-1 != mesh().cellZones().findIndex(zoneNames_))
)
{
const polyBoundaryMesh& patches = mesh().boundaryMesh();
const label exposedPatchi =
mesh().boundaryMesh().findPatchID(exposedPatchName_);
// Patch to put exposed internal faces into
const label exposedPatchi = patches.findPatchID(exposedPatchName_);
bitSet cellsToSelect = mesh().cellZones().selection(zoneNames_);
bitSet cellsToSelect(mesh().cellZones().selection(zoneNames_));
DebugInfo
<< "Allocating subset of size "
<< cellsToSelect.count()
<< " with exposed faces into patch "
<< patches[exposedPatchi].name() << endl;
<< exposedPatchi << endl;
// If we will use a fvMeshSubset so can apply bounds as well to make
......@@ -161,7 +159,10 @@ void Foam::sampledCuttingPlane::createGeometry()
<< cellsToSelect.count() << endl;
}
subMeshPtr_.reset(new fvMeshSubset(fvm, cellsToSelect, exposedPatchi));
subMeshPtr_.reset
(
new fvMeshSubset(fvm, cellsToSelect, exposedPatchi)
);
}
......@@ -209,11 +210,11 @@ void Foam::sampledCuttingPlane::createGeometry()
}
}
volScalarField::Boundary& cellDistanceBf =
cellDistance.boundaryFieldRef();
// Patch fields
{
volScalarField::Boundary& cellDistanceBf =
cellDistance.boundaryFieldRef();
forAll(cellDistanceBf, patchi)
{
if
......@@ -387,21 +388,12 @@ Foam::sampledCuttingPlane::sampledCuttingPlane