Merge branch 'master' of github.com-OpenFOAM:OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry Weller
2019-07-15 12:51:00 +01:00
9 changed files with 140 additions and 121 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -467,7 +467,7 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::Eo
)
{
const vector dir = U/max(mag(U), rootVSmall);
return mag(g & dir)*(rho - rhoc)*sqr(d)/max(sigma, rootVSmall);
return mag(g & dir)*mag(rho - rhoc)*sqr(d)/max(sigma, rootVSmall);
}

View File

@ -269,7 +269,7 @@ void Foam::sampledSurfaces::distanceSurface::createGeometry()
cellDistance,
pointDistance_,
distance_,
regularise_ ? isoSurface::DIAGCELL : isoSurface::NONE
filter_
)
);
@ -310,7 +310,12 @@ Foam::sampledSurfaces::distanceSurface::distanceSurface
),
distance_(readScalar(dict.lookup("distance"))),
signed_(readBool(dict.lookup("signed"))),
regularise_(dict.lookupOrDefault("regularise", true)),
filter_
(
dict.found("filtering")
? isoSurface::filterTypeNames_.read(dict.lookup("filtering"))
: isoSurface::filterType::full
),
average_(dict.lookupOrDefault("average", false)),
zoneKey_(keyType::null),
needsUpdate_(true),
@ -327,7 +332,7 @@ Foam::sampledSurfaces::distanceSurface::distanceSurface
const word& surfaceName,
const scalar distance,
const bool signedDistance,
const Switch regularise,
const isoSurface::filterType filter,
const Switch average
)
:
@ -351,7 +356,7 @@ Foam::sampledSurfaces::distanceSurface::distanceSurface
),
distance_(distance),
signed_(signedDistance),
regularise_(regularise),
filter_(filter),
average_(average),
zoneKey_(keyType::null),
needsUpdate_(true),

View File

@ -66,7 +66,7 @@ class distanceSurface
const bool signed_;
//- Whether to coarsen
const Switch regularise_;
const isoSurface::filterType filter_;
//- Whether to recalculate cell values as average of point values
const Switch average_;
@ -132,7 +132,7 @@ public:
const word& surfaceName,
const scalar distance,
const bool signedDistance,
const Switch regularise,
const isoSurface::filterType filter,
const Switch average
);

View File

