Skip to content
GitLab
Commit f0f677e3 authored by Mark Olesen's avatar Mark Olesen
Browse files

Merge commit 'OpenCFD/master' into olesenm

parents 55d4fb2e bdb4349e
Hide whitespace changes
Inline Side-by-side
applications/utilities/postProcessing/sampling/sample/sampleDict
View file @ f0f677e3
......@@ -161,6 +161,10 @@ surfaces
isoField rho;
isoValue 0.5;
interpolate true;
//zone ABC; // Optional: zone only
//exposedPatchName fixedWalls; // Optional: zone only
// regularise false; // Optional: do not simplify
}
constantIso
......@@ -171,7 +175,7 @@ surfaces
isoField rho;
isoValue 0.5;
interpolate false;
// regularise false; // Optional: do not simplify
regularise false; // do not simplify
}
);
......
src/sampling/sampledSurface/isoSurface/isoSurface.C
View file @ f0f677e3
......@@ -206,7 +206,8 @@ void Foam::isoSurface::calcCutTypes
if (debug)
{
Pout<< "isoSurface : detected " << nCutCells_
<< " candidate cut cells." << endl;
<< " candidate cut cells (out of " << mesh_.nCells()
<< ")." << endl;
}
}
......@@ -814,65 +815,53 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
tris.transfer(dynTris);
}
if (debug)
{
Pout<< "isoSurface : merged from " << nTris
<< " down to " << tris.size() << " triangles." << endl;
}
// Use face centres to determine 'flat hole' situation (see RMT paper).
// Determine 'flat hole' situation (see RMT paper).
// Two unconnected triangles get connected because (some of) the edges
// separating them get collapsed. Below only checks for duplicate triangles,
// not non-manifold edge connectivity.
if (checkDuplicates)
{
if (debug)
{
Pout<< "isoSurface : merged from " << nTris
<< " down to " << tris.size() << " triangles." << endl;
}
labelListList pointFaces;
invertManyToMany(newPoints.size(), tris, pointFaces);
// Filter out duplicates.
DynamicList<label> newToOldTri(tris.size());
pointField centres(tris.size());
forAll(tris, triI)
{
centres[triI] = tris[triI].centre(newPoints);
}
pointField mergedCentres;
labelList oldToMerged;
bool hasMerged = mergePoints
(
centres,
mergeDistance_,
false,
oldToMerged,
mergedCentres
);
const labelledTri& tri = tris[triI];
if (debug)
{
Pout<< "isoSurface : detected "
<< centres.size()-mergedCentres.size()
<< " duplicate triangles." << endl;
}
const labelList& pFaces = pointFaces[tri[0]];
if (hasMerged)
{
// Filter out duplicates.
label newTriI = 0;
DynamicList<label> newToOldTri(tris.size());
labelList newToMaster(mergedCentres.size(), -1);
forAll(tris, triI)
// Find the minimum of any duplicates
label dupTriI = -1;
forAll(pFaces, i)
{
label mergedI = oldToMerged[triI];
if (newToMaster[mergedI] == -1)
if (pFaces[i] < triI && tris[pFaces[i]] == tri)
{
newToMaster[mergedI] = triI;
newToOldTri.append(triMap[triI]);
tris[newTriI++] = tris[triI];
dupTriI = pFaces[i];
}
}
triMap.transfer(newToOldTri);
tris.setSize(newTriI);
if (dupTriI == -1)
{
// There is no lower triangle
label newTriI = newToOldTri.size();
newToOldTri.append(triI);
tris[newTriI] = tris[triI];
}
}
triMap.transfer(newToOldTri);
tris.setSize(triMap.size());
}
return triSurface(tris, geometricSurfacePatchList(0), newPoints, true);
......@@ -986,7 +975,7 @@ void Foam::isoSurface::calcAddressing
{
Pout<< "isoSurface : detected "
<< mergedCentres.size()
<< " edges on " << surf.size() << " triangles." << endl;
<< " geometric edges on " << surf.size() << " triangles." << endl;
}
if (!hasMerged)
......@@ -1013,6 +1002,10 @@ void Foam::isoSurface::calcAddressing
edgeFace1 = -1;
edgeFacesRest.clear();
// Overflow edge faces for geometric shared edges that turned
// out to be different anyway.
EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
forAll(oldToMerged, oldEdgeI)
{
label triI = oldEdgeI / 3;
......@@ -1020,33 +1013,127 @@ void Foam::isoSurface::calcAddressing
if (edgeFace0[edgeI] == -1)
{
// First triangle for edge
edgeFace0[edgeI] = triI;
}
else if (edgeFace1[edgeI] == -1)
{
edgeFace1[edgeI] = triI;
}
else
{
//WarningIn("orientSurface(triSurface&)")
// << "Edge " << edgeI << " with centre " << mergedCentres[edgeI]
// << " used by more than two triangles: " << edgeFace0[edgeI]
// << ", "
// << edgeFace1[edgeI] << " and " << triI << endl;
Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
if (iter != edgeFacesRest.end())
//- Check that the two triangles actually topologically
// share an edge
const labelledTri& prevTri = surf[edgeFace0[edgeI]];
const labelledTri& tri = surf[triI];
label fp = oldEdgeI % 3;
edge e(tri[fp], tri[tri.fcIndex(fp)]);
label prevTriIndex = -1;
forAll(prevTri, i)
{
labelList& eFaces = iter();
label sz = eFaces.size();
eFaces.setSize(sz+1);
eFaces[sz] = triI;
if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
{
prevTriIndex = i;
break;
}
}
if (prevTriIndex == -1)
{
// Different edge. Store for later.
EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
if (iter != extraEdgeFaces.end())
{
labelList& eFaces = iter();
label sz = eFaces.size();
eFaces.setSize(sz+1);
eFaces[sz] = triI;
}
else
{
extraEdgeFaces.insert(e, labelList(1, triI));
}
}
else if (edgeFace1[edgeI] == -1)
{
edgeFace1[edgeI] = triI;
}
else
{
edgeFacesRest.insert(edgeI, labelList(1, triI));
//WarningIn("orientSurface(triSurface&)")
// << "Edge " << edgeI << " with centre "
// << mergedCentres[edgeI]
// << " used by more than two triangles: "
// << edgeFace0[edgeI] << ", "
// << edgeFace1[edgeI] << " and " << triI << endl;
Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
if (iter != edgeFacesRest.end())
{
labelList& eFaces = iter();
label sz = eFaces.size();
eFaces.setSize(sz+1);
eFaces[sz] = triI;
}
else
{
edgeFacesRest.insert(edgeI, labelList(1, triI));
}
}
}
}
// Add extraEdgeFaces
edgeI = edgeFace0.size();
edgeFace0.setSize(edgeI + extraEdgeFaces.size());
edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
{
const labelList& eFaces = iter();
// The current edge will become edgeI. Replace all occurrences in
// faceEdges
forAll(eFaces, i)
{
label triI = eFaces[i];
const labelledTri& tri = surf[triI];
FixedList<label, 3>& fEdges = faceEdges[triI];
forAll(tri, fp)
{
edge e(tri[fp], tri[tri.fcIndex(fp)]);
if (e == iter.key())
{
fEdges[fp] = edgeI;
break;
}
}
}
// Add face to edgeFaces
edgeFace0[edgeI] = eFaces[0];
if (eFaces.size() >= 2)
{
edgeFace1[edgeI] = eFaces[1];
if (eFaces.size() > 2)
{
edgeFacesRest.insert
(
edgeI,
SubList<label>(eFaces, eFaces.size()-2, 2)
);
}
}
edgeI++;
}
}
......@@ -1097,6 +1184,24 @@ void Foam::isoSurface::walkOrientation
// nbr points
label nbrFp = findIndex(nbrTri, p0);
if (nbrFp == -1)
{
FatalErrorIn("isoSurface::walkOrientation(..)")
<< "triI:" << triI
<< " tri:" << tri
<< " p0:" << p0
<< " p1:" << p1
<< " fEdges:" << fEdges
<< " edgeI:" << edgeI
<< " edgeFace0:" << edgeFace0[edgeI]
<< " edgeFace1:" << edgeFace1[edgeI]
<< " nbrI:" << nbrI
<< " nbrTri:" << nbrTri
<< abort(FatalError);
}
label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
bool sameOrientation = (p1 == nbrP1);
......@@ -1352,8 +1457,17 @@ Foam::isoSurface::isoSurface
iso_(iso),
mergeDistance_(mergeTol*mesh_.bounds().mag())
{
if (debug)
{
Pout<< "isoSurface :" << nl
<< " isoField : " << cVals.name() << nl
<< " isoValue : " << iso << nl
<< " regularise : " << regularise << nl
<< " mergeTol : " << mergeTol << nl
<< endl;
}
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const labelList& own = mesh_.faceOwner();
// Check
forAll(patches, patchI)
......@@ -1442,24 +1556,6 @@ Foam::isoSurface::isoSurface
snappedCc = -1;
}
// Determine neighbouring snap status
labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
label faceI = pp.start();
forAll(pp, i)
{
neiSnappedCc[faceI-mesh_.nInternalFaces()] =
snappedCc[own[faceI]];
faceI++;
}
}
}
syncTools::swapBoundaryFaceList(mesh_, neiSnappedCc, false);
if (debug)
......
src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
View file @ f0f677e3
......@@ -299,7 +299,16 @@ void Foam::isoSurface::generateTriPoints
)
{
FatalErrorIn("isoSurface::generateTriPoints(..)")
<< "Incorrect size." << abort(FatalError);
<< "Incorrect size." << endl
<< "mesh: nCells:" << mesh_.nCells()
<< " points:" << mesh_.nPoints() << endl
<< "cVals:" << cVals.size() << endl
<< "cCoords:" << cCoords.size() << endl
<< "snappedCc:" << snappedCc.size() << endl
<< "pVals:" << pVals.size() << endl
<< "pCoords:" << pCoords.size() << endl
<< "snappedPoint:" << snappedPoint.size() << endl
<< abort(FatalError);
}
// Determine neighbouring snap status
......
src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
View file @ f0f677e3
......@@ -29,7 +29,6 @@ License
#include "volFields.H"
#include "volPointInterpolation.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -45,8 +44,6 @@ void Foam::sampledIsoSurface::getIsoFields() const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
word pointFldName = "volPointInterpolate(" + isoField_ + ')';
// Get volField
// ~~~~~~~~~~~~
......@@ -54,7 +51,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() : lookup "
Info<< "sampledIsoSurface::getIsoField() : lookup volField "
<< isoField_ << endl;
}
storedVolFieldPtr_.clear();
......@@ -79,7 +76,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() : reading "
Info<< "sampledIsoSurface::getIsoField() : reading volField "
<< isoField_ << " from time " << fvm.time().timeName()
<< endl;
}
......@@ -101,78 +98,212 @@ void Foam::sampledIsoSurface::getIsoFields() const
)
);
volFieldPtr_ = storedVolFieldPtr_.operator->();
}
}
// Interpolate to get pointField
// Get pointField
// ~~~~~~~~~~~~~~
if (!subMeshPtr_.valid())
{
word pointFldName = "volPointInterpolate(" + isoField_ + ')';
if (fvm.foundObject<pointScalarField>(pointFldName))
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() : interpolating "
Info<< "sampledIsoSurface::getIsoField() : lookup pointField "
<< pointFldName << endl;
}
storedPointFieldPtr_.reset
(
volPointInterpolation::New(fvm).interpolate(*volFieldPtr_).ptr()
);
pointFieldPtr_ = storedPointFieldPtr_.operator->();
pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
}
}
else
{
// Not in registry. Interpolate.
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() : checking pointField "
<< pointFldName << " for same time "
<< fvm.time().timeName() << endl;
}
if
(
storedPointFieldPtr_.empty()
|| (fvm.time().timeName() != storedPointFieldPtr_().instance())
)
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() :"
<< " interpolating volField " << volFieldPtr_->name()
<< " to get pointField " << pointFldName << endl;
}
// Get pointField
// ~~~~~~~~~~~~~~
storedPointFieldPtr_.reset
(
volPointInterpolation::New(fvm)
.interpolate(*volFieldPtr_).ptr()
);
storedPointFieldPtr_->checkOut();
pointFieldPtr_ = storedPointFieldPtr_.operator->();
}
}
if (fvm.foundObject<pointScalarField>(pointFldName))
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() : lookup "
<< pointFldName << endl;
Info<< "sampledIsoSurface::getIsoField() : volField "
<< volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
<< " max:" << max(*volFieldPtr_).value() << endl;
Info<< "sampledIsoSurface::getIsoField() : pointField "
<< pointFieldPtr_->name()
<< " min:" << gMin(pointFieldPtr_->internalField())
<< " max:" << gMax(pointFieldPtr_->internalField()) << endl;
}
pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
}
else
{
// Not in registry. Interpolate.
// Get subMesh variants
const fvMesh& subFvm = subMeshPtr_().subMesh();
if (debug)
// Either lookup on the submesh or subset the whole-mesh volField
if (subFvm.foundObject<volScalarField>(isoField_))
{
Info<< "sampledIsoSurface::getIsoField() : checking interpolate "
<< isoField_ << " for same time " << fvm.time().timeName()
<< endl;
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() :"
<< " submesh lookup volField "
<< isoField_ << endl;
}
storedVolSubFieldPtr_.clear();
volSubFieldPtr_ = &subFvm.lookupObject<volScalarField>(isoField_);
}
if
(
storedPointFieldPtr_.empty()
|| (fvm.time().timeName() != storedPointFieldPtr_().instance())
)
else
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() : interpolating "
<< pointFldName << endl;
Info<< "sampledIsoSurface::getIsoField() : subsetting volField "
<< isoField_ << endl;
}
storedVolSubFieldPtr_.reset
(
subMeshPtr_().interpolate
(
*volFieldPtr_
).ptr()
);
storedVolSubFieldPtr_->checkOut();
volSubFieldPtr_ = storedVolSubFieldPtr_.operator->();
}
storedPointFieldPtr_.reset
// Pointfield on submesh
word pointFldName =
"volPointInterpolate("
+ volSubFieldPtr_->name()
+ ')';
if (subFvm.foundObject<pointScalarField>(pointFldName))
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() :"
<< " submesh lookup pointField " << pointFldName << endl;
}
storedPointSubFieldPtr_.clear();
pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
(
pointFldName
);
}
else
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() :"
<< " interpolating submesh volField "
<< volSubFieldPtr_->name()
<< " to get submesh pointField " << pointFldName << endl;
}
storedPointSubFieldPtr_.reset
(
volPointInterpolation::New(fvm).interpolate(*volFieldPtr_).ptr()
volPointInterpolation::New
(
subFvm
).interpolate(*volSubFieldPtr_).ptr()
);
pointFieldPtr_ = storedPointFieldPtr_.operator->();
storedPointSubFieldPtr_->checkOut();
pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
}
if (debug)
{
Info<< "sampledIsoSurface::getIsoField() : volSubField "
<< volSubFieldPtr_->name()
<< " min:" << min(*volSubFieldPtr_).value()
<< " max:" << max(*volSubFieldPtr_).value() << endl;
Info<< "sampledIsoSurface::getIsoField() : pointSubField "
<< pointSubFieldPtr_->name()
<< " min:" << gMin(pointSubFieldPtr_->internalField())
<< " max:" << gMax(pointSubFieldPtr_->internalField()) << endl;
}
}
}