From 7943603bf0a8f04794791f84d73e3138f7aa77ce Mon Sep 17 00:00:00 2001
From: Mattijs Janssens <m.janssens@opencfd.co.uk>
Date: Thu, 6 Dec 2018 15:56:26 +0000
Subject: [PATCH] ENH: isoSurfaceTopo: replacement for isoSurfaceCell.
---
src/sampling/Make/files | 2 +
.../isoSurface/sampledIsoSurfaceTopo.C | 328 +++++
.../isoSurface/sampledIsoSurfaceTopo.H | 291 ++++
.../sampledIsoSurfaceTopoTemplates.C | 95 ++
.../surface/isoSurface/isoSurfaceTopo.C | 1234 +++++++++++++++++
.../surface/isoSurface/isoSurfaceTopo.H | 266 ++++
.../isoSurface/isoSurfaceTopoTemplates.C | 91 ++
7 files changed, 2307 insertions(+)
create mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
create mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
create mode 100644 src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
create mode 100644 src/sampling/surface/isoSurface/isoSurfaceTopo.C
create mode 100644 src/sampling/surface/isoSurface/isoSurfaceTopo.H
create mode 100644 src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index 3e42081cc06..dbef3b39873 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -29,6 +29,7 @@ surface/cutting/cuttingSurfaceBaseSelection.C
surface/distanceSurface/distanceSurface.C
surface/isoSurface/isoSurface.C
surface/isoSurface/isoSurfaceCell.C
+surface/isoSurface/isoSurfaceTopo.C
surface/thresholdCellFaces/thresholdCellFaces.C
surface/triSurfaceMesh/discreteSurface.C
@@ -44,6 +45,7 @@ sampledSurface/sampledPatchInternalField/sampledPatchInternalField.C
sampledSurface/sampledPlane/sampledPlane.C
sampledSurface/isoSurface/sampledIsoSurface.C
sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
sampledSurface/distanceSurface/sampledDistanceSurface.C
sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
sampledSurface/sampledCuttingSurface/sampledCuttingSurface.C
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
new file mode 100644
index 00000000000..6f3ff764f0c
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.C
@@ -0,0 +1,328 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2018 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 "dictionary.H"
+#include "volFields.H"
+#include "volPointInterpolation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvMesh.H"
+#include "isoSurfaceTopo.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(sampledIsoSurfaceTopo, 0);
+ addNamedToRunTimeSelectionTable
+ (
+ sampledSurface,
+ sampledIsoSurfaceTopo,
+ word,
+ isoSurfaceTopo
+ );
+}
+
+// * * * * * * * * * * * * * 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 derived data
+ sampledSurface::clearGeom();
+
+ // 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.valid() ? *fieldReadPtr : *cellFldPtr);
+
+ auto tpointFld = volPointInterpolation::New(fvm).interpolate(cellFld);
+
+ //- Direct from cell field and point field. Gives bad continuity.
+ isoSurfaceTopo surf
+ (
+ fvm,
+ cellFld.primitiveField(),
+ tpointFld().primitiveField(),
+ isoVal_,
+ (regularise_ ? isoSurfaceTopo::DIAGCELL : isoSurfaceTopo::NONE)
+ );
+
+ MeshedSurface<face>& mySurface = const_cast<sampledIsoSurfaceTopo&>(*this);
+
+ mySurface.transfer(static_cast<meshedSurface&>(surf));
+ meshCells_ = std::move(surf.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_)
+ {
+ labelList faceMap;
+ mySurface.triangulate(faceMap);
+ meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+ }
+
+ if (debug)
+ {
+ Pout<< "sampledIsoSurfaceTopo::updateGeometry() : constructed iso:"
+ << nl
+ << " regularise : " << regularise_ << nl
+ << " triangulate : " << triangulate_ << nl
+ << " isoField : " << isoField_ << nl
+ << " isoValue : " << isoVal_ << nl
+ << " points : " << points().size() << nl
+ << " faces : " << MeshStorage::size() << nl
+ << " cut cells : " << meshCells_.size() << endl;
+ }
+
+ return true;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::sampledIsoSurfaceTopo::sampledIsoSurfaceTopo
+(
+ const word& name,
+ const polyMesh& mesh,
+ const dictionary& dict
+)
+:
+ sampledSurface(name, mesh, dict),
+ MeshStorage(),
+ isoField_(dict.get<word>("isoField")),
+ isoVal_(dict.get<scalar>("isoValue")),
+ regularise_(dict.lookupOrDefault("regularise", true)),
+ triangulate_(dict.lookupOrDefault("triangulate", false)),
+ prevTimeIndex_(-1),
+ meshCells_()
+{
+ if (triangulate_ && !regularise_)
+ {
+ FatalIOErrorInFunction(dict) << "Cannot both use regularise"
+ << " and triangulate" << exit(FatalIOError);
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * 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()
+{
+ // Clear derived data
+ sampledSurface::clearGeom();
+
+ // 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 << "sampledIsoSurfaceTopo: " << name() << " :"
+ << " field:" << isoField_
+ << " value:" << isoVal_;
+ //<< " faces:" << faces().size() // possibly no geom yet
+ //<< " points:" << points().size();
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
new file mode 100644
index 00000000000..31612b92311
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopo.H
@@ -0,0 +1,291 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2018 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::sampledIsoSurfaceTopo
+
+Description
+ A sampledSurface defined by a surface of iso value.
+ To be used in sampleSurfaces / functionObjects. Recalculates iso surface
+ only if time changes.
+
+ This is often embedded as part of a sampled surfaces function object.
+
+Usage
+ Example of function object partial specification:
+ \verbatim
+ surfaces
+ (
+ surface1
+ {
+ type isoSurfaceTopo;
+ isoField p;
+ isoValue 0.0;
+ }
+ );
+ \endverbatim
+
+ Where the sub-entries comprise:
+ \table
+ Property | Description | Required | Default
+ type | isoSurfaceTopo | yes |
+ isoField | field name for obtaining iso-surface | yes |
+ isoValue | value of iso-surface | yes |
+ regularise | filter faces | no | true
+ triangulate | triangulate faces (if regularise) | no | false
+ \endtable
+
+Note
+ Does not currently support cell zones.
+
+SourceFiles
+ sampledIsoSurfaceTopo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledIsoSurfaceTopo_H
+#define sampledIsoSurfaceTopo_H
+
+#include "sampledSurface.H"
+#include "MeshedSurface.H"
+#include "MeshedSurfacesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class sampledIsoSurfaceTopo Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledIsoSurfaceTopo
+:
+ public sampledSurface,
+ public MeshedSurface<face>
+{
+ // Private typedefs for convenience
+ typedef MeshedSurface<face> MeshStorage;
+
+ // Private data
+
+ //- Field to get isoSurface of
+ const word isoField_;
+
+ //- Iso value
+ const scalar isoVal_;
+
+ //- Whether to coarse
+ const bool regularise_;
+
+ //- Whether to triangulate
+ const bool triangulate_;
+
+ // 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
+ mutable labelList meshCells_;
+
+
+ // 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;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("sampledIsoSurfaceTopo");
+
+
+ // Constructors
+
+ //- Construct from dictionary
+ sampledIsoSurfaceTopo
+ (
+ const word& name,
+ const polyMesh& mesh,
+ const dictionary& dict
+ );
+
+
+ //- 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();
+
+
+ //- Points of surface
+ virtual const pointField& points() const
+ {
+ return MeshStorage::points();
+ }
+
+ //- Faces of surface
+ virtual const faceList& faces() const
+ {
+ return *this;
+ }
+
+ //- Const access to per-face zone/region information
+ virtual const labelList& zoneIds() const
+ {
+ return Foam::emptyLabelList;
+ }
+
+ //- Face area magnitudes
+ virtual const vectorField& Sf() const
+ {
+ return MeshStorage::Sf();
+ }
+
+ //- Face area magnitudes
+ virtual const scalarField& magSf() const
+ {
+ return MeshStorage::magSf();
+ }
+
+ //- Face centres
+ virtual const vectorField& Cf() const
+ {
+ return MeshStorage::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&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "sampledIsoSurfaceTopoTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
new file mode 100644
index 00000000000..5517417f476
--- /dev/null
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTopoTemplates.C
@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2018 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2018 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 "isoSurface.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
+
+ const labelList& elements = meshCells_;
+
+ // One value per point
+ auto tvalues = tmp<Field<Type>>::New(points().size());
+ auto& values = tvalues.ref();
+
+ const faceList& fcs = faces();
+ const pointField& pts = points();
+
+ bitSet pointDone(points().size());
+
+ forAll(faces(), cutFacei)
+ {
+ const face& f = fcs[cutFacei];
+ const label celli = elements[cutFacei];
+
+ for (const label pointi : f)
+ {
+ if (pointDone.set(pointi))
+ {
+ values[pointi] = interpolator.interpolate
+ (
+ pts[pointi],
+ celli
+ );
+ }
+ }
+ }
+
+ return tvalues;
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
new file mode 100644
index 00000000000..fb0531e423f
--- /dev/null
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -0,0 +1,1234 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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 "isoSurfaceTopo.H"
+#include "polyMesh.H"
+#include "tetMatcher.H"
+#include "tetPointRef.H"
+#include "DynamicField.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(isoSurfaceTopo, 0);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+bool Foam::isoSurfaceTopo::isTriCut
+(
+ const triFace& tri,
+ const scalarField& pointValues
+) const
+{
+ const bool aLower = (pointValues[tri[0]] < iso_);
+ const bool bLower = (pointValues[tri[1]] < iso_);
+ const bool cLower = (pointValues[tri[2]] < iso_);
+
+ return !(aLower == bLower && aLower == cLower);
+}
+
+
+Foam::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
+(
+ const bool isTet,
+ const label celli
+) const
+{
+ const cell& cFaces = mesh_.cells()[celli];
+
+ if (isTet)
+ {
+ for (const label facei : cFaces)
+ {
+ const face& f = mesh_.faces()[facei];
+
+ for (label fp = 1; fp < f.size() - 1; ++fp)
+ {
+ const triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
+
+ if (isTriCut(tri, pVals_))
+ {
+ return CUT;
+ }
+ }
+ }
+ return NOTCUT;
+ }
+ else
+ {
+ const bool cellLower = (cVals_[celli] < iso_);
+
+ // First check if there is any cut in cell
+ bool edgeCut = false;
+
+ for (const label facei : cFaces)
+ {
+ const face& f = mesh_.faces()[facei];
+
+ // Check pyramids cut
+ for (const label pointi : f)
+ {
+ if ((pVals_[pointi] < iso_) != cellLower)
+ {
+ edgeCut = true;
+ break;
+ }
+ }
+
+ if (edgeCut)
+ {
+ break;
+ }
+
+ const label fp0 = mesh_.tetBasePtIs()[facei];
+ label fp = f.fcIndex(fp0);
+ for (label i = 2; i < f.size(); ++i)
+ {
+ label nextFp = f.fcIndex(fp);
+
+ if (isTriCut(triFace(f[fp0], f[fp], f[nextFp]), pVals_))
+ {
+ edgeCut = true;
+ break;
+ }
+
+ fp = nextFp;
+ }
+
+ if (edgeCut)
+ {
+ break;
+ }
+ }
+
+ if (edgeCut)
+ {
+ // Count actual cuts (expensive since addressing needed)
+ // Note: not needed if you don't want to preserve maxima/minima
+ // centred around cellcentre. In that case just always return CUT
+
+ const labelList& cPoints = mesh_.cellPoints(celli);
+
+ label nPyrEdgeCuts = 0;
+
+ for (const label pointi : cPoints)
+ {
+ if ((pVals_[pointi] < iso_) != cellLower)
+ {
+ ++nPyrEdgeCuts;
+ }
+ }
+
+ if (nPyrEdgeCuts == cPoints.size())
+ {
+ return SPHERE;
+ }
+ else if (nPyrEdgeCuts)
+ {
+ return CUT;
+ }
+ }
+ return NOTCUT;
+ }
+}
+
+
+Foam::label Foam::isoSurfaceTopo::calcCutTypes
+(
+ tetMatcher& tet,
+ List<cellCutType>& cellCutTypes
+)
+{
+ const cellList& cells = mesh_.cells();
+
+ cellCutTypes.setSize(cells.size());
+ label nCutCells = 0;
+ forAll(cells, celli)
+ {
+ cellCutTypes[celli] = calcCutType(tet.isA(mesh_, celli), celli);
+
+ if (cellCutTypes[celli] == CUT)
+ {
+ ++nCutCells;
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "isoSurfaceCell : candidate cut cells "
+ << nCutCells << " / " << mesh_.nCells() << endl;
+ }
+ return nCutCells;
+}
+
+
+Foam::label Foam::isoSurfaceTopo::generatePoint
+(
+ const label facei,
+ const bool edgeIsDiag,
+ const edge& vertices,
+
+ DynamicList<edge>& pointToVerts,
+ DynamicList<label>& pointToFace,
+ DynamicList<bool>& pointFromDiag,
+ EdgeMap<label>& vertsToPoint
+) const
+{
+ EdgeMap<label>::const_iterator edgeFnd = vertsToPoint.find(vertices);
+ if (edgeFnd != vertsToPoint.end())
+ {
+ return edgeFnd();
+ }
+ else
+ {
+ // Generate new point
+ label pointi = pointToVerts.size();
+
+ pointToVerts.append(vertices);
+ pointToFace.append(facei);
+ pointFromDiag.append(edgeIsDiag);
+ vertsToPoint.insert(vertices, pointi);
+ return pointi;
+ }
+}
+
+
+void Foam::isoSurfaceTopo::generateTriPoints
+(
+ const label facei,
+ const FixedList<scalar, 4>& s,
+ const FixedList<point, 4>& p,
+ const FixedList<label, 4>& pIndex,
+ const FixedList<bool, 6>& edgeIsDiag,// Per tet edge whether is face diag
+
+ DynamicList<edge>& pointToVerts,
+ DynamicList<label>& pointToFace,
+ DynamicList<bool>& pointFromDiag,
+
+ EdgeMap<label>& vertsToPoint,
+ DynamicList<label>& verts // Every three verts is new triangle
+) const
+{
+ int triIndex = 0;
+ if (s[0] < iso_)
+ {
+ triIndex |= 1;
+ }
+ if (s[1] < iso_)
+ {
+ triIndex |= 2;
+ }
+ if (s[2] < iso_)
+ {
+ triIndex |= 4;
+ }
+ if (s[3] < iso_)
+ {
+ triIndex |= 8;
+ }
+
+ // Form the vertices of the triangles for each case
+ switch (triIndex)
+ {
+ case 0x00:
+ case 0x0F:
+ break;
+
+ case 0x01:
+ case 0x0E:
+ {
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[0],
+ edge(pIndex[0], pIndex[1]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[1],
+ edge(pIndex[0], pIndex[2]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[2],
+ edge(pIndex[0], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+
+ if (triIndex == 0x0E)
+ {
+ // Flip normals
+ label sz = verts.size();
+ Swap(verts[sz-2], verts[sz-1]);
+ }
+ }
+ break;
+
+ case 0x02:
+ case 0x0D:
+ {
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[0],
+ edge(pIndex[1], pIndex[0]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[3],
+ edge(pIndex[1], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[4],
+ edge(pIndex[1], pIndex[2]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+
+ if (triIndex == 0x0D)
+ {
+ // Flip normals
+ label sz = verts.size();
+ Swap(verts[sz-2], verts[sz-1]);
+ }
+ }
+ break;
+
+ case 0x03:
+ case 0x0C:
+ {
+ label p0p2
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[1],
+ edge(pIndex[0], pIndex[2]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ label p1p3
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[3],
+ edge(pIndex[1], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[2],
+ edge(pIndex[0], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append(p1p3);
+ verts.append(p0p2);
+ verts.append(p1p3);
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[4],
+ edge(pIndex[1], pIndex[2]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append(p0p2);
+
+ if (triIndex == 0x0C)
+ {
+ // Flip normals
+ label sz = verts.size();
+ Swap(verts[sz-5], verts[sz-4]);
+ Swap(verts[sz-2], verts[sz-1]);
+ }
+ }
+ break;
+
+ case 0x04:
+ case 0x0B:
+ {
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[1],
+ edge(pIndex[2], pIndex[0]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[4],
+ edge(pIndex[2], pIndex[1]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[5],
+ edge(pIndex[2], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+
+ if (triIndex == 0x0B)
+ {
+ // Flip normals
+ label sz = verts.size();
+ Swap(verts[sz-2], verts[sz-1]);
+ }
+ }
+ break;
+
+ case 0x05:
+ case 0x0A:
+ {
+ label p0p1
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[0],
+ edge(pIndex[0], pIndex[1]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ label p2p3
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[5],
+ edge(pIndex[2], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+
+ verts.append(p0p1);
+ verts.append(p2p3);
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[2],
+ edge(pIndex[0], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append(p0p1);
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[4],
+ edge(pIndex[1], pIndex[2]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append(p2p3);
+
+ if (triIndex == 0x0A)
+ {
+ // Flip normals
+ label sz = verts.size();
+ Swap(verts[sz-5], verts[sz-4]);
+ Swap(verts[sz-2], verts[sz-1]);
+ }
+ }
+ break;
+
+ case 0x06:
+ case 0x09:
+ {
+ label p0p1
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[0],
+ edge(pIndex[0], pIndex[1]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ label p2p3
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[5],
+ edge(pIndex[2], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+
+ verts.append(p0p1);
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[3],
+ edge(pIndex[1], pIndex[3]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append(p2p3);
+ verts.append(p0p1);
+ verts.append(p2p3);
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[1],
+ edge(pIndex[0], pIndex[2]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+
+ if (triIndex == 0x09)
+ {
+ // Flip normals
+ label sz = verts.size();
+ Swap(verts[sz-5], verts[sz-4]);
+ Swap(verts[sz-2], verts[sz-1]);
+ }
+ }
+ break;
+
+ case 0x08:
+ case 0x07:
+ {
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[2],
+ edge(pIndex[3], pIndex[0]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[5],
+ edge(pIndex[3], pIndex[2]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ verts.append
+ (
+ generatePoint
+ (
+ facei,
+ edgeIsDiag[3],
+ edge(pIndex[3], pIndex[1]),
+ pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ )
+ );
+ if (triIndex == 0x07)
+ {
+ // Flip normals
+ label sz = verts.size();
+ Swap(verts[sz-2], verts[sz-1]);
+ }
+ }
+ break;
+ }
+}
+
+
+void Foam::isoSurfaceTopo::generateTriPoints
+(
+ const polyMesh& mesh,
+ const label celli,
+ const bool isTet,
+
+ DynamicList<edge>& pointToVerts,
+ DynamicList<label>& pointToFace,
+ DynamicList<bool>& pointFromDiag,
+
+ EdgeMap<label>& vertsToPoint,
+ DynamicList<label>& verts,
+ DynamicList<label>& faceLabels
+) const
+{
+ const cell& cFaces = mesh.cells()[celli];
+
+ if (isTet)
+ {
+ // For tets don't do cell-centre decomposition, just use the
+ // tet points and values
+
+ label facei = cFaces[0];
+ const face& f0 = mesh_.faces()[facei];
+
+ // Get the other point
+ const face& f1 = mesh_.faces()[cFaces[1]];
+ label oppositeI = -1;
+ forAll(f1, fp)
+ {
+ oppositeI = f1[fp];
+ if (findIndex(f0, oppositeI) == -1)
+ {
+ break;
+ }
+ }
+
+
+ label p0 = f0[0];
+ label p1 = f0[1];
+ label p2 = f0[2];
+ FixedList<bool, 6> edgeIsDiag(false);
+
+ if (mesh.faceOwner()[facei] == celli)
+ {
+ Swap(p1, p2);
+ }
+
+ tetPointRef tet
+ (
+ mesh.points()[p0],
+ mesh.points()[p1],
+ mesh.points()[p2],
+ mesh.points()[oppositeI]
+ );
+
+ label startTrii = verts.size();
+ generateTriPoints
+ (
+ facei,
+ FixedList<scalar, 4>
+ ({
+ pVals_[p0],
+ pVals_[p1],
+ pVals_[p2],
+ pVals_[oppositeI]
+ }),
+ FixedList<point, 4>
+ ({
+ mesh.points()[p0],
+ mesh.points()[p1],
+ mesh.points()[p2],
+ mesh.points()[oppositeI]
+ }),
+ FixedList<label, 4>
+ ({
+ p0,
+ p1,
+ p2,
+ oppositeI
+ }),
+ edgeIsDiag,
+
+ pointToVerts,
+ pointToFace,
+ pointFromDiag,
+ vertsToPoint,
+ verts // Every three verts is new triangle
+ );
+
+ label nTris = (verts.size()-startTrii)/3;
+ for (label i = 0; i < nTris; ++i)
+ {
+ faceLabels.append(facei);
+ }
+ }
+ else
+ {
+ for (const label facei : cFaces)
+ {
+ const face& f = mesh.faces()[facei];
+
+ label fp0 = mesh.tetBasePtIs()[facei];
+
+ label startTrii = verts.size();
+
+ // Skip undefined tets
+ if (fp0 < 0)
+ {
+ fp0 = 0;
+ }
+
+ label fp = f.fcIndex(fp0);
+ for (label i = 2; i < f.size(); ++i)
+ {
+ label nextFp = f.fcIndex(fp);
+
+ FixedList<bool, 6> edgeIsDiag(false);
+
+ label p0 = f[fp0];
+ label p1 = f[fp];
+ label p2 = f[nextFp];
+ if (mesh.faceOwner()[facei] == celli)
+ {
+ Swap(p1, p2);
+ if (i != 2) edgeIsDiag[1] = true;
+ if (i != f.size()-1) edgeIsDiag[0] = true;
+ }
+ else
+ {
+ if (i != 2) edgeIsDiag[0] = true;
+ if (i != f.size()-1) edgeIsDiag[1] = true;
+ }
+
+ tetPointRef tet
+ (
+ mesh.points()[p0],
+ mesh.points()[p1],
+ mesh.points()[p2],
+ mesh.cellCentres()[celli]
+ );
+
+ generateTriPoints
+ (
+ facei,
+ FixedList<scalar, 4>
+ ({
+ pVals_[p0],
+ pVals_[p1],
+ pVals_[p2],
+ cVals_[celli]
+ }),
+ FixedList<point, 4>
+ ({
+ mesh.points()[p0],
+ mesh.points()[p1],
+ mesh.points()[p2],
+ mesh.cellCentres()[celli]
+ }),
+ FixedList<label, 4>
+ ({
+ p0,
+ p1,
+ p2,
+ mesh.nPoints()+celli
+ }),
+ edgeIsDiag,
+
+ pointToVerts,
+ pointToFace,
+ pointFromDiag,
+ vertsToPoint,
+ verts // Every three verts is new triangle
+ );
+
+ fp = nextFp;
+ }
+
+ label nTris = (verts.size()-startTrii)/3;
+ for (label i = 0; i < nTris; ++i)
+ {
+ faceLabels.append(facei);
+ }
+ }
+ }
+}
+
+
+void Foam::isoSurfaceTopo::triangulateOutside
+(
+ const bool filterDiag,
+ const primitivePatch& pp,
+ const boolList& pointFromDiag,
+ const labelList& pointToFace,
+ const label cellID,
+
+ DynamicList<face>& compactFaces,
+ DynamicList<label>& compactCellIDs
+) const
+{
+ // We can form pockets:
+ // - 1. triangle on face
+ // - 2. multiple triangles on interior (from diag edges)
+ // - the edge loop will be pocket since it is only the diag
+ // edges that give it volume?
+
+ // Retriangulate the exterior loops
+ const labelListList& edgeLoops = pp.edgeLoops();
+ const labelList& mp = pp.meshPoints();
+
+ for (const labelList& loop : edgeLoops)
+ {
+ if (loop.size() > 2)
+ {
+ compactFaces.append(face(0));
+ face& f = compactFaces.last();
+
+ f.setSize(loop.size());
+ label fpi = 0;
+ forAll(f, i)
+ {
+ label pointi = mp[loop[i]];
+ if (filterDiag && pointFromDiag[pointi])
+ {
+ label prevPointi = mp[loop[loop.fcIndex(i)]];
+ if
+ (
+ pointFromDiag[prevPointi]
+ && (pointToFace[pointi] != pointToFace[prevPointi])
+ )
+ {
+ f[fpi++] = pointi;
+ }
+ else
+ {
+ // Filter out diagonal point
+ }
+ }
+ else
+ {
+ f[fpi++] = pointi;
+ }
+ }
+
+ if (fpi > 2)
+ {
+ f.setSize(fpi);
+ }
+ else
+ {
+ // Keep original face
+ forAll(f, i)
+ {
+ label pointi = mp[loop[i]];
+ f[i] = pointi;
+ }
+ }
+ compactCellIDs.append(cellID);
+ }
+ }
+}
+
+
+Foam::MeshedSurface<Foam::face> Foam::isoSurfaceTopo::removeInsidePoints
+(
+ const bool filterDiag,
+ const MeshStorage& s,
+ const boolList& pointFromDiag,
+ const labelList& pointToFace,
+ const labelList& start, // Per cell the starting triangle
+ DynamicList<label>& pointCompactMap, // Per returned point the original
+ DynamicList<label>& compactCellIDs // Per returned tri the cellID
+) const
+{
+ const pointField& points = s.points();
+
+ pointCompactMap.clear();
+ compactCellIDs.clear();
+
+ DynamicList<face> compactFaces(s.size()/8);
+
+ for (label celli = 0; celli < start.size()-1; ++celli)
+ {
+ // All triangles of the current cell
+
+ label nTris = start[celli+1]-start[celli];
+
+ if (nTris)
+ {
+ const SubList<face> cellFaces(s, nTris, start[celli]);
+ const primitivePatch pp(cellFaces, points);
+
+ triangulateOutside
+ (
+ filterDiag,
+ pp,
+ pointFromDiag,
+ pointToFace,
+ //protectedFace,
+ celli,
+
+ compactFaces,
+ compactCellIDs
+ );
+ }
+ }
+
+
+ // Compact out unused points
+ // Pick up the used vertices
+ labelList oldToCompact(points.size(), -1);
+ DynamicField<point> compactPoints(points.size());
+ pointCompactMap.clear();
+
+ for (face& f : compactFaces)
+ {
+ forAll(f, fp)
+ {
+ label pointi = f[fp];
+ label compacti = oldToCompact[pointi];
+ if (compacti == -1)
+ {
+ compacti = compactPoints.size();
+ oldToCompact[pointi] = compacti;
+ compactPoints.append(points[pointi]);
+ pointCompactMap.append(pointi);
+ }
+ f[fp] = compacti;
+ }
+ }
+
+
+ MeshStorage cpSurface
+ (
+ std::move(compactPoints),
+ std::move(compactFaces),
+ s.surfZones()
+ );
+
+ return cpSurface;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::isoSurfaceTopo::isoSurfaceTopo
+(
+ const polyMesh& mesh,
+ const scalarField& cVals,
+ const scalarField& pVals,
+ const scalar iso,
+ const filterType filter
+)
+:
+ mesh_(mesh),
+ cVals_(cVals),
+ pVals_(pVals),
+ iso_(iso)
+{
+ if (debug)
+ {
+ Pout<< "isoSurfaceTopo : iso:" << iso_ << " filter:" << filter << endl;
+ }
+
+ tetMatcher tet;
+
+ // Determine if any cut through cell
+ List<cellCutType> cellCutTypes;
+ const label nCutCells = calcCutTypes(tet, cellCutTypes);
+
+ // Per cell: 5 pyramids cut, each generating 2 triangles
+ // - pointToVerts : from generated iso point to originating mesh verts
+ DynamicList<edge> pointToVerts(10*nCutCells);
+ // - pointToFace : from generated iso point to originating mesh face
+ DynamicList<label> pointToFace(10*nCutCells);
+ // - pointToFace : from generated iso point whether is on face diagonal
+ DynamicList<bool> pointFromDiag(10*nCutCells);
+
+ // Per cell: number of intersected edges:
+ // - four faces cut so 4 mesh edges + 4 face-diagonal edges
+ // - 4 of the pyramid edges
+ EdgeMap<label> vertsToPoint(12*nCutCells);
+ DynamicList<label> verts(12*nCutCells);
+ // Per cell: 5 pyramids cut (since only one pyramid not cut)
+ DynamicList<label> faceLabels(5*nCutCells);
+ DynamicList<label> cellLabels(5*nCutCells);
+
+
+ labelList startTri(mesh_.nCells()+1, 0);
+
+ for (label celli = 0; celli < mesh_.nCells(); ++celli)
+ {
+ startTri[celli] = faceLabels.size();
+ if (cellCutTypes[celli] != NOTCUT)
+ {
+ generateTriPoints
+ (
+ mesh,
+ celli,
+ tet.isA(mesh_, celli),
+
+ pointToVerts,
+ pointToFace,
+ pointFromDiag,
+
+ vertsToPoint,
+ verts,
+ faceLabels
+ );
+
+ for (label i = startTri[celli]; i < faceLabels.size(); ++i)
+ {
+ cellLabels.append(celli);
+ }
+ }
+ }
+ startTri[mesh_.nCells()] = faceLabels.size();
+
+
+ pointToVerts_.transfer(pointToVerts);
+ meshCells_.transfer(cellLabels);
+ pointToFace_.transfer(pointToFace);
+
+ tmp<pointField> allPoints
+ (
+ interpolate
+ (
+ mesh_.cellCentres(),
+ mesh_.points()
+ )
+ );
+
+
+ // Assign to MeshedSurface
+ faceList allTris(faceLabels.size());
+ label verti = 0;
+ for (face& allTri : allTris)
+ {
+ allTri.setSize(3);
+ allTri[0] = verts[verti++];
+ allTri[1] = verts[verti++];
+ allTri[2] = verts[verti++];
+ }
+
+
+ surfZoneList allZones(1);
+ allZones[0] = surfZone
+ (
+ "allFaces",
+ allTris.size(), // Size
+ 0, // Start
+ 0 // Index
+ );
+
+ MeshStorage updated
+ (
+ std::move(allPoints),
+ std::move(allTris),
+ std::move(allZones)
+ );
+ MeshStorage::transfer(updated);
+
+ // Now:
+ // - generated faces and points are assigned to *this
+ // - per point we know:
+ // - pointOnDiag: whether it is on a face-diagonal edge
+ // - pointToFace_: from what pyramid (cell+face) it was produced
+ // (note that the pyramid faces are shared between multiple mesh faces)
+ // - pointToVerts_ : originating mesh vertex or cell centre
+
+
+ if (debug)
+ {
+ Pout<< "isoSurfaceTopo : generated " << size() << " faces." << endl;
+ }
+
+
+ if (filter != NONE)
+ {
+ // Triangulate outside (filter edges to cell centres and optionally
+ // face diagonals)
+ DynamicList<label> pointCompactMap(size()); // Back to original point
+ DynamicList<label> compactCellIDs(size()); // Per tri the cell
+ MeshStorage::operator=
+ (
+ removeInsidePoints
+ (
+ (filter == DIAGCELL ? true : false),
+ *this,
+ pointFromDiag,
+ pointToFace_,
+ startTri,
+ pointCompactMap,
+ compactCellIDs
+ )
+ );
+
+ pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointCompactMap)();
+ pointToFace_ = UIndirectList<label>(pointToFace_, pointCompactMap)();
+ pointFromDiag = UIndirectList<bool>(pointFromDiag, pointCompactMap)();
+ meshCells_.transfer(compactCellIDs);
+
+ if (debug)
+ {
+ Pout<< "isoSurfaceTopo :"
+ << " after removing cell centre and face-diag triangles : "
+ << size() << endl;
+ }
+
+
+ if (filter == DIAGCELL)
+ {
+ // We remove verts on face diagonals. This is in fact just
+ // straightening the edges of the face through the cell. This can
+ // close off 'pockets' of triangles and create open or
+ // multiply-connected triangles
+
+ // Solved by eroding open-edges
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+ // Mark points on mesh outside. Note that we extend with nCells
+ // so we can easily index with pointToVerts_.
+ PackedBoolList isBoundaryPoint(mesh.nPoints() + mesh.nCells());
+ for
+ (
+ label facei = mesh.nInternalFaces();
+ facei < mesh.nFaces();
+ ++facei
+ )
+ {
+ isBoundaryPoint.set(mesh.faces()[facei]);
+ }
+
+
+ while (true)
+ {
+ const labelList& mp = meshPoints();
+
+ PackedBoolList removeFace(this->size());
+ label nFaces = 0;
+ {
+ const labelListList& edgeFaces =
+ MeshStorage::edgeFaces();
+ forAll(edgeFaces, edgei)
+ {
+ const labelList& eFaces = edgeFaces[edgei];
+ if (eFaces.size() == 1)
+ {
+ // Open edge. Check that vertices do not originate
+ // from a boundary face
+ const edge& e = edges()[edgei];
+ const edge& verts0 = pointToVerts_[mp[e[0]]];
+ const edge& verts1 = pointToVerts_[mp[e[1]]];
+ if
+ (
+ isBoundaryPoint[verts0[0]]
+ && isBoundaryPoint[verts0[1]]
+ && isBoundaryPoint[verts1[0]]
+ && isBoundaryPoint[verts1[1]]
+ )
+ {
+ // Open edge on boundary face. Keep
+ }
+ else
+ {
+ // Open edge. Mark for erosion
+ if (removeFace.set(eFaces[0]))
+ {
+ ++nFaces;
+ }
+ }
+ }
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "isoSurfaceTopo :"
+ << " removing " << nFaces
+ << " faces since on open edges" << endl;
+ }
+
+ if (returnReduce(nFaces, sumOp<label>()) == 0)
+ {
+ break;
+ }
+
+ // Remove the faces
+ labelHashSet keepFaces(2*size());
+ forAll(removeFace, facei)
+ {
+ if (!removeFace[facei])
+ {
+ keepFaces.insert(facei);
+ }
+ }
+
+ labelList pointMap;
+ labelList faceMap;
+ MeshStorage filteredSurf
+ (
+ MeshStorage::subsetMesh
+ (
+ keepFaces,
+ pointMap,
+ faceMap
+ )
+ );
+ MeshStorage::transfer(filteredSurf);
+
+ pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
+ pointToFace_ = UIndirectList<label>(pointToFace_, pointMap)();
+ pointFromDiag = UIndirectList<bool>(pointFromDiag, pointMap)();
+ meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.H b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
new file mode 100644
index 00000000000..322eab8d7c4
--- /dev/null
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
@@ -0,0 +1,266 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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::isoSurfaceTopo
+
+Description
+ Marching tet iso surface algorithm with optional filtering to keep only
+ points originating from mesh edges.
+
+SourceFiles
+ isoSurfaceTopo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef isoSurfaceTopo_H
+#define isoSurfaceTopo_H
+
+#include "labelPair.H"
+#include "pointIndexHit.H"
+#include "PackedBoolList.H"
+#include "MeshedSurface.H"
+#include "edgeList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+class tetMatcher;
+
+/*---------------------------------------------------------------------------*\
+ Class isoSurfaceTopo Declaration
+\*---------------------------------------------------------------------------*/
+
+class isoSurfaceTopo
+:
+ public MeshedSurface<face>
+{
+ // Private typedefs for convenience
+ typedef MeshedSurface<face> MeshStorage;
+
+public:
+
+ enum filterType
+ {
+ NONE, // No filtering
+ DIAGCELL, // Remove points from face-diagonal and pyramid
+ // (vertex to cell-centre) edges
+ CELL // Only remove points from pyramid edges
+ };
+
+
+private:
+
+ // Private data
+
+ enum cellCutType
+ {
+ NOTCUT, // Not cut
+ SPHERE, // All edges to cell centre cut
+ CUT // Normal cut
+ };
+
+
+ //- Reference to mesh
+ const polyMesh& mesh_;
+
+ const scalarField& cVals_;
+
+ const scalarField& pVals_;
+
+ //- Iso value
+ const scalar iso_;
+
+ //- Per point: originating mesh vertex/cc. See encoding above
+ edgeList pointToVerts_;
+
+ //- For every face the original cell in mesh
+ labelList meshCells_;
+
+ //- For every point the originating face in mesh
+ labelList pointToFace_;
+
+
+ // Private Member Functions
+
+ //- Does any edge of triangle cross iso value?
+ bool isTriCut
+ (
+ const triFace& tri,
+ const scalarField& pointValues
+ ) const;
+
+ //- Determine whether cell is cut
+ cellCutType calcCutType
+ (
+ const bool isTet,
+ const label
+ ) const;
+
+ //- Determine for all mesh whether cell is cut
+ label calcCutTypes
+ (
+ tetMatcher& tet,
+ List<cellCutType>& cellCutTypes
+ );
+
+ //- Generate single point on edge
+ label generatePoint
+ (
+ const label facei,
+ const bool edgeIsDiag,
+ const edge& vertices,
+
+ DynamicList<edge>& pointToVerts,
+ DynamicList<label>& pointToFace,
+ DynamicList<bool>& pointFromDiag,
+ EdgeMap<label>& vertsToPoint
+ ) const;
+
+ //- Generate triangles from tet
+ void generateTriPoints
+ (
+ const label facei,
+ const FixedList<scalar, 4>& s,
+ const FixedList<point, 4>& p,
+ const FixedList<label, 4>& pIndex,
+ const FixedList<bool, 6>& edgeIsDiag,
+
+ DynamicList<edge>& pointToVerts,
+ DynamicList<label>& pointToFace,
+ DynamicList<bool>& pointFromDiag,
+
+ EdgeMap<label>& vertsToPoint,
+ DynamicList<label>& verts
+ ) const;
+
+ //- Generate triangles from cell
+ void generateTriPoints
+ (
+ const polyMesh& mesh,
+ const label celli,
+ const bool isTet,
+
+ DynamicList<edge>& pointToVerts,
+ DynamicList<label>& pointToFace,
+ DynamicList<bool>& pointFromDiag,
+
+ EdgeMap<label>& vertsToPoint,
+ DynamicList<label>& verts,
+ DynamicList<label>& faceLabels
+ ) const;
+
+
+ // Simplification
+
+ void triangulateOutside
+ (
+ const bool filterDiag,
+ const PrimitivePatch<face, SubList, const pointField&>& pp,
+ const boolList& pointFromDiag,
+ const labelList& pointToFace,
+ const label cellID,
+
+ DynamicList<face>& compactFaces,
+ DynamicList<label>& compactCellIDs
+ ) const;
+
+ MeshStorage removeInsidePoints
+ (
+ const bool filterDiag,
+ const MeshStorage& s,
+ const boolList& pointFromDiag,
+ const labelList& pointToFace,
+ const labelList& start, // Per cell:starting tri
+ DynamicList<label>& pointCompactMap, // Per point the original
+ DynamicList<label>& compactCellIDs // Per face the cellID
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("isoSurfaceTopo");
+
+
+ // Constructors
+
+ //- Construct from dictionary
+ isoSurfaceTopo
+ (
+ const polyMesh& mesh,
+ const scalarField& cellValues,
+ const scalarField& pointValues,
+ const scalar iso,
+ const filterType filter = DIAGCELL
+ );
+
+
+ // Member Functions
+
+ //- For every face original cell in mesh
+ const labelList& meshCells() const
+ {
+ return meshCells_;
+ }
+
+ //- For every point originating face (pyramid) in mesh
+ const labelList& pointToFace() const
+ {
+ return pointToFace_;
+ }
+
+ //- Per point: originating mesh vertex/cc. See encoding above�
+ const edgeList& pointToVerts() const
+ {
+ return pointToVerts_;
+ }
+
+ //- Interpolates cCoords,pCoords.
+ template<class Type>
+ tmp<Field<Type>> interpolate
+ (
+ const Field<Type>& cCoords,
+ const Field<Type>& pCoords
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "isoSurfaceTopoTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
new file mode 100644
index 00000000000..631d948e4c2
--- /dev/null
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::isoSurfaceTopo::interpolate
+(
+ const Field<Type>& cellCoords,
+ const Field<Type>& pointCoords
+) const
+{
+ tmp<Field<Type>> tfld(new Field<Type>(pointToVerts_.size()));
+ Field<Type>& fld = tfld.ref();
+
+ forAll(pointToVerts_, i)
+ {
+ scalar s0;
+ Type p0;
+ {
+ label v0 = pointToVerts_[i][0];
+ if (v0 < mesh_.nPoints())
+ {
+ s0 = pVals_[v0];
+ p0 = pointCoords[v0];
+ }
+ else
+ {
+ label celli = v0-mesh_.nPoints();
+ s0 = cVals_[celli];
+ p0 = cellCoords[celli];
+ }
+ }
+
+ scalar s1;
+ Type p1;
+ {
+ label v1 = pointToVerts_[i][1];
+ if (v1 < mesh_.nPoints())
+ {
+ s1 = pVals_[v1];
+ p1 = pointCoords[v1];
+ }
+ else
+ {
+ label celli = v1-mesh_.nPoints();
+ s1 = cVals_[celli];
+ p1 = cellCoords[celli];
+ }
+ }
+
+ scalar d = s1-s0;
+ if (mag(d) > VSMALL)
+ {
+ scalar s = (iso_-s0)/d;
+ fld[i] = s*p1+(1.0-s)*p0;
+ }
+ else
+ {
+ fld[i] = 0.5*(p0+p1);
+ }
+ }
+
+ return tfld;
+}
+
+
+// ************************************************************************* //
--
GitLab