zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: unify sampling code for isoSurfaces, support multiple offsets

Browse Source 
- support multiple offsets for cutting plane samples
  and multiple iso-values
This commit is contained in:
Mark Olesen
2020-12-08 02:27:35 +01:00
committed by Andrew Heather
parent 8fe0d1bacd
commit 7a99866db0
13 changed files with 829 additions and 1553 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/sampling/Make/files
View File

@ -43,6 +43,7 @@ sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
sampledSurface/sampledPlane/sampledPlane.C
sampledSurface/isoSurface/sampledIsoSurface.C
sampledSurface/isoSurface/sampledIsoSurfaceCell.C
sampledSurface/isoSurface/sampledIsoSurfacePoint.C
sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
sampledSurface/distanceSurface/sampledDistanceSurface.C
sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C

346
src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
View File

@ -27,12 +27,12 @@ License
\*---------------------------------------------------------------------------*/
#include "sampledIsoSurface.H"
#include "isoSurfacePoint.H"
#include "dictionary.H"
#include "fvMesh.H"
#include "volFields.H"
#include "volPointInterpolation.H"
#include "addToRunTimeSelectionTable.H"
#include "PtrList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -46,13 +46,6 @@ namespace Foam
word,
isoSurface
);
addNamedToRunTimeSelectionTable
(
sampledSurface,
sampledIsoSurface,
word,
isoSurfacePoint
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -312,6 +305,111 @@ void Foam::sampledIsoSurface::getIsoFields() const
}
void Foam::sampledIsoSurface::combineSurfaces
(
PtrList<isoSurfaceBase>& isoSurfPtrs
)
{
isoSurfacePtr_.reset(nullptr);
// Already checked previously for ALGO_POINT, but do it again
// - ALGO_POINT still needs fields (for interpolate)
// The others can do straight transfer
if
(
isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT
&& isoSurfPtrs.size() == 1
)
{
// Shift from list to autoPtr
isoSurfacePtr_.reset(isoSurfPtrs.release(0));
}
else if (isoSurfPtrs.size() == 1)
{
autoPtr<isoSurfaceBase> ptr(isoSurfPtrs.release(0));
auto& surf = *ptr;
surface_.transfer(static_cast<meshedSurface&>(surf));
meshCells_.transfer(surf.meshCells());
}
else
{
// Combine faces with point offsets
//
// Note: use points().size() from surface, not nPoints()
// since there may be uncompacted dangling nodes
label nFaces = 0, nPoints = 0;
for (const auto& surf : isoSurfPtrs)
{
nFaces += surf.size();
nPoints += surf.points().size();
}
faceList newFaces(nFaces);
pointField newPoints(nPoints);
meshCells_.resize(nFaces);
surfZoneList newZones(isoSurfPtrs.size());
nFaces = 0;
nPoints = 0;
forAll(isoSurfPtrs, surfi)
{
autoPtr<isoSurfaceBase> ptr(isoSurfPtrs.release(surfi));
auto& surf = *ptr;
SubList<face> subFaces(newFaces, surf.size(), nFaces);
SubList<point> subPoints(newPoints, surf.points().size(), nPoints);
SubList<label> subCells(meshCells_, surf.size(), nFaces);
newZones[surfi] = surfZone
(
surfZoneIdentifier::defaultName(surfi),
subFaces.size(), // size
nFaces, // start
surfi // index
);
subFaces = surf.surfFaces();
subPoints = surf.points();
subCells = surf.meshCells();
if (nPoints)
{
for (face& f : subFaces)
{
for (label& pointi : f)
{
pointi += nPoints;
}
}
}
nFaces += subFaces.size();
nPoints += subPoints.size();
}
meshedSurface combined
(
std::move(newPoints),
std::move(newFaces),
newZones
);
surface_.transfer(combined);
}
// Addressing into the full mesh
if (subMeshPtr_ && meshCells_.size())
{
meshCells_ =
UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
}
}
bool Foam::sampledIsoSurface::updateGeometry() const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
@ -332,34 +430,76 @@ bool Foam::sampledIsoSurface::updateGeometry() const
// Clear derived data
sampledSurface::clearGeom();
const bool hasCellZones =
(-1 != mesh().cellZones().findIndex(zoneNames_));
// Get sub-mesh if any
// Geometry
if
(
!subMeshPtr_
&& (-1 != mesh().cellZones().findIndex(zoneNames_))
simpleSubMesh_
&& isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
)
{
const label exposedPatchi =
mesh().boundaryMesh().findPatchID(exposedPatchName_);
subMeshPtr_.reset(nullptr);
DebugInfo
<< "Allocating subset of size "
<< mesh().cellZones().selection(zoneNames_).count()
<< " with exposed faces into patch "
<< exposedPatchi << endl;
// Handle cell zones as inverse (blocked) selection
if (!ignoreCellsPtr_)
{
ignoreCellsPtr_.reset(new bitSet);
subMeshPtr_.reset
(
new fvMeshSubset
if (hasCellZones)
{
bitSet select(mesh().cellZones().selection(zoneNames_));
if (select.any() && !select.all())
{
// From selection to blocking
select.flip();
*ignoreCellsPtr_ = std::move(select);
}
}
}
}
else
{
// A standard subMesh treatment
if (ignoreCellsPtr_)
{
ignoreCellsPtr_->clearStorage();
}
else
{
ignoreCellsPtr_.reset(new bitSet);
}
// Get sub-mesh if any
if (!subMeshPtr_ && hasCellZones)
{
const label exposedPatchi =
mesh().boundaryMesh().findPatchID(exposedPatchName_);
DebugInfo
<< "Allocating subset of size "
<< mesh().cellZones().selection(zoneNames_).count()
<< " with exposed faces into patch "
<< exposedPatchi << endl;
subMeshPtr_.reset
(
fvm,
mesh().cellZones().selection(zoneNames_),
exposedPatchi
)
);
new fvMeshSubset
(
fvm,
mesh().cellZones().selection(zoneNames_),
exposedPatchi
)
);
}
}
// The fields
getIsoFields();
refPtr<volScalarField> tvolFld(*volFieldPtr_);
@ -371,23 +511,53 @@ bool Foam::sampledIsoSurface::updateGeometry() const
tpointFld.cref(*pointSubFieldPtr_);
}
isoSurfacePtr_.reset
(
new isoSurfacePoint
(
tvolFld(),
tpointFld(),
isoVal_,
isoParams_
)
);
// Create surfaces for each iso level
PtrList<isoSurfaceBase> isoSurfPtrs(isoValues_.size());
forAll(isoValues_, surfi)
{
isoSurfPtrs.set
(
surfi,
isoSurfaceBase::New
(
isoParams_,
tvolFld(),
tpointFld().primitiveField(),
isoValues_[surfi],
*ignoreCellsPtr_
)
);
}
// And flatten
const_cast<sampledIsoSurface&>(*this)
.combineSurfaces(isoSurfPtrs);
// triangulate uses remapFaces()
// - this is somewhat less efficient since it recopies the faces
// that we just created, but we probably don't want to do this
// too often anyhow.
if
(
triangulate_
&& surface_.size()
&& (isoParams_.algorithm() == isoSurfaceParams::ALGO_TOPO)
)
{
labelList faceMap;
surface_.triangulate(faceMap);
meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
}
if (debug)
{
Pout<< "isoSurfacePoint::updateGeometry() : constructed iso:" << nl
<< " isoField : " << isoField_ << nl
<< " isoValue : " << isoVal_ << nl
Pout<< "isoSurface::updateGeometry() : constructed iso:" << nl
<< " field : " << isoField_ << nl
<< " value : " << flatOutput(isoValues_) << nl
<< " average : " << Switch(average_) << nl
<< " filter : "
<< Switch(bool(isoParams_.filter())) << nl
@ -411,6 +581,7 @@ bool Foam::sampledIsoSurface::updateGeometry() const
Foam::sampledIsoSurface::sampledIsoSurface
(
const isoSurfaceParams& params,
const word& name,
const polyMesh& mesh,
const dictionary& dict
@ -418,33 +589,103 @@ Foam::sampledIsoSurface::sampledIsoSurface
:
sampledSurface(name, mesh, dict),
isoField_(dict.get<word>("isoField")),
isoVal_(dict.get<scalar>("isoValue")),
isoParams_(dict),
isoValues_(),
isoParams_(dict, params),
average_(dict.getOrDefault("average", false)),
triangulate_(dict.getOrDefault("triangulate", false)),
simpleSubMesh_(dict.getOrDefault("simpleSubMesh", false)),
zoneNames_(),
exposedPatchName_(),
prevTimeIndex_(-1),
surface_(),
meshCells_(),
isoSurfacePtr_(nullptr),
subMeshPtr_(nullptr),
ignoreCellsPtr_(nullptr),
storedVolFieldPtr_(nullptr),
volFieldPtr_(nullptr),
pointFieldPtr_(nullptr),
subMeshPtr_(nullptr),
storedVolSubFieldPtr_(nullptr),
volSubFieldPtr_(nullptr),
pointSubFieldPtr_(nullptr)
{
isoParams_.algorithm(isoSurfaceParams::ALGO_POINT); // Force
if (params.algorithm() != isoSurfaceParams::ALGO_DEFAULT)
{
// Forced use of specified algorithm (ignore dictionary entry)
isoParams_.algorithm(params.algorithm());
}
if (!sampledSurface::interpolate())
// The isoValues or isoValue
if (!dict.readIfPresent("isoValues", isoValues_))
{
isoValues_.resize(1);
dict.readEntry("isoValue", isoValues_.first());
}
if (isoValues_.empty())
{
FatalIOErrorInFunction(dict)
<< "Non-interpolated iso surface not supported since triangles"
<< " span across cells." << nl
<< "No isoValue or isoValues specified." << nl
<< exit(FatalIOError);
}
if (isoValues_.size() > 1)
{
const label nOrig = isoValues_.size();
inplaceUniqueSort(isoValues_);
if (nOrig != isoValues_.size())
{
IOWarningInFunction(dict)
<< "Removed non-unique isoValues" << nl;
}
}
if (isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT)
{
// Not possible for ALGO_POINT
simpleSubMesh_ = false;
if (!sampledSurface::interpolate())
{
FatalIOErrorInFunction(dict)
<< "Non-interpolated iso-surface (point) not supported"
<< " since triangles span across cells." << nl
<< exit(FatalIOError);
}
if (isoValues_.size() > 1)
{
FatalIOErrorInFunction(dict)
<< "Multiple values on iso-surface (point) not supported"
<< " since needs original interpolators." << nl
<< exit(FatalIOError);
}
}
if (isoParams_.algorithm() == isoSurfaceParams::ALGO_TOPO)
{
if
(
triangulate_
&& (isoParams_.filter() == isoSurfaceParams::filterType::NONE)
)
{
FatalIOErrorInFunction(dict)
<< "Cannot triangulate without a regularise filter" << nl
<< exit(FatalIOError);
}
}
// Zones
if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
{
zoneNames_.resize(1);
@ -463,6 +704,17 @@ Foam::sampledIsoSurface::sampledIsoSurface
}
Foam::sampledIsoSurface::sampledIsoSurface
(
const word& name,
const polyMesh& mesh,
const dictionary& dict
)
:
sampledIsoSurface(isoSurfaceParams(), name, mesh, dict)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sampledIsoSurface::~sampledIsoSurface()
@ -610,7 +862,7 @@ void Foam::sampledIsoSurface::print(Ostream& os) const
{
os << "isoSurfacePoint: " << name() << " :"
<< " field :" << isoField_
<< " value :" << isoVal_;
<< " value :" << flatOutput(isoValues_);
}

55
src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
View File

@ -28,8 +28,7 @@ Class
Foam::sampledIsoSurface
Description
A sampledSurface defined by a surface of iso value using a
\em point algorithm (always triangulated!).
A sampledSurface defined by a surface of iso value.
It only recalculates the iso-surface if time changes.
To be used in sampleSurfaces / functionObjects.
@ -40,9 +39,10 @@ Usage
{
surface1
{
type isoSurfacePoint;
type isoSurface;
isoField T;
isoValue 373;
isoMethod topo;
}
}
\endverbatim
@ -50,18 +50,26 @@ Usage
Where the sub-entries comprise:
\table
Property | Description | Required | Default
type | isoSurfacePoint / isoSurface | yes |
type | isoSurface | yes |
isoField | field name for obtaining iso-surface | yes |
isoValue | value of iso-surface | yes |
isoValues| values for iso-surfaces | yes |
isoMethod | Iso-algorithm (cell/topo/point) | no | topo
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 |
simpleSubMesh | Simple sub-meshing in algorithm itself | no | false
exposedPatchName | name for zone subset | optional |
regularise | point snapping (bool or enum) | no | true
triangulate | triangulate faces (if regularise) | no | false
mergeTol | tolerance for merging points | no | 1e-6
\endtable
Some options are limited to particular algorithms.
- triangulate is topo-only
- simpleSubMesh and multiple isoValues are not available for point.
SourceFiles
sampledIsoSurface.C
sampledIsoSurfaceTemplates.C
@ -71,9 +79,9 @@ SourceFiles
#ifndef sampledIsoSurface_H
#define sampledIsoSurface_H
#include "isoSurfacePoint.H"
#include "sampledSurface.H"
#include "fvMeshSubset.H"
#include "isoSurfaceBase.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -93,8 +101,8 @@ class sampledIsoSurface
//- Field to get isoSurface of
const word isoField_;
//- Iso value
const scalar isoVal_;
//- The iso-value(s)
List<scalar> isoValues_;
//- Parameters (filtering etc) for iso-surface
isoSurfaceParams isoParams_;
@ -102,6 +110,12 @@ class sampledIsoSurface
//- Whether to recalculate cell values as average of point values
bool average_;
//- Whether to triangulate ALGO_TOPO (after filtering)
bool triangulate_;
//- Use simple sub-meshing in algorithm itself
bool simpleSubMesh_;
//- The zone or zones for the iso-surface
wordRes zoneNames_;
@ -121,7 +135,16 @@ class sampledIsoSurface
mutable labelList meshCells_;
//- Extracted iso-surface, for interpolators
mutable autoPtr<isoSurfacePoint> isoSurfacePtr_;
mutable autoPtr<isoSurfaceBase> isoSurfacePtr_;
// Mesh Subsetting
//- Cached subMesh for (pre-)subset of cell zones
mutable autoPtr<fvMeshSubset> subMeshPtr_;
//- Cached ignore cells for (post-)subset of cell zones
mutable autoPtr<bitSet> ignoreCellsPtr_;
// Fields
@ -133,10 +156,8 @@ class sampledIsoSurface
//- Cached pointfield
mutable const pointScalarField* pointFieldPtr_;
// And on subsetted mesh
//- Cached submesh
mutable autoPtr<fvMeshSubset> subMeshPtr_;
// And on (pre-)subsetted mesh
//- Cached volfield
mutable autoPtr<volScalarField> storedVolSubFieldPtr_;
@ -151,6 +172,9 @@ class sampledIsoSurface
//- Get fields needed to recreate iso surface.
void getIsoFields() const;
//- Collect iso-surfaces into a single surface (No point merging)
void combineSurfaces(PtrList<isoSurfaceBase>& isoSurfPtrs);
//- Create iso surface (if time has changed)
// Do nothing (and return false) if no update was needed
bool updateGeometry() const;
@ -196,6 +220,15 @@ public:
// Constructors
//- Construct from dictionary
sampledIsoSurface
(
const isoSurfaceParams& params,
const word& name,
const polyMesh& mesh,
const dictionary& dict
);
//- Construct from dictionary
sampledIsoSurface
(

351
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
View File

@ -5,8 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2020 OpenCFD Ltd.
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,19 +26,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"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(sampledIsoSurfaceCell, 0);
defineTypeName(sampledIsoSurfaceCell);
addNamedToRunTimeSelectionTable
(
sampledSurface,
@ -49,161 +42,6 @@ namespace Foam
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::sampledIsoSurfaceCell::updateGeometry() const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
// No update needed
if (fvm.time().timeIndex() == prevTimeIndex_)
{
return false;
}
prevTimeIndex_ = fvm.time().timeIndex();
// Clear any previously stored topologies
surface_.clear();
meshCells_.clear();
isoSurfacePtr_.reset(nullptr);
// Clear derived data
sampledSurface::clearGeom();
// 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)
{
if (cellFldPtr)
{
InfoInFunction << "Lookup " << isoField_ << endl;
}
else
{
InfoInFunction
<< "Reading " << isoField_
<< " from time " << fvm.time().timeName()
<< endl;
}
}
// For holding the volScalarField read in.
autoPtr<volScalarField> fieldReadPtr;
if (!cellFldPtr)
{
// Bit of a hack. Read field and store.
fieldReadPtr = autoPtr<volScalarField>::New
(
IOobject
(
isoField_,
fvm.time().timeName(),
fvm,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
fvm
);
}
const volScalarField& cellFld =
(fieldReadPtr ? *fieldReadPtr : *cellFldPtr);
auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
// Field reference (assuming non-averaged)
tmp<scalarField> tcellValues(cellFld.primitiveField());
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];
}
}
{
isoSurfaceCell surf
(
fvm,
tcellValues(), // A primitiveField
tpointFld().primitiveField(),
isoVal_,
isoParams_,
*ignoreCellsPtr_
);
surface_.transfer(static_cast<meshedSurface&>(surf));
meshCells_.transfer(surf.meshCells());
}
// if (subMeshPtr_ && meshCells_.size())
// {
// // With the correct addressing into the full mesh
// meshCells_ =
// UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
// }
if (debug)
{
Pout<< "isoSurfaceCell::updateGeometry() : constructed iso:" << nl
<< " isoField : " << isoField_ << nl
<< " isoValue : " << isoVal_ << nl
<< " average : " << Switch(average_) << nl
<< " filter : "
<< Switch(bool(isoParams_.filter())) << nl
<< " bounds : " << isoParams_.getClipBounds() << nl
<< " points : " << points().size() << nl
<< " faces : " << surface().size() << nl
<< " cut cells : " << meshCells().size() << endl;
}
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -214,185 +52,14 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
const dictionary& dict
)
:
sampledSurface(name, mesh, dict),
isoField_(dict.get<word>("isoField")),
isoVal_(dict.get<scalar>("isoValue")),
isoParams_(dict),
average_(dict.getOrDefault("average", true)),
zoneNames_(),
prevTimeIndex_(-1),
surface_(),
meshCells_(),
isoSurfacePtr_(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;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sampledIsoSurfaceCell::~sampledIsoSurfaceCell()
sampledIsoSurface
(
isoSurfaceParams(isoSurfaceParams::ALGO_CELL),
name,
mesh,
dict
)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::sampledIsoSurfaceCell::needsUpdate() const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
return fvm.time().timeIndex() != prevTimeIndex_;
}
bool Foam::sampledIsoSurfaceCell::expire()
{
surface_.clear();
meshCells_.clear();
isoSurfacePtr_.reset(nullptr);
// Clear derived data
sampledSurface::clearGeom();
ignoreCellsPtr_.reset(nullptr);
// Already marked as expired
if (prevTimeIndex_ == -1)
{
return false;
}
// Force update
prevTimeIndex_ = -1;
return true;
}
bool Foam::sampledIsoSurfaceCell::update()
{
return updateGeometry();
}
Foam::tmp<Foam::scalarField>
Foam::sampledIsoSurfaceCell::sample
(
const interpolation<scalar>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::vectorField>
Foam::sampledIsoSurfaceCell::sample
(
const interpolation<vector>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::sphericalTensorField>
Foam::sampledIsoSurfaceCell::sample
(
const interpolation<sphericalTensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::symmTensorField>
Foam::sampledIsoSurfaceCell::sample
(
const interpolation<symmTensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::tensorField>
Foam::sampledIsoSurfaceCell::sample
(
const interpolation<tensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::scalarField>
Foam::sampledIsoSurfaceCell::interpolate
(
const interpolation<scalar>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::vectorField>
Foam::sampledIsoSurfaceCell::interpolate
(
const interpolation<vector>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::sphericalTensorField>
Foam::sampledIsoSurfaceCell::interpolate
(
const interpolation<sphericalTensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::symmTensorField>
Foam::sampledIsoSurfaceCell::interpolate
(
const interpolation<symmTensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::tensorField>
Foam::sampledIsoSurfaceCell::interpolate
(
const interpolation<tensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
void Foam::sampledIsoSurfaceCell::print(Ostream& os) const
{
os << "isoSurfaceCell: " << name() << " :"
<< " field:" << isoField_
<< " value:" << isoVal_;
//<< " faces:" << faces().size() // possibly no geom yet
//<< " points:" << points().size();
}
// ************************************************************************* //

240
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
View File

@ -63,17 +63,13 @@ Usage
SourceFiles
sampledIsoSurfaceCell.C
sampledIsoSurfaceCellTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef sampledIsoSurfaceCell_H
#define sampledIsoSurfaceCell_H
#include "sampledSurface.H"
#include "MeshedSurface.H"
#include "MeshedSurfacesFwd.H"
#include "isoSurfaceCell.H"
#include "sampledIsoSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -86,94 +82,12 @@ namespace Foam
class sampledIsoSurfaceCell
:
public sampledSurface
public sampledIsoSurface
{
// Private typedefs for convenience
typedef meshedSurface Mesh;
// Private Data
//- Field to get isoSurface of
const word isoField_;
//- Iso value
const scalar isoVal_;
//- Parameters (filtering etc) for iso-surface
isoSurfaceParams isoParams_;
//- Recalculate cell values as average of point values
bool average_;
//- The zone or zones for the iso-surface
wordRes zoneNames_;
// Sampling geometry. Directly stored or via an iso-surface (ALGO_POINT)
//- Time at last call, also track if surface needs an update
mutable label prevTimeIndex_;
//- The extracted surface (direct storage)
mutable meshedSurface surface_;
//- For every face the original cell in mesh (direct storage)
mutable labelList meshCells_;
//- Extracted iso-surface, for interpolators
mutable autoPtr<isoSurfaceCell> isoSurfacePtr_;
// Mesh subsetting
//- Cached ignore cells for sub-mesh (zoned)
mutable autoPtr<bitSet> ignoreCellsPtr_;
// Private Member Functions
//- Create iso surface (if time has changed)
// Do nothing (and return false) if no update was needed
bool updateGeometry() const;
//- Sample volume field onto surface faces
template<class Type>
tmp<Field<Type>> sampleOnFaces
(
const interpolation<Type>& sampler
) const;
//- Interpolate volume field onto surface points
template<class Type>
tmp<Field<Type>> sampleOnPoints
(
const interpolation<Type>& interpolator
) const;
//- Use isoSurfacePtr_ for point interpolation
template<class Type>
tmp<Field<Type>> sampleOnIsoSurfacePoints
(
const interpolation<Type>& interpolator
) const;
protected:
// Protected Member Functions
//- Is currently backed by an isoSurfacePtr_
bool hasIsoSurface() const
{
return bool(isoSurfacePtr_);
}
public:
//- Runtime type information
TypeName("sampledIsoSurfaceCell");
TypeNameNoDebug("sampledIsoSurfaceCell");
// Constructors
@ -188,147 +102,7 @@ public:
//- Destructor
virtual ~sampledIsoSurfaceCell();
// Member Functions
//- Does the surface need an update?
virtual bool needsUpdate() const;
//- Mark the surface as needing an update.
// May also free up unneeded data.
// Return false if surface was already marked as expired.
virtual bool expire();
//- Update the surface as required.
// Do nothing (and return false) if no update was needed
virtual bool update();
//- The currently created surface geometry
const meshedSurface& surface() const
{
if (isoSurfacePtr_)
{
return *isoSurfacePtr_;
}
return surface_;
}
//- For each face, the original cell in mesh
const labelList& meshCells() const
{
if (isoSurfacePtr_)
{
return isoSurfacePtr_->meshCells();
}
return meshCells_;
}
//- Points of surface
virtual const pointField& points() const
{
return surface().points();
}
//- Faces of surface
virtual const faceList& faces() const
{
return surface().surfFaces();
}
//- Per-face zone/region information
virtual const labelList& zoneIds() const
{
return labelList::null();
}
//- Face area magnitudes
virtual const vectorField& Sf() const
{
return surface().Sf();
}
//- Face area magnitudes
virtual const scalarField& magSf() const
{
return surface().magSf();
}
//- Face centres
virtual const vectorField& Cf() const
{
return surface().Cf();
}
// Sample
//- Sample volume field onto surface faces
virtual tmp<scalarField> sample
(
const interpolation<scalar>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<vectorField> sample
(
const interpolation<vector>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<sphericalTensorField> sample
(
const interpolation<sphericalTensor>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<symmTensorField> sample
(
const interpolation<symmTensor>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<tensorField> sample
(
const interpolation<tensor>& sampler
) const;
// Interpolate
//- Interpolate volume field onto surface points
virtual tmp<scalarField> interpolate
(
const interpolation<scalar>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<vectorField> interpolate
(
const interpolation<vector>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<sphericalTensorField> interpolate
(
const interpolation<sphericalTensor>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<symmTensorField> interpolate
(
const interpolation<symmTensor>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<tensorField> interpolate
(
const interpolation<tensor>& interpolator
) const;
//- Write
virtual void print(Ostream& os) const;
virtual ~sampledIsoSurfaceCell() = default;
};
@ -338,12 +112,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "sampledIsoSurfaceCellTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

121
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C
View File

@ -1,120 +1 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2018-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "sampledIsoSurfaceCell.H"
#include "volFieldsFwd.H"
#include "pointFields.H"
#include "volPointInterpolation.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledIsoSurfaceCell::sampleOnFaces
(
const interpolation<Type>& sampler
) const
{
updateGeometry(); // Recreate geometry if time has changed
return sampledSurface::sampleOnFaces
(
sampler,
meshCells(),
surface(),
points()
);
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledIsoSurfaceCell::sampleOnPoints
(
const interpolation<Type>& interpolator
) const
{
updateGeometry(); // Recreate geometry if time has changed
if (isoSurfacePtr_)
{
return this->sampleOnIsoSurfacePoints(interpolator);
}
return sampledSurface::sampleOnPoints
(
interpolator,
meshCells(),
faces(),
points()
);
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledIsoSurfaceCell::sampleOnIsoSurfacePoints
(
const interpolation<Type>& interpolator
) const
{
if (!isoSurfacePtr_)
{
FatalErrorInFunction
<< "cannot call without an iso-surface" << nl
<< exit(FatalError);
}
// 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_)
// {
// // Replace with subset
// tvolFld.reset(subMeshPtr_->interpolate(volFld));
// }
// Interpolated point field
tpointFld.reset
(
volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
);
if (average_)
{
tvolFld.reset(pointAverage(tpointFld()));
}
return isoSurfacePtr_->interpolate(tvolFld(), tpointFld());
}
// ************************************************************************* //
#warning File removed - left for old dependency check only

64
src/sampling/sampledSurface/isoSurface/sampledIsoSurfacePoint.C Normal file
View File

@ -0,0 +1,64 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "sampledIsoSurfacePoint.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeName(sampledIsoSurfacePoint);
addNamedToRunTimeSelectionTable
(
sampledSurface,
sampledIsoSurfacePoint,
word,
isoSurfacePoint
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sampledIsoSurfacePoint::sampledIsoSurfacePoint
(
const word& name,
const polyMesh& mesh,
const dictionary& dict
)
:
sampledIsoSurface
(
isoSurfaceParams(isoSurfaceParams::ALGO_POINT),
name,
mesh,
dict
)
{}
// ************************************************************************* //

116
src/sampling/sampledSurface/isoSurface/sampledIsoSurfacePoint.H Normal file
View File

@ -0,0 +1,116 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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/>.
Class
Foam::sampledIsoSurfacePoint
Description
A sampledSurface defined by a surface of iso value using a
\em point algorithm (always triangulated!).
It only recalculates the iso-surface if time changes.
To be used in sampleSurfaces / functionObjects.
Usage
Example of function object partial specification:
\verbatim
surfaces
{
surface1
{
type isoSurfacePoint;
isoField T;
isoValue 373;
}
}
\endverbatim
Where the sub-entries comprise:
\table
Property | Description | Required | Default
type | isoSurfacePoint | 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 | point snapping | yes |
mergeTol | tolerance for merging points | no | 1e-6
\endtable
SourceFiles
sampledIsoSurfacePoint.C
\*---------------------------------------------------------------------------*/
#ifndef sampledIsoSurfacePoint_H
#define sampledIsoSurfacePoint_H
#include "sampledIsoSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class sampledIsoSurfacePoint Declaration
\*---------------------------------------------------------------------------*/
class sampledIsoSurfacePoint
:
public sampledIsoSurface
{
public:
//- Runtime type information
TypeNameNoDebug("sampledIsoSurfacePoint");
// Constructors
//- Construct from dictionary
sampledIsoSurfacePoint
(
const word& name,
const polyMesh& mesh,
const dictionary& dict
);
//- Destructor
virtual ~sampledIsoSurfacePoint() = default;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

372
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
View File

@ -5,8 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2018 OpenFOAM Foundation
Copyright (C) 2018-2020 OpenCFD Ltd.
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,18 +26,13 @@ License
\*---------------------------------------------------------------------------*/
#include "sampledIsoSurfaceTopo.H"
#include "isoSurfaceTopo.H"
#include "dictionary.H"
#include "fvMesh.H"
#include "volFields.H"
#include "volPointInterpolation.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(sampledIsoSurfaceTopo, 0);
defineTypeName(sampledIsoSurfaceTopo);
addNamedToRunTimeSelectionTable
(
sampledSurface,
@ -48,175 +42,6 @@ namespace Foam
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::sampledIsoSurfaceTopo::updateGeometry() const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
// No update needed
if (fvm.time().timeIndex() == prevTimeIndex_)
{
return false;
}
prevTimeIndex_ = fvm.time().timeIndex();
// Clear any previously stored topologies
surface_.clear();
meshCells_.clear();
isoSurfacePtr_.reset(nullptr);
// Clear derived data
sampledSurface::clearGeom();
// 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)
{
if (cellFldPtr)
{
InfoInFunction << "Lookup " << isoField_ << endl;
}
else
{
InfoInFunction
<< "Reading " << isoField_
<< " from time " << fvm.time().timeName()
<< endl;
}
}
// For holding the volScalarField read in.
autoPtr<volScalarField> fieldReadPtr;
if (!cellFldPtr)
{
// Bit of a hack. Read field and store.
fieldReadPtr = autoPtr<volScalarField>::New
(
IOobject
(
isoField_,
fvm.time().timeName(),
fvm,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
fvm
);
}
const volScalarField& cellFld =
(fieldReadPtr ? *fieldReadPtr : *cellFldPtr);
auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
// Field reference (assuming non-averaged)
tmp<scalarField> tcellValues(cellFld.primitiveField());
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];
}
}
{
isoSurfaceTopo surf
(
fvm,
cellFld.primitiveField(),
tpointFld().primitiveField(),
isoVal_,
isoParams_,
*ignoreCellsPtr_
);
surface_.transfer(static_cast<meshedSurface&>(surf));
meshCells_.transfer(surf.meshCells());
}
// if (subMeshPtr_ && meshCells_.size())
// {
// // With the correct addressing into the full mesh
// meshCells_ =
// UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
// }
// triangulate uses remapFaces()
// - this is somewhat less efficient since it recopies the faces
// that we just created, but we probably don't want to do this
// too often anyhow.
if (triangulate_ && surface_.size())
{
labelList faceMap;
surface_.triangulate(faceMap);
meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
}
if (debug)
{
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
<< " bounds : " << isoParams_.getClipBounds() << nl
<< " points : " << points().size() << nl
<< " faces : " << surface().size() << nl
<< " cut cells : " << meshCells().size() << endl;
}
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
@ -226,197 +51,14 @@ Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
const dictionary& dict
)
:
sampledSurface(name, mesh, dict),
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),
surface_(),
meshCells_(),
isoSurfacePtr_(nullptr)
{
isoParams_.algorithm(isoSurfaceParams::ALGO_TOPO); // Force
if
sampledIsoSurface
(
triangulate_
&& (isoParams_.filter() == isoSurfaceParams::filterType::NONE)
isoSurfaceParams(isoSurfaceParams::ALGO_TOPO),
name,
mesh,
dict
)
{
FatalIOErrorInFunction(dict)
<< "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;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sampledIsoSurfaceTopo::~sampledIsoSurfaceTopo()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::sampledIsoSurfaceTopo::needsUpdate() const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
return fvm.time().timeIndex() != prevTimeIndex_;
}
bool Foam::sampledIsoSurfaceTopo::expire()
{
surface_.clear();
meshCells_.clear();
isoSurfacePtr_.reset(nullptr);
// Clear derived data
sampledSurface::clearGeom();
ignoreCellsPtr_.reset(nullptr);
// Already marked as expired
if (prevTimeIndex_ == -1)
{
return false;
}
// Force update
prevTimeIndex_ = -1;
return true;
}
bool Foam::sampledIsoSurfaceTopo::update()
{
return updateGeometry();
}
Foam::tmp<Foam::scalarField>
Foam::sampledIsoSurfaceTopo::sample
(
const interpolation<scalar>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::vectorField>
Foam::sampledIsoSurfaceTopo::sample
(
const interpolation<vector>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::sphericalTensorField>
Foam::sampledIsoSurfaceTopo::sample
(
const interpolation<sphericalTensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::symmTensorField>
Foam::sampledIsoSurfaceTopo::sample
(
const interpolation<symmTensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::tensorField>
Foam::sampledIsoSurfaceTopo::sample
(
const interpolation<tensor>& sampler
) const
{
return sampleOnFaces(sampler);
}
Foam::tmp<Foam::scalarField>
Foam::sampledIsoSurfaceTopo::interpolate
(
const interpolation<scalar>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::vectorField>
Foam::sampledIsoSurfaceTopo::interpolate
(
const interpolation<vector>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::sphericalTensorField>
Foam::sampledIsoSurfaceTopo::interpolate
(
const interpolation<sphericalTensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::symmTensorField>
Foam::sampledIsoSurfaceTopo::interpolate
(
const interpolation<symmTensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
Foam::tmp<Foam::tensorField>
Foam::sampledIsoSurfaceTopo::interpolate
(
const interpolation<tensor>& interpolator
) const
{
return sampleOnPoints(interpolator);
}
void Foam::sampledIsoSurfaceTopo::print(Ostream& os) const
{
os << "isoSurfaceTopo: " << name() << " :"
<< " field:" << isoField_
<< " value:" << isoVal_;
//<< " faces:" << faces().size() // possibly no geom yet
//<< " points:" << points().size();
}
// ************************************************************************* //

243
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
View File

@ -63,17 +63,13 @@ Usage
SourceFiles
sampledIsoSurfaceTopo.C
sampledIsoSurfaceTopoTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef sampledIsoSurfaceTopo_H
#define sampledIsoSurfaceTopo_H
#include "sampledSurface.H"
#include "MeshedSurface.H"
#include "MeshedSurfacesFwd.H"
#include "isoSurfaceTopo.H"
#include "sampledIsoSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -86,97 +82,12 @@ namespace Foam
class sampledIsoSurfaceTopo
:
public sampledSurface
public sampledIsoSurface
{
// Private typedefs for convenience
typedef meshedSurface Mesh;
// Private Data
//- Field to get isoSurface of
const word isoField_;
//- Iso value
const scalar isoVal_;
//- 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_;
// Sampling geometry. Directly stored or via an iso-surface (ALGO_POINT)
//- Time at last call, also track it surface needs an update
mutable label prevTimeIndex_;
//- The extracted surface (direct storage)
mutable meshedSurface surface_;
//- For every face the original cell in mesh (direct storage)
mutable labelList meshCells_;
//- Extracted iso-surface, for interpolators
mutable autoPtr<isoSurfaceTopo> isoSurfacePtr_;
// Mesh subsetting
//- Cached ignore cells for sub-mesh (zoned)
mutable autoPtr<bitSet> ignoreCellsPtr_;
// Private Member Functions
//- Create iso surface (if time has changed)
// Do nothing (and return false) if no update was needed
bool updateGeometry() const;
//- Sample volume field onto surface faces
template<class Type>
tmp<Field<Type>> sampleOnFaces
(
const interpolation<Type>& sampler
) const;
//- Interpolate volume field onto surface points
template<class Type>
tmp<Field<Type>> sampleOnPoints
(
const interpolation<Type>& interpolator
) const;
//- Use isoSurfacePtr_ for point interpolation
template<class Type>
tmp<Field<Type>> sampleOnIsoSurfacePoints
(
const interpolation<Type>& interpolator
) const;
protected:
// Protected Member Functions
//- Is currently backed by an isoSurfacePtr_
bool hasIsoSurface() const
{
return bool(isoSurfacePtr_);
}
public:
//- Runtime type information
TypeName("sampledIsoSurfaceTopo");
TypeNameNoDebug("sampledIsoSurfaceTopo");
// Constructors
@ -191,147 +102,7 @@ public:
//- Destructor
virtual ~sampledIsoSurfaceTopo();
// Member Functions
//- Does the surface need an update?
virtual bool needsUpdate() const;
//- Mark the surface as needing an update.
// May also free up unneeded data.
// Return false if surface was already marked as expired.
virtual bool expire();
//- Update the surface as required.
// Do nothing (and return false) if no update was needed
virtual bool update();
//- The currently created surface geometry
const meshedSurface& surface() const
{
if (isoSurfacePtr_)
{
return *isoSurfacePtr_;
}
return surface_;
}
//- For each face, the original cell in mesh
const labelList& meshCells() const
{
if (isoSurfacePtr_)
{
return isoSurfacePtr_->meshCells();
}
return meshCells_;
}
//- Points of surface
virtual const pointField& points() const
{
return surface().points();
}
//- Faces of surface
virtual const faceList& faces() const
{
return surface().surfFaces();
}
//- Per-face zone/region information
virtual const labelList& zoneIds() const
{
return labelList::null();
}
//- Face area magnitudes
virtual const vectorField& Sf() const
{
return surface().Sf();
}
//- Face area magnitudes
virtual const scalarField& magSf() const
{
return surface().magSf();
}
//- Face centres
virtual const vectorField& Cf() const
{
return surface().Cf();
}
// Sample
//- Sample volume field onto surface faces
virtual tmp<scalarField> sample
(
const interpolation<scalar>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<vectorField> sample
(
const interpolation<vector>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<sphericalTensorField> sample
(
const interpolation<sphericalTensor>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<symmTensorField> sample
(
const interpolation<symmTensor>& sampler
) const;
//- Sample volume field onto surface faces
virtual tmp<tensorField> sample
(
const interpolation<tensor>& sampler
) const;
// Interpolate
//- Interpolate volume field onto surface points
virtual tmp<scalarField> interpolate
(
const interpolation<scalar>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<vectorField> interpolate
(
const interpolation<vector>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<sphericalTensorField> interpolate
(
const interpolation<sphericalTensor>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<symmTensorField> interpolate
(
const interpolation<symmTensor>& interpolator
) const;
//- Interpolate volume field onto surface points
virtual tmp<tensorField> interpolate
(
const interpolation<tensor>& interpolator
) const;
//- Write
virtual void print(Ostream& os) const;
virtual ~sampledIsoSurfaceTopo() = default;
};
@ -341,12 +112,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "sampledIsoSurfaceTopoTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

121
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
View File

@ -1,120 +1 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2018 OpenFOAM Foundation
Copyright (C) 2018-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 "sampledIsoSurfaceTopo.H"
#include "volFieldsFwd.H"
#include "pointFields.H"
#include "volPointInterpolation.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledIsoSurfaceTopo::sampleOnFaces
(
const interpolation<Type>& sampler
) const
{
updateGeometry(); // Recreate geometry if time has changed
return sampledSurface::sampleOnFaces
(
sampler,
meshCells(),
faces(),
points()
);
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledIsoSurfaceTopo::sampleOnPoints
(
const interpolation<Type>& interpolator
) const
{
updateGeometry(); // Recreate geometry if time has changed
if (isoSurfacePtr_)
{
return this->sampleOnIsoSurfacePoints(interpolator);
}
return sampledSurface::sampleOnPoints
(
interpolator,
meshCells(),
faces(),
points()
);
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledIsoSurfaceTopo::sampleOnIsoSurfacePoints
(
const interpolation<Type>& interpolator
) const
{
if (!isoSurfacePtr_)
{
FatalErrorInFunction
<< "cannot call without an iso-surface" << nl
<< exit(FatalError);
}
// 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_)
// {
// // Replace with subset
// tvolFld.reset(subMeshPtr_->interpolate(volFld));
// }
// Interpolated point field
tpointFld.reset
(
volPointInterpolation::New(tvolFld().mesh()).interpolate(tvolFld())
);
if (average_)
{
tvolFld.reset(pointAverage(tpointFld()));
}
return isoSurfacePtr_->interpolate(tvolFld(), tpointFld());
}
// ************************************************************************* //
#warning File removed - left for old dependency check only

314
src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
View File

@ -28,13 +28,11 @@ License
#include "sampledCuttingPlane.H"
#include "dictionary.H"
#include "fvMesh.H"
#include "volFields.H"
#include "volPointInterpolation.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
#include "isoSurfaceCell.H"
#include "isoSurfacePoint.H"
#include "isoSurfaceTopo.H"
#include "PtrList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -189,6 +187,111 @@ void Foam::sampledCuttingPlane::setDistanceFields(const plane& pln)
}
void Foam::sampledCuttingPlane::combineSurfaces
(
PtrList<isoSurfaceBase>& isoSurfPtrs
)
{
isoSurfacePtr_.reset(nullptr);
// Already checked previously for ALGO_POINT, but do it again
// - ALGO_POINT still needs fields (for interpolate)
// The others can do straight transfer
if
(
isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT
&& isoSurfPtrs.size() == 1
)
{
// Shift ownership from list to autoPtr
isoSurfacePtr_.reset(isoSurfPtrs.release(0));
}
else if (isoSurfPtrs.size() == 1)
{
autoPtr<isoSurfaceBase> ptr(isoSurfPtrs.release(0));
auto& surf = *ptr;
surface_.transfer(static_cast<meshedSurface&>(surf));
meshCells_.transfer(surf.meshCells());
}
else
{
// Combine faces with point offsets
//
// Note: use points().size() from surface, not nPoints()
// since there may be uncompacted dangling nodes
label nFaces = 0, nPoints = 0;
for (const auto& surf : isoSurfPtrs)
{
nFaces += surf.size();
nPoints += surf.points().size();
}
faceList newFaces(nFaces);
pointField newPoints(nPoints);
meshCells_.resize(nFaces);
surfZoneList newZones(isoSurfPtrs.size());
nFaces = 0;
nPoints = 0;
forAll(isoSurfPtrs, surfi)
{
autoPtr<isoSurfaceBase> ptr(isoSurfPtrs.release(surfi));
auto& surf = *ptr;
SubList<face> subFaces(newFaces, surf.size(), nFaces);
SubList<point> subPoints(newPoints, surf.points().size(), nPoints);
SubList<label> subCells(meshCells_, surf.size(), nFaces);
newZones[surfi] = surfZone
(
surfZoneIdentifier::defaultName(surfi),
subFaces.size(), // size
nFaces, // start
surfi // index
);
subFaces = surf.surfFaces();
subPoints = surf.points();
subCells = surf.meshCells();
if (nPoints)
{
for (face& f : subFaces)
{
for (label& pointi : f)
{
pointi += nPoints;
}
}
}
nFaces += subFaces.size();
nPoints += subPoints.size();
}
meshedSurface combined
(
std::move(newPoints),
std::move(newFaces),
newZones
);
surface_.transfer(combined);
}
// Addressing into the full mesh
if (subMeshPtr_ && meshCells_.size())
{
meshCells_ =
UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
}
}
void Foam::sampledCuttingPlane::createGeometry()
{
if (debug)
@ -209,54 +312,94 @@ void Foam::sampledCuttingPlane::createGeometry()
pointDistance_.clear();
cellDistancePtr_.clear();
const bool hasCellZones =
(-1 != mesh().cellZones().findIndex(zoneNames_));
const fvMesh& fvm = static_cast<const fvMesh&>(this->mesh());
// Get sub-mesh if any
// Geometry
if
(
!subMeshPtr_
&& (-1 != mesh().cellZones().findIndex(zoneNames_))
simpleSubMesh_
&& isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
)
{
const label exposedPatchi =
mesh().boundaryMesh().findPatchID(exposedPatchName_);
subMeshPtr_.reset(nullptr);
bitSet cellsToSelect(mesh().cellZones().selection(zoneNames_));
DebugInfo
<< "Allocating subset of size "
<< cellsToSelect.count()
<< " with exposed faces into patch "
<< exposedPatchi << endl;
// If we will use a fvMeshSubset so can apply bounds as well to make
// the initial selection smaller.
const boundBox& clipBb = isoParams_.getClipBounds();
if (clipBb.valid() && cellsToSelect.any())
// Handle cell zones as inverse (blocked) selection
if (!ignoreCellsPtr_)
{
const auto& cellCentres = fvm.C();
ignoreCellsPtr_.reset(new bitSet);
for (const label celli : cellsToSelect)
if (hasCellZones)
{
const point& cc = cellCentres[celli];
bitSet select(mesh().cellZones().selection(zoneNames_));
if (!clipBb.contains(cc))
if (select.any() && !select.all())
{
cellsToSelect.unset(celli);
// From selection to blocking
select.flip();
*ignoreCellsPtr_ = std::move(select);
}
}
}
}
else
{
// A standard subMesh treatment
DebugInfo
<< "Bounded subset of size "
<< cellsToSelect.count() << endl;
if (ignoreCellsPtr_)
{
ignoreCellsPtr_->clearStorage();
}
else
{
ignoreCellsPtr_.reset(new bitSet);
}
subMeshPtr_.reset
(
new fvMeshSubset(fvm, cellsToSelect, exposedPatchi)
);
// Get sub-mesh if any
if (!subMeshPtr_ && hasCellZones)
{
const label exposedPatchi =
mesh().boundaryMesh().findPatchID(exposedPatchName_);
bitSet cellsToSelect(mesh().cellZones().selection(zoneNames_));
DebugInfo
<< "Allocating subset of size "
<< cellsToSelect.count() << " with exposed faces into patch "
<< exposedPatchi << endl;
// If we will use a fvMeshSubset so can apply bounds as well to make
// the initial selection smaller.
const boundBox& clipBb = isoParams_.getClipBounds();
if (clipBb.valid() && cellsToSelect.any())
{
const auto& cellCentres = fvm.C();
for (const label celli : cellsToSelect)
{
const point& cc = cellCentres[celli];
if (!clipBb.contains(cc))
{
cellsToSelect.unset(celli);
}
}
DebugInfo
<< "Bounded subset of size "
<< cellsToSelect.count() << endl;
}
subMeshPtr_.reset
(
new fvMeshSubset(fvm, cellsToSelect, exposedPatchi)
);
}
}
@ -291,7 +434,6 @@ void Foam::sampledCuttingPlane::createGeometry()
dimensionedScalar(dimLength, Zero)
)
);
const volScalarField& cellDistance = cellDistancePtr_();
setDistanceFields(plane_);
@ -321,35 +463,35 @@ void Foam::sampledCuttingPlane::createGeometry()
}
isoSurfacePtr_.reset
(
isoSurfaceBase::New
(
isoParams_,
cellDistance,
pointDistance_,
scalar(0)
// nothing ignored: ignoreCells
)
);
// Create surfaces for each offset
// ALGO_POINT uses cell field and point field
// The others can do straight transfer
if (isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT)
PtrList<isoSurfaceBase> isoSurfPtrs(offsets_.size());
forAll(offsets_, surfi)
{
surface_.transfer(static_cast<meshedSurface&>(*isoSurfacePtr_));
meshCells_.transfer(isoSurfacePtr_->meshCells());
isoSurfacePtr_.reset(nullptr);
cellDistancePtr_.reset(nullptr);
pointDistance_.clear();
isoSurfPtrs.set
(
surfi,
isoSurfaceBase::New
(
isoParams_,
cellDistance,
pointDistance_,
offsets_[surfi],
*ignoreCellsPtr_
)
);
}
if (subMeshPtr_ && meshCells_.size())
// And flatten
combineSurfaces(isoSurfPtrs);
// Discard fields if not required by an iso-surface
if (!isoSurfacePtr_)
{
// With the correct addressing into the full mesh
meshCells_ =
UIndirectList<label>(subMeshPtr_->cellMap(), meshCells_);
cellDistancePtr_.reset(nullptr);
pointDistance_.clear();
}
if (debug)
@ -371,6 +513,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
:
sampledSurface(name, mesh, dict),
plane_(dict),
offsets_(),
isoParams_
(
dict,
@ -378,15 +521,59 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
isoSurfaceParams::filterType::DIAGCELL
),
average_(dict.getOrDefault("average", false)),
simpleSubMesh_(dict.getOrDefault("simpleSubMesh", false)),
zoneNames_(),
exposedPatchName_(),
needsUpdate_(true),
subMeshPtr_(nullptr),
cellDistancePtr_(nullptr),
surface_(),
meshCells_(),
isoSurfacePtr_(nullptr)
isoSurfacePtr_(nullptr),
subMeshPtr_(nullptr),
ignoreCellsPtr_(nullptr),
cellDistancePtr_(nullptr),
pointDistance_()
{
dict.readIfPresent("offsets", offsets_);
if (offsets_.empty())
{
offsets_.resize(1);
offsets_.first() = Zero;
}
if (offsets_.size() > 1)
{
const label nOrig = offsets_.size();
inplaceUniqueSort(offsets_);
if (nOrig != offsets_.size())
{
IOWarningInFunction(dict)
<< "Removed non-unique offsets" << nl;
}
}
if (isoParams_.algorithm() == isoSurfaceParams::ALGO_POINT)
{
// Not possible for ALGO_POINT
simpleSubMesh_ = false;
// Not possible for ALGO_POINT
if (offsets_.size() > 1)
{
FatalIOErrorInFunction(dict)
<< "Multiple offsets with iso-surface (point) not supported"
<< " since needs original interpolators." << nl
<< exit(FatalIOError);
}
}
// Zones
if (!dict.readIfPresent("zones", zoneNames_) && dict.found("zone"))
{
zoneNames_.resize(1);
@ -563,6 +750,7 @@ void Foam::sampledCuttingPlane::print(Ostream& os) const
{
os << "sampledCuttingPlane: " << name() << " :"
<< " plane:" << plane_
<< " offsets:" << flatOutput(offsets_)
<< " faces:" << faces().size()
<< " points:" << points().size();
}

38
src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
View File

@ -55,6 +55,7 @@ Usage
Property | Description | Required | Default
type | cuttingPlane | yes |
planeType | plane description (pointAndNormal etc) | yes |
offsets | Offsets of the origin in the normal direction | no | (0)
isoMethod | Iso-algorithm (cell/topo/point) | no | topo
bounds | limit with bounding box | no |
zone | limit to cell zone (name or regex) | no |
@ -101,12 +102,18 @@ class sampledCuttingPlane
//- Plane
const plane plane_;
//- The offsets to the plane - defaults to (0).
List<scalar> offsets_;
//- Parameters (filtering etc) for iso-surface
isoSurfaceParams isoParams_;
//- Whether to recalculate cell values as average of point values
bool average_;
//- Use simple sub-meshing in algorithm itself
bool simpleSubMesh_;
//- The zone or zones in which cutting is to occur
wordRes zoneNames_;
@ -117,16 +124,6 @@ class sampledCuttingPlane
mutable bool needsUpdate_;
//- Mesh subset (optional: only used with zones)
autoPtr<fvMeshSubset> subMeshPtr_;
//- Distance to cell centres
autoPtr<volScalarField> cellDistancePtr_;
//- Distance to points
scalarField pointDistance_;
// Sampling geometry. Directly stored or via an iso-surface (ALGO_POINT)
//- The extracted surface (direct storage)
@ -139,6 +136,24 @@ class sampledCuttingPlane
autoPtr<isoSurfaceBase> isoSurfacePtr_;
// Mesh Subsetting
//- Cached subMesh for (pre-)subset of cell zones
mutable autoPtr<fvMeshSubset> subMeshPtr_;
//- Cached ignore cells for (post-)subset of cell zones
mutable autoPtr<bitSet> ignoreCellsPtr_;
// Fields
//- Distance to cell centres
autoPtr<volScalarField> cellDistancePtr_;
//- Distance to points
scalarField pointDistance_;
// Private Member Functions
//- Check and warn if bounding box does not intersect mesh or plane
@ -151,6 +166,9 @@ class sampledCuttingPlane
//- Fill cellDistance, pointDistance fields for the specified plane
void setDistanceFields(const plane& pln);
//- Collect iso-surfaces into a single surface (No point merging)
void combineSurfaces(PtrList<isoSurfaceBase>& isoSurfPtrs);
//- Create iso surface
void createGeometry();