@ -38,6 +38,15 @@ namespace Foam
defineTypeNameAndDebug(isoSurface, 0);
}
namespace Foam
{
template<>
const char* NamedEnum<isoSurface::filterType, 3>::names[] =
{"none", "partial", "full"};
}
const Foam::NamedEnum<Foam::isoSurface::filterType, 3>
Foam::isoSurface::filterTypeNames_;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -75,11 +84,11 @@ Foam::isoSurface::cellCutType Foam::isoSurface::calcCutType
if (isTriCut(tri, pVals_))
{
return CUT;
return cellCutType::cut;
}
}
}
return NOTCUT;
return cellCutType::notCut;
}
else
{
@ -108,13 +117,7 @@ Foam::isoSurface::cellCutType Foam::isoSurface::calcCutType
break;
}
label fp0 = tetBasePtIs_[facei];
// Fall back for problem decompositions
if (fp0 < 0)
{
fp0 = 0;
}
const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
label fp = f.fcIndex(fp0);
for (label i = 2; i < f.size(); i++)
@ -140,7 +143,7 @@ Foam::isoSurface::cellCutType Foam::isoSurface::calcCutType
{
// 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
// centred around cellcentre. In that case just always return cut
const labelList& cPoints = mesh_.cellPoints(celli);
@ -156,16 +159,16 @@ Foam::isoSurface::cellCutType Foam::isoSurface::calcCutType
if (nPyrCuts == cPoints.size())
{
return SPHERE;
return cellCutType::sphere;
}
else
{
return CUT;
return cellCutType::cut;
}
}
else
{
return NOTCUT;
return cellCutType::notCut;
}
}
}
@ -185,7 +188,7 @@ Foam::label Foam::isoSurface::calcCutTypes
{
cellCutTypes[celli] = calcCutType(tet.isA(mesh_, celli), celli);
if (cellCutTypes[celli] == CUT)
if (cellCutTypes[celli] == cellCutType::cut)
{
nCutCells++;
}
@ -247,37 +250,31 @@ void Foam::isoSurface::fixTetBasePtIs()
tetBasePtIs_ = mesh_.tetBasePtIs();
// Pre-filter: mark all cells with illegal base points
labelHashSet problemCells(cells.size()/128);
// Mark all cells with illegal base points as potentially problematic
PackedBoolList problemCells(cells.size(), false);
forAll(tetBasePtIs_, facei)
{
if (tetBasePtIs_[facei] == -1)
{
problemCells.insert(faceOwner[facei]);
problemCells[faceOwner[facei]] = true;
if (mesh_.isInternalFace(facei))
{
problemCells.insert(faceNeighbour[facei]);
problemCells[faceNeighbour[facei]] = true;
}
}
}
label nAdapted = 0;
// Number of times a point occurs in a cell. Used to detect dangling
// vertices (count = 2)
Map<label> pointCount;
// Analyse problem cells for points shared by two faces only
// Mark all points which are shared by just two faces within an adjacent
// problem cell as problematic
PackedBoolList problemPoints(mesh_.points().size(), false);
forAll(cells, celli)
{
if (problemCells.found(celli))
if (problemCells[celli])
{
const cell& cFaces = cells[celli];
pointCount.clear();
Map<label> pointCount;
forAll(cFaces, i)
{
const label facei = cFaces[i];
@ -293,13 +290,11 @@ void Foam::isoSurface::fixTetBasePtIs()
}
else
{
++pointFnd();
++ pointFnd();
}
}
}
// Check for any points with count 2
bool haveDangling = false;
forAllConstIter(Map<label>, pointCount, iter)
{
if (iter() == 1)
@ -308,63 +303,75 @@ void Foam::isoSurface::fixTetBasePtIs()
<< " at:" << mesh_.points()[iter.key()]
<< " only used by one face" << exit(FatalError);
}
else if (iter() == 2)
if (iter() == 2)
{
haveDangling = true;
break;
problemPoints[iter.key()] = true;
}
}
}
}
if (haveDangling)
// For all faces which form a part of a problem-cell, check if the base
// point is adjacent to any problem points. If it is, re-calculate the base
// point so that it is not.
label nAdapted = 0;
forAll(tetBasePtIs_, facei)
{
if
(
problemCells[faceOwner[facei]]
|| (mesh_.isInternalFace(facei) && problemCells[faceNeighbour[facei]])
)
{
const face& f = faces[facei];
// Check if either of the points adjacent to the base point is a
// problem point. If not, the existing base point can be retained.
const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
const bool prevPointIsProblem = problemPoints[f[f.rcIndex(fp0)]];
const bool nextPointIsProblem = problemPoints[f[f.fcIndex(fp0)]];
if (!prevPointIsProblem && !nextPointIsProblem)
{
// Any point next to a dangling point should not be used
// as the fan base since this would cause two duplicate
// triangles.
forAll(cFaces, i)
continue;
}
// A new base point is required. Pick the point that results in the
// least-worst-worst tet, and which is not adjacent to any problem
// points.
scalar maxQ = -GREAT;
label maxFp = -1;
forAll(f, fp)
{
const bool prevPointIsProblem = problemPoints[f[f.rcIndex(fp)]];
const bool nextPointIsProblem = problemPoints[f[f.fcIndex(fp)]];
if (!prevPointIsProblem && !nextPointIsProblem)
{
const label facei = cFaces[i];
if (tetBasePtIs_[facei] == -1)
const scalar q = minTetQ(facei, fp);
if (q > maxQ)
{
const face& f = faces[facei];
// All the possible base points cause negative tets.
// Choose the least-worst one
scalar maxQ = -GREAT;
label maxFp = -1;
label prevCount = pointCount[f.last()];
forAll(f, fp)
{
label nextCount = pointCount[f[f.fcIndex(fp)]];
if (prevCount > 2 && nextCount > 2)
{
const scalar q = minTetQ(facei, fp);
if (q > maxQ)
{
maxQ = q;
maxFp = fp;
}
}
prevCount = pointCount[f[fp]];
}
if (maxFp != -1)
{
// Least worst base point
tetBasePtIs_[facei] = maxFp;
}
else
{
// No point found on face that would not result
// in some duplicate triangle. Very rare. Do what?
tetBasePtIs_[facei] = 0;
}
nAdapted++;
maxQ = q;
maxFp = fp;
}
}
}
if (maxFp != -1)
{
// Success! Set the new base point
tetBasePtIs_[facei] = maxFp;
}
else
{
// No point was found on face that would not result in some
// duplicate triangle. Do what? Continue and hope? Spit an
// error? Silently or noisily reduce the filtering level?
}
++ nAdapted;
}
}
@ -373,8 +380,6 @@ void Foam::isoSurface::fixTetBasePtIs()
Pout<< "isoSurface : adapted starting point of triangulation on "
<< nAdapted << " faces." << endl;
}
syncTools::syncFaceList(mesh_, tetBasePtIs_, maxEqOp<label>());
}
@ -908,16 +913,10 @@ void Foam::isoSurface::generateTriPoints
label facei = cFaces[cFacei];
const face& f = faces[facei];
label fp0 = tetBasePtIs_[facei];
const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
label startTrii = verts.size();
// Fallback
if (fp0 < 0)
{
fp0 = 0;
}
label fp = f.fcIndex(fp0);
for (label i = 2; i < f.size(); i++)
{
@ -1171,7 +1170,8 @@ Foam::isoSurface::isoSurface
{
if (debug)
{
Pout<< "isoSurface : iso:" << iso_ << " filter:" << filter << endl;
Pout<< "isoSurface : iso:" << iso_
<< " filter:" << filterTypeNames_[filter] << endl;
}
fixTetBasePtIs();
@ -1205,7 +1205,7 @@ Foam::isoSurface::isoSurface
for (label celli = 0; celli < mesh_.nCells(); celli++)
{
startTri[celli] = faceLabels.size();
if (cellCutTypes[celli] != NOTCUT)
if (cellCutTypes[celli] != cellCutType::notCut)
{
generateTriPoints
(
@ -1288,7 +1288,7 @@ Foam::isoSurface::isoSurface
}
if (filter != NONE)
if (filter != filterType::none)
{
// Triangulate outside (filter edges to cell centres and optionally
// face diagonals)
@ -1298,7 +1298,7 @@ Foam::isoSurface::isoSurface
(
removeInsidePoints
(
(filter == DIAGCELL ? true : false),
(filter == filterType::full ? true : false),
*this,
pointFromDiag,
pointToFace_,
@ -1321,7 +1321,7 @@ Foam::isoSurface::isoSurface
}
if (filter == DIAGCELL)
if (filter == filterType::full)
{
// We remove verts on face diagonals. This is in fact just
// straightening the edges of the face through the cell. This can

View File

@ -41,6 +41,7 @@ SourceFiles
#include "PackedBoolList.H"
#include "MeshedSurface.H"
#include "edgeList.H"
#include "NamedEnum.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -60,24 +61,28 @@ class isoSurface
{
public:
enum filterType
enum class filterType
{
NONE, // No filtering
DIAGCELL, // Remove points from face-diagonal and pyramid
// (vertex to cell-centre) edges
CELL // Only remove points from pyramid edges
none, // No filtering
partial, // Remove points from vertex to cell-centre edges and
// merge triangles that form a contiguous cut through a
// single cell
full // Also remove points from face-diagonals and merge
// edges that originate from the same face
};
static const NamedEnum<filterType, 3> filterTypeNames_;
private:
// Private Data
enum cellCutType
enum class cellCutType
{
NOTCUT, // Not cut
SPHERE, // All edges to cell centre cut
CUT // Normal cut
notCut, // Not cut
sphere, // All edges to cell centre cut
cut // Normal cut
};
@ -222,7 +227,7 @@ public:
const scalarField& cellValues,
const scalarField& pointValues,
const scalar iso,
const filterType filter = DIAGCELL
const filterType filter
);

View File

@ -24,7 +24,6 @@ License
\*---------------------------------------------------------------------------*/
#include "sampledIsoSurface.H"
#include "isoSurface.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -121,9 +120,7 @@ bool Foam::sampledSurfaces::isoSurface::updateGeometry() const
cellFld.primitiveField(),
pointFld().primitiveField(),
isoVals_[isoi],
regularise_
? Foam::isoSurface::DIAGCELL
: Foam::isoSurface::NONE
filter_
)
);
}
@ -215,7 +212,8 @@ bool Foam::sampledSurfaces::isoSurface::updateGeometry() const
Pout<< "sampledSurfaces::isoSurface::updateGeometry() : "
"constructed iso:"
<< nl
<< " regularise : " << regularise_ << nl
<< " filtering : "
<< Foam::isoSurface::filterTypeNames_[filter_] << nl
<< " isoField : " << isoField_ << nl;
if (isoVals_.size() == 1)
{
@ -251,7 +249,12 @@ Foam::sampledSurfaces::isoSurface::isoSurface
? scalarField(dict.lookup("isoValues"))
: scalarField(1, readScalar(dict.lookup("isoValue")))
),
regularise_(dict.lookupOrDefault("regularise", true)),
filter_
(
dict.found("filtering")
? Foam::isoSurface::filterTypeNames_.read(dict.lookup("filtering"))
: Foam::isoSurface::filterType::full
),
zoneKey_(keyType::null),
prevTimeIndex_(-1),
meshCells_(0)

View File

@ -38,6 +38,7 @@ SourceFiles
#define sampledIsoSurface_H
#include "sampledSurface.H"
#include "isoSurface.H"
#include "MeshedSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -65,7 +66,7 @@ class isoSurface
const scalarField isoVals_;
//- Whether to coarsen
const Switch regularise_;
const Foam::isoSurface::filterType filter_;
//- If restricted to zones, name of this zone or a regular expression
keyType zoneKey_;

View File

@ -225,7 +225,7 @@ void Foam::sampledSurfaces::cuttingPlane::createGeometry()
cellDistance,
pointDistance_,
0,
regularise_ ? isoSurface::DIAGCELL : isoSurface::NONE
filter_
)
);
@ -248,7 +248,12 @@ Foam::sampledSurfaces::cuttingPlane::cuttingPlane
:
sampledSurface(name, mesh, dict),
plane_(dict),
regularise_(dict.lookupOrDefault("regularise", true)),
filter_
(
dict.found("filtering")
? isoSurface::filterTypeNames_.read(dict.lookup("filtering"))
: isoSurface::filterType::full
),
average_(dict.lookupOrDefault("average", false)),
zoneID_(dict.lookupOrDefault("zone", word::null), mesh.cellZones()),
exposedPatchName_(word::null),

View File

@ -63,7 +63,7 @@ class cuttingPlane
const plane plane_;
//- Whether to coarsen
const Switch regularise_;
const isoSurface::filterType filter_;
//- Whether to recalculate cell values as average of point values
const Switch average_;