From 6e2e761d37f1d16f09cdc80a5b871d0ef80fb286 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sat, 12 Mar 2016 11:50:04 +0000
Subject: [PATCH] src/sampling/sampledSurface/isoSurface: Correct referencing
to temporary fields
Patch provided by Mattijs Janssens
Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1487
---
.../sampledSurface/isoSurface/isoSurface.C | 691 ++----------------
.../sampledSurface/isoSurface/isoSurface.H | 87 +--
.../isoSurface/isoSurfaceCell.C | 232 ++----
.../isoSurface/isoSurfaceCell.H | 49 +-
.../isoSurface/isoSurfaceCellTemplates.C | 361 ++++-----
.../isoSurface/isoSurfaceTemplates.C | 293 ++++----
.../isoSurface/sampledIsoSurface.C | 72 +-
.../isoSurface/sampledIsoSurface.H | 1 -
8 files changed, 479 insertions(+), 1307 deletions(-)
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index a9d57ad5000..bff64af5cf1 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -37,7 +37,7 @@ License
namespace Foam
{
-defineTypeNameAndDebug(isoSurface, 0);
+ defineTypeNameAndDebug(isoSurface, 0);
}
@@ -58,7 +58,6 @@ bool Foam::isoSurface::noTransform(const tensor& tt) const
}
-// Calculates per face whether couple is collocated.
bool Foam::isoSurface::collocatedPatch(const polyPatch& pp)
{
const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
@@ -66,7 +65,6 @@ bool Foam::isoSurface::collocatedPatch(const polyPatch& pp)
}
-// Calculates per face whether couple is collocated.
Foam::PackedBoolList Foam::isoSurface::collocatedFaces
(
const coupledPolyPatch& pp
@@ -297,7 +295,6 @@ bool Foam::isoSurface::isEdgeOfFaceCut
}
-// Get neighbour value and position.
void Foam::isoSurface::getNeighbour
(
const labelList& boundaryRegion,
@@ -331,7 +328,6 @@ void Foam::isoSurface::getNeighbour
}
-// Determine for every face/cell whether it (possibly) generates triangles.
void Foam::isoSurface::calcCutTypes
(
const labelList& boundaryRegion,
@@ -471,43 +467,6 @@ void Foam::isoSurface::calcCutTypes
}
-// Return the two common points between two triangles
-Foam::labelPair Foam::isoSurface::findCommonPoints
-(
- const labelledTri& tri0,
- const labelledTri& tri1
-)
-{
- labelPair common(-1, -1);
-
- label fp0 = 0;
- label fp1 = findIndex(tri1, tri0[fp0]);
-
- if (fp1 == -1)
- {
- fp0 = 1;
- fp1 = findIndex(tri1, tri0[fp0]);
- }
-
- if (fp1 != -1)
- {
- // So tri0[fp0] is tri1[fp1]
-
- // Find next common point
- label fp0p1 = tri0.fcIndex(fp0);
- label fp1p1 = tri1.fcIndex(fp1);
- label fp1m1 = tri1.rcIndex(fp1);
-
- if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
- {
- common[0] = tri0[fp0];
- common[1] = tri0[fp0p1];
- }
- }
- return common;
-}
-
-
// Caculate centre of surface.
Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
{
@@ -521,75 +480,6 @@ Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
}
-// Replace surface (localPoints, localTris) with single point. Returns
-// point. Destructs arguments.
-Foam::pointIndexHit Foam::isoSurface::collapseSurface
-(
- pointField& localPoints,
- DynamicList<labelledTri, 64>& localTris
-)
-{
- pointIndexHit info(false, vector::zero, localTris.size());
-
- if (localTris.size() == 1)
- {
- const labelledTri& tri = localTris[0];
- info.setPoint(tri.centre(localPoints));
- info.setHit();
- }
- else if (localTris.size() == 2)
- {
- // Check if the two triangles share an edge.
- const labelledTri& tri0 = localTris[0];
- const labelledTri& tri1 = localTris[0];
-
- labelPair shared = findCommonPoints(tri0, tri1);
-
- if (shared[0] != -1)
- {
- info.setPoint
- (
- 0.5
- * (
- tri0.centre(localPoints)
- + tri1.centre(localPoints)
- )
- );
- info.setHit();
- }
- }
- else if (localTris.size())
- {
- // Check if single region. Rare situation.
- triSurface surf
- (
- localTris,
- geometricSurfacePatchList(0),
- localPoints,
- true
- );
- localTris.clearStorage();
-
- labelList faceZone;
- label nZones = surf.markZones
- (
- boolList(surf.nEdges(), false),
- faceZone
- );
-
- if (nZones == 1)
- {
- info.setPoint(calcCentre(surf));
- info.setHit();
- }
- }
-
- return info;
-}
-
-
-// Determine per cell centre whether all the intersections get collapsed
-// to a single point
void Foam::isoSurface::calcSnappedCc
(
const labelList& boundaryRegion,
@@ -755,8 +645,6 @@ void Foam::isoSurface::calcSnappedCc
}
-// Determine per meshpoint whether all the intersections get collapsed
-// to a single point
void Foam::isoSurface::calcSnappedPoint
(
const PackedBoolList& isBoundaryPoint,
@@ -1131,7 +1019,6 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
}
-// Does face use valid vertices?
bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
{
// Simple check on indices ok.
@@ -1202,435 +1089,6 @@ bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
}
-void Foam::isoSurface::calcAddressing
-(
- const triSurface& surf,
- List<FixedList<label, 3>>& faceEdges,
- labelList& edgeFace0,
- labelList& edgeFace1,
- Map<labelList>& edgeFacesRest
-) const
-{
- const pointField& points = surf.points();
-
- pointField edgeCentres(3*surf.size());
- label edgeI = 0;
- forAll(surf, triI)
- {
- const labelledTri& tri = surf[triI];
- edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
- edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
- edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
- }
-
- pointField mergedCentres;
- labelList oldToMerged;
- bool hasMerged = mergePoints
- (
- edgeCentres,
- mergeDistance_,
- false,
- oldToMerged,
- mergedCentres
- );
-
- if (debug)
- {
- Pout<< "isoSurface : detected "
- << mergedCentres.size()
- << " geometric edges on " << surf.size() << " triangles." << endl;
- }
-
- if (!hasMerged)
- {
- return;
- }
-
-
- // Determine faceEdges
- faceEdges.setSize(surf.size());
- edgeI = 0;
- forAll(surf, triI)
- {
- faceEdges[triI][0] = oldToMerged[edgeI++];
- faceEdges[triI][1] = oldToMerged[edgeI++];
- faceEdges[triI][2] = oldToMerged[edgeI++];
- }
-
-
- // Determine edgeFaces
- edgeFace0.setSize(mergedCentres.size());
- edgeFace0 = -1;
- edgeFace1.setSize(mergedCentres.size());
- 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;
- label edgeI = oldToMerged[oldEdgeI];
-
- if (edgeFace0[edgeI] == -1)
- {
- // First triangle for edge
- edgeFace0[edgeI] = triI;
- }
- else
- {
- //- 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)
- {
- 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
- {
- //WarningInFunction
- // << "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++;
- }
-}
-
-
-void Foam::isoSurface::walkOrientation
-(
- const triSurface& surf,
- const List<FixedList<label, 3>>& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const label seedTriI,
- labelList& flipState
-)
-{
- // Do walk for consistent orientation.
- DynamicList<label> changedFaces(surf.size());
-
- changedFaces.append(seedTriI);
-
- while (changedFaces.size())
- {
- DynamicList<label> newChangedFaces(changedFaces.size());
-
- forAll(changedFaces, i)
- {
- label triI = changedFaces[i];
- const labelledTri& tri = surf[triI];
- const FixedList<label, 3>& fEdges = faceEdges[triI];
-
- forAll(fEdges, fp)
- {
- label edgeI = fEdges[fp];
-
- // my points:
- label p0 = tri[fp];
- label p1 = tri[tri.fcIndex(fp)];
-
- label nbrI =
- (
- edgeFace0[edgeI] != triI
- ? edgeFace0[edgeI]
- : edgeFace1[edgeI]
- );
-
- if (nbrI != -1 && flipState[nbrI] == -1)
- {
- const labelledTri& nbrTri = surf[nbrI];
-
- // nbr points
- label nbrFp = findIndex(nbrTri, p0);
-
- if (nbrFp == -1)
- {
- FatalErrorInFunction
- << "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);
-
- if (flipState[triI] == 0)
- {
- flipState[nbrI] = (sameOrientation ? 0 : 1);
- }
- else
- {
- flipState[nbrI] = (sameOrientation ? 1 : 0);
- }
- newChangedFaces.append(nbrI);
- }
- }
- }
-
- changedFaces.transfer(newChangedFaces);
- }
-}
-
-
-void Foam::isoSurface::orientSurface
-(
- triSurface& surf,
- const List<FixedList<label, 3>>& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const Map<labelList>& edgeFacesRest
-)
-{
- // -1 : unvisited
- // 0 : leave as is
- // 1 : flip
- labelList flipState(surf.size(), -1);
-
- label seedTriI = 0;
-
- while (true)
- {
- // Find first unvisited triangle
- for
- (
- ;
- seedTriI < surf.size() && flipState[seedTriI] != -1;
- seedTriI++
- )
- {}
-
- if (seedTriI == surf.size())
- {
- break;
- }
-
- // Note: Determine orientation of seedTriI?
- // for now assume it is ok
- flipState[seedTriI] = 0;
-
- walkOrientation
- (
- surf,
- faceEdges,
- edgeFace0,
- edgeFace1,
- seedTriI,
- flipState
- );
- }
-
- // Do actual flipping
- surf.clearOut();
- forAll(surf, triI)
- {
- if (flipState[triI] == 1)
- {
- labelledTri tri(surf[triI]);
-
- surf[triI][0] = tri[0];
- surf[triI][1] = tri[2];
- surf[triI][2] = tri[1];
- }
- else if (flipState[triI] == -1)
- {
- FatalErrorInFunction
- << "problem" << abort(FatalError);
- }
- }
-}
-
-
-// Checks if triangle is connected through edgeI only.
-bool Foam::isoSurface::danglingTriangle
-(
- const FixedList<label, 3>& fEdges,
- const labelList& edgeFace1
-)
-{
- label nOpen = 0;
- forAll(fEdges, i)
- {
- if (edgeFace1[fEdges[i]] == -1)
- {
- nOpen++;
- }
- }
-
- if (nOpen == 1 || nOpen == 2 || nOpen == 3)
- {
- return true;
- }
- else
- {
- return false;
- }
-}
-
-
-// Mark triangles to keep. Returns number of dangling triangles.
-Foam::label Foam::isoSurface::markDanglingTriangles
-(
- const List<FixedList<label, 3>>& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const Map<labelList>& edgeFacesRest,
- boolList& keepTriangles
-)
-{
- keepTriangles.setSize(faceEdges.size());
- keepTriangles = true;
-
- label nDangling = 0;
-
- // Remove any dangling triangles
- forAllConstIter(Map<labelList>, edgeFacesRest, iter)
- {
- // These are all the non-manifold edges. Filter out all triangles
- // with only one connected edge (= this edge)
-
- label edgeI = iter.key();
- const labelList& otherEdgeFaces = iter();
-
- // Remove all dangling triangles
- if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
- {
- keepTriangles[edgeFace0[edgeI]] = false;
- nDangling++;
- }
- if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
- {
- keepTriangles[edgeFace1[edgeI]] = false;
- nDangling++;
- }
- forAll(otherEdgeFaces, i)
- {
- label triI = otherEdgeFaces[i];
- if (danglingTriangle(faceEdges[triI], edgeFace1))
- {
- keepTriangles[triI] = false;
- nDangling++;
- }
- }
- }
- return nDangling;
-}
-
-
Foam::triSurface Foam::isoSurface::subsetMesh
(
const triSurface& s,
@@ -1953,57 +1411,58 @@ Foam::isoSurface::isoSurface
}
+ {
+ DynamicList<point> triPoints(3*nCutCells_);
+ DynamicList<label> triMeshCells(nCutCells_);
- DynamicList<point> triPoints(nCutCells_);
- DynamicList<label> triMeshCells(nCutCells_);
-
- generateTriPoints
- (
- cValsPtr_(),
- pVals_,
+ generateTriPoints
+ (
+ cValsPtr_(),
+ pVals_,
- meshC,
- mesh_.points(),
+ meshC,
+ mesh_.points(),
- snappedPoints,
- snappedCc,
- snappedPoint,
+ snappedPoints,
+ snappedCc,
+ snappedPoint,
- triPoints,
- triMeshCells
- );
+ triPoints, // 3 points of the triangle
+ triMeshCells // per triangle the originating cell
+ );
- if (debug)
- {
- Pout<< "isoSurface : generated " << triMeshCells.size()
- << " unmerged triangles from " << triPoints.size()
- << " unmerged points." << endl;
- }
+ if (debug)
+ {
+ Pout<< "isoSurface : generated " << triMeshCells.size()
+ << " unmerged triangles from " << triPoints.size()
+ << " unmerged points." << endl;
+ }
- // Merge points and compact out non-valid triangles
- labelList triMap; // merged to unmerged triangle
- triSurface::operator=
- (
- stitchTriPoints
+ // Merge points and compact out non-valid triangles
+ labelList triMap; // merged to unmerged triangle
+ triSurface::operator=
(
- true, // check for duplicate tris
- triPoints,
- triPointMergeMap_, // unmerged to merged point
- triMap
- )
- );
+ stitchTriPoints
+ (
+ true, // check for duplicate tris
+ triPoints,
+ triPointMergeMap_, // unmerged to merged point
+ triMap
+ )
+ );
- if (debug)
- {
- Pout<< "isoSurface : generated " << triMap.size()
- << " merged triangles." << endl;
- }
+ if (debug)
+ {
+ Pout<< "isoSurface : generated " << triMap.size()
+ << " merged triangles." << endl;
+ }
- meshCells_.setSize(triMap.size());
- forAll(triMap, i)
- {
- meshCells_[i] = triMeshCells[triMap[i]];
+ meshCells_.setSize(triMap.size());
+ forAll(triMap, i)
+ {
+ meshCells_[i] = triMeshCells[triMap[i]];
+ }
}
if (debug)
@@ -2019,70 +1478,6 @@ Foam::isoSurface::isoSurface
}
- if (false)
- {
- List<FixedList<label, 3>> faceEdges;
- labelList edgeFace0, edgeFace1;
- Map<labelList> edgeFacesRest;
-
-
- while (true)
- {
- // Calculate addressing
- calcAddressing
- (
- *this,
- faceEdges,
- edgeFace0,
- edgeFace1,
- edgeFacesRest
- );
-
- // See if any dangling triangles
- boolList keepTriangles;
- label nDangling = markDanglingTriangles
- (
- faceEdges,
- edgeFace0,
- edgeFace1,
- edgeFacesRest,
- keepTriangles
- );
-
- if (debug)
- {
- Pout<< "isoSurface : detected " << nDangling
- << " dangling triangles." << endl;
- }
-
- if (nDangling == 0)
- {
- break;
- }
-
- // Create face map (new to old)
- labelList subsetTriMap(findIndices(keepTriangles, true));
-
- labelList subsetPointMap;
- labelList reversePointMap;
- triSurface::operator=
- (
- subsetMesh
- (
- *this,
- subsetTriMap,
- reversePointMap,
- subsetPointMap
- )
- );
- meshCells_ = labelField(meshCells_, subsetTriMap);
- inplaceRenumber(reversePointMap, triPointMergeMap_);
- }
-
- orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
- }
-
-
if (debug)
{
fileName stlFile = mesh_.time().path() + ".stl";
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 967536b7c6d..1493c2ccf94 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -173,6 +173,7 @@ class isoSurface
const bool neiLower
) const;
+ //- Get neighbour value and position.
void getNeighbour
(
const labelList& boundaryRegion,
@@ -184,7 +185,8 @@ class isoSurface
point& nbrPoint
) const;
- //- Set faceCutType,cellCutType.
+ //- Determine for every face/cell whether it (possibly) generates
+ // triangles.
void calcCutTypes
(
const labelList& boundaryRegion,
@@ -193,20 +195,8 @@ class isoSurface
const scalarField& pVals
);
- static labelPair findCommonPoints
- (
- const labelledTri&,
- const labelledTri&
- );
-
static point calcCentre(const triSurface&);
- static pointIndexHit collapseSurface
- (
- pointField& localPoints,
- DynamicList<labelledTri, 64>& localTris
- );
-
//- Determine per cc whether all near cuts can be snapped to single
// point.
void calcSnappedCc
@@ -257,6 +247,9 @@ class isoSurface
const Type& snapP1
) const;
+
+ //- Note: cannot use simpler isoSurfaceCell::generateTriPoints since
+ // the need here to sometimes pass in remote 'snappedPoints'
template<class Type>
void generateTriPoints
(
@@ -323,6 +316,14 @@ class isoSurface
DynamicList<label>& triMeshCells
) const;
+ template<class Type>
+ static tmp<Field<Type>> interpolate
+ (
+ const label nPoints,
+ const labelList& triPointMergeMap,
+ const DynamicList<Type>& unmergedValues
+ );
+
triSurface stitchTriPoints
(
const bool checkDuplicates,
@@ -334,54 +335,6 @@ class isoSurface
//- Check single triangle for (topological) validity
static bool validTri(const triSurface&, const label);
- //- Determine edge-face addressing
- void calcAddressing
- (
- const triSurface& surf,
- List<FixedList<label, 3>>& faceEdges,
- labelList& edgeFace0,
- labelList& edgeFace1,
- Map<labelList>& edgeFacesRest
- ) const;
-
- //- Determine orientation
- static void walkOrientation
- (
- const triSurface& surf,
- const List<FixedList<label, 3>>& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const label seedTriI,
- labelList& flipState
- );
-
- //- Orient surface
- static void orientSurface
- (
- triSurface&,
- const List<FixedList<label, 3>>& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const Map<labelList>& edgeFacesRest
- );
-
- //- Is triangle (given by 3 edges) not fully connected?
- static bool danglingTriangle
- (
- const FixedList<label, 3>& fEdges,
- const labelList& edgeFace1
- );
-
- //- Mark all non-fully connected triangles
- static label markDanglingTriangles
- (
- const List<FixedList<label, 3>>& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const Map<labelList>& edgeFacesRest,
- boolList& keepTriangles
- );
-
static triSurface subsetMesh
(
const triSurface& s,
@@ -392,6 +345,10 @@ class isoSurface
public:
+ //- Declare friendship with isoSurfaceCell to share some functionality
+ friend class isoSurfaceCell;
+
+
//- Runtime type information
TypeName("isoSurface");
@@ -413,18 +370,12 @@ public:
// Member Functions
- //- For every face original cell in mesh
+ //- For every triangle the original cell in mesh
const labelList& meshCells() const
{
return meshCells_;
}
- //- For every unmerged triangle point the point in the triSurface
- const labelList& triPointMergeMap() const
- {
- return triPointMergeMap_;
- }
-
//- Interpolates cCoords,pCoords. Uses the references to the original
// fields used to create the iso surface.
template<class Type>
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
index c2fec197418..afd756ac3ce 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
@@ -35,7 +35,7 @@ License
namespace Foam
{
-defineTypeNameAndDebug(isoSurfaceCell, 0);
+ defineTypeNameAndDebug(isoSurfaceCell, 0);
}
@@ -212,8 +212,6 @@ void Foam::isoSurfaceCell::calcCutTypes
}
-
-// Return the two common points between two triangles
Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
(
const labelledTri& tri0,
@@ -250,7 +248,6 @@ Foam::labelPair Foam::isoSurfaceCell::findCommonPoints
}
-// Caculate centre of surface.
Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
{
vector sum = vector::zero;
@@ -263,8 +260,6 @@ Foam::point Foam::isoSurfaceCell::calcCentre(const triSurface& s)
}
-// Replace surface (localPoints, localTris) with single point. Returns
-// point. Destructs arguments.
Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface
(
const label cellI,
@@ -532,7 +527,6 @@ void Foam::isoSurfaceCell::calcSnappedCc
}
-// Generate triangles for face connected to pointI
void Foam::isoSurfaceCell::genPointTris
(
const scalarField& cellValues,
@@ -605,7 +599,6 @@ void Foam::isoSurfaceCell::genPointTris
}
-// Generate triangle for tet connected to pointI
void Foam::isoSurfaceCell::genPointTris
(
const scalarField& pointValues,
@@ -1053,7 +1046,6 @@ Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
}
-// Does face use valid vertices?
bool Foam::isoSurfaceCell::validTri(const triSurface& surf, const label faceI)
{
// Simple check on indices ok.
@@ -1222,143 +1214,6 @@ void Foam::isoSurfaceCell::calcAddressing
}
-//void Foam::isoSurfaceCell::walkOrientation
-//(
-// const triSurface& surf,
-// const List<FixedList<label, 3>>& faceEdges,
-// const labelList& edgeFace0,
-// const labelList& edgeFace1,
-// const label seedTriI,
-// labelList& flipState
-//)
-//{
-// // Do walk for consistent orientation.
-// DynamicList<label> changedFaces(surf.size());
-//
-// changedFaces.append(seedTriI);
-//
-// while (changedFaces.size())
-// {
-// DynamicList<label> newChangedFaces(changedFaces.size());
-//
-// forAll(changedFaces, i)
-// {
-// label triI = changedFaces[i];
-// const labelledTri& tri = surf[triI];
-// const FixedList<label, 3>& fEdges = faceEdges[triI];
-//
-// forAll(fEdges, fp)
-// {
-// label edgeI = fEdges[fp];
-//
-// // my points:
-// label p0 = tri[fp];
-// label p1 = tri[tri.fcIndex(fp)];
-//
-// label nbrI =
-// (
-// edgeFace0[edgeI] != triI
-// ? edgeFace0[edgeI]
-// : edgeFace1[edgeI]
-// );
-//
-// if (nbrI != -1 && flipState[nbrI] == -1)
-// {
-// const labelledTri& nbrTri = surf[nbrI];
-//
-// // nbr points
-// label nbrFp = findIndex(nbrTri, p0);
-// label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-//
-// bool sameOrientation = (p1 == nbrP1);
-//
-// if (flipState[triI] == 0)
-// {
-// flipState[nbrI] = (sameOrientation ? 0 : 1);
-// }
-// else
-// {
-// flipState[nbrI] = (sameOrientation ? 1 : 0);
-// }
-// newChangedFaces.append(nbrI);
-// }
-// }
-// }
-//
-// changedFaces.transfer(newChangedFaces);
-// }
-//}
-//
-//
-//void Foam::isoSurfaceCell::orientSurface
-//(
-// triSurface& surf,
-// const List<FixedList<label, 3>>& faceEdges,
-// const labelList& edgeFace0,
-// const labelList& edgeFace1,
-// const Map<labelList>& edgeFacesRest
-//)
-//{
-// // -1 : unvisited
-// // 0 : leave as is
-// // 1 : flip
-// labelList flipState(surf.size(), -1);
-//
-// label seedTriI = 0;
-//
-// while (true)
-// {
-// // Find first unvisited triangle
-// for
-// (
-// ;
-// seedTriI < surf.size() && flipState[seedTriI] != -1;
-// seedTriI++
-// )
-// {}
-//
-// if (seedTriI == surf.size())
-// {
-// break;
-// }
-//
-// // Note: Determine orientation of seedTriI?
-// // for now assume it is ok
-// flipState[seedTriI] = 0;
-//
-// walkOrientation
-// (
-// surf,
-// faceEdges,
-// edgeFace0,
-// edgeFace1,
-// seedTriI,
-// flipState
-// );
-// }
-//
-// // Do actual flipping
-// surf.clearOut();
-// forAll(surf, triI)
-// {
-// if (flipState[triI] == 1)
-// {
-// labelledTri tri(surf[triI]);
-//
-// surf[triI][0] = tri[0];
-// surf[triI][1] = tri[2];
-// surf[triI][2] = tri[1];
-// }
-// else if (flipState[triI] == -1)
-// {
-// FatalErrorInFunction
-// << "problem" << abort(FatalError);
-// }
-// }
-//}
-
-
-// Checks if triangle is connected through edgeI only.
bool Foam::isoSurfaceCell::danglingTriangle
(
const FixedList<label, 3>& fEdges,
@@ -1385,7 +1240,6 @@ bool Foam::isoSurfaceCell::danglingTriangle
}
-// Mark triangles to keep. Returns number of dangling triangles.
Foam::label Foam::isoSurfaceCell::markDanglingTriangles
(
const List<FixedList<label, 3>>& faceEdges,
@@ -1605,57 +1459,59 @@ Foam::isoSurfaceCell::isoSurfaceCell
}
+ {
+ DynamicList<point> triPoints(nCutCells_);
+ DynamicList<label> triMeshCells(nCutCells_);
- DynamicList<point> triPoints(nCutCells_);
- DynamicList<label> triMeshCells(nCutCells_);
-
- generateTriPoints
- (
- cVals,
- pVals,
+ generateTriPoints
+ (
+ cVals,
+ pVals,
- mesh_.cellCentres(),
- mesh_.points(),
+ mesh_.cellCentres(),
+ mesh_.points(),
- snappedPoints,
- snappedCc,
- snappedPoint,
+ snappedPoints,
+ snappedCc,
+ snappedPoint,
- triPoints,
- triMeshCells
- );
+ triPoints,
+ triMeshCells
+ );
- if (debug)
- {
- Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
- << " unmerged triangles." << endl;
- }
+ if (debug)
+ {
+ Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
+ << " unmerged triangles." << endl;
+ }
- // Merge points and compact out non-valid triangles
- labelList triMap; // merged to unmerged triangle
- triSurface::operator=
- (
- stitchTriPoints
+ // Merge points and compact out non-valid triangles
+ labelList triMap;
+ triSurface::operator=
(
- regularise, // check for duplicate tris
- triPoints,
- triPointMergeMap_, // unmerged to merged point
- triMap
- )
- );
+ stitchTriPoints
+ (
+ regularise, // check for duplicate tris
+ triPoints,
+ triPointMergeMap_, // unmerged to merged point
+ triMap // merged to unmerged triangle
+ )
+ );
- if (debug)
- {
- Pout<< "isoSurfaceCell : generated " << triMap.size()
- << " merged triangles." << endl;
- }
+ if (debug)
+ {
+ Pout<< "isoSurfaceCell : generated " << triMap.size()
+ << " merged triangles." << endl;
+ }
- meshCells_.setSize(triMap.size());
- forAll(triMap, i)
- {
- meshCells_[i] = triMeshCells[triMap[i]];
+ meshCells_.setSize(triMap.size());
+ forAll(triMap, i)
+ {
+ meshCells_[i] = triMeshCells[triMap[i]];
+ }
}
+
if (debug)
{
Pout<< "isoSurfaceCell : checking " << size()
@@ -1728,8 +1584,6 @@ Foam::isoSurfaceCell::isoSurfaceCell
meshCells_ = labelField(meshCells_, subsetTriMap);
inplaceRenumber(reversePointMap, triPointMergeMap_);
}
-
- //orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
}
}
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
index 034fbd31712..91de7fc8e3d 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
@@ -139,14 +139,18 @@ class isoSurfaceCell
const scalarField& pointValues
);
+ //- Return the two common points between two triangles
static labelPair findCommonPoints
(
const labelledTri&,
const labelledTri&
);
+ //- Caculate centre of surface.
static point calcCentre(const triSurface&);
+ //- Replace surface (localPoints, localTris) with single point.
+ // Returns point. Destroys arguments.
pointIndexHit collapseSurface
(
const label cellI,
@@ -214,19 +218,15 @@ class isoSurfaceCell
void generateTriPoints
(
const DynamicList<Type>& snapped,
- const scalar isoVal0,
const scalar s0,
const Type& p0,
const label p0Index,
- const scalar isoVal1,
const scalar s1,
const Type& p1,
const label p1Index,
- const scalar isoVal2,
const scalar s2,
const Type& p2,
const label p2Index,
- const scalar isoVal3,
const scalar s3,
const Type& p3,
const label p3Index,
@@ -271,27 +271,6 @@ class isoSurfaceCell
Map<labelList>& edgeFacesRest
) const;
- ////- Determine orientation
- //static void walkOrientation
- //(
- // const triSurface& surf,
- // const List<FixedList<label, 3>>& faceEdges,
- // const labelList& edgeFace0,
- // const labelList& edgeFace1,
- // const label seedTriI,
- // labelList& flipState
- //);
-
- ////- Orient surface
- //static void orientSurface
- //(
- // triSurface&,
- // const List<FixedList<label, 3>>& faceEdges,
- // const labelList& edgeFace0,
- // const labelList& edgeFace1,
- // const Map<labelList>& edgeFacesRest
- //);
-
//- Is triangle (given by 3 edges) not fully connected?
static bool danglingTriangle
(
@@ -317,8 +296,6 @@ class isoSurfaceCell
labelList& newToOldPoints
);
- //- Combine all triangles inside a cell into a minimal triangulation
- void combineCellTriangles();
public:
@@ -348,24 +325,6 @@ public:
return meshCells_;
}
- //- For every unmerged triangle point the point in the triSurface
- const labelList triPointMergeMap() const
- {
- return triPointMergeMap_;
- }
-
-
- //- Interpolates cCoords,pCoords. Takes the original fields
- // used to create the iso surface.
- template<class Type>
- tmp<Field<Type>> interpolate
- (
- const scalarField& cVals,
- const scalarField& pVals,
- const Field<Type>& cCoords,
- const Field<Type>& pCoords
- ) const;
-
//- Interpolates cCoords,pCoords.
template<class Type>
tmp<Field<Type>> interpolate
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index 1eba2063c40..ea26461e5ef 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -26,6 +26,7 @@ License
#include "isoSurfaceCell.H"
#include "polyMesh.H"
#include "tetMatcher.H"
+#include "isoSurface.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -76,22 +77,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
const DynamicList<Type>& snapped,
- const scalar isoVal0,
const scalar s0,
const Type& p0,
const label p0Index,
- const scalar isoVal1,
const scalar s1,
const Type& p1,
const label p1Index,
- const scalar isoVal2,
const scalar s2,
const Type& p2,
const label p2Index,
- const scalar isoVal3,
const scalar s3,
const Type& p3,
const label p3Index,
@@ -117,167 +114,203 @@ void Foam::isoSurfaceCell::generateTriPoints
triIndex |= 8;
}
- // Form the vertices of the triangles for each case
+ /* Form the vertices of the triangles for each case */
switch (triIndex)
{
case 0x00:
case 0x0F:
break;
- case 0x0E:
case 0x01:
+ case 0x0E:
{
- // 0 is common point. Orient such that normal points in positive
- // gradient direction
- if (isoVal0 >= isoVal1)
- {
- pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
- pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- }
- else
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+ );
+ if (triIndex == 0x0E)
{
- pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
- pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x0D:
case 0x02:
+ case 0x0D:
{
- // 1 is common point
- if (isoVal1 >= isoVal0)
- {
- pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- }
- else
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+ );
+
+ if (triIndex == 0x0D)
{
- pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x0C:
case 0x03:
+ case 0x0C:
{
- Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
- Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
-
- if (isoVal0 >= isoVal3)
- {
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- pts.append(s02);
- pts.append(s13);
- pts.append(s13);
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- pts.append(s02);
- }
- else
+ Type p0p2 =
+ generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+ Type p1p3 =
+ generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+ );
+ pts.append(p1p3);
+ pts.append(p0p2);
+
+ pts.append(p1p3);
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+ );
+ pts.append(p0p2);
+
+ if (triIndex == 0x0C)
{
- pts.append(s02);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- pts.append(s13);
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- pts.append(s13);
- pts.append(s02);
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-5], pts[sz-4]);
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x0B:
case 0x04:
+ case 0x0B:
{
- // 2 is common point
- if (isoVal2 >= isoVal0)
+ pts.append
+ (
+ generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)
+ );
+
+ if (triIndex == 0x0B)
{
- pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
- pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
- }
- else
- {
- pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
- pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x0A:
case 0x05:
+ case 0x0A:
{
- Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
- Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
- if (isoVal3 >= isoVal0)
- {
- pts.append(s01);
- pts.append(s23);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- pts.append(s01);
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- pts.append(s23);
- }
- else
+ Type p0p1 =
+ generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+ Type p2p3 =
+ generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+ pts.append(p0p1);
+ pts.append(p2p3);
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+ );
+
+ pts.append(p0p1);
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+ );
+ pts.append(p2p3);
+
+ if (triIndex == 0x0A)
{
- pts.append(s23);
- pts.append(s01);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- pts.append(s01);
- pts.append(s23);
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-5], pts[sz-4]);
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x09:
case 0x06:
+ case 0x09:
{
- Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
- Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
- if (isoVal3 >= isoVal1)
- {
- pts.append(s01);
- pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
- pts.append(s23);
- pts.append(s01);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
- pts.append(s23);
- }
- else
+ Type p0p1 =
+ generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+ Type p2p3 =
+ generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+ pts.append(p0p1);
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+ );
+ pts.append(p2p3);
+
+ pts.append(p0p1);
+ pts.append(p2p3);
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+ );
+
+ if (triIndex == 0x09)
{
- pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
- pts.append(s01);
- pts.append(s23);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
- pts.append(s01);
- pts.append(s23);
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-5], pts[sz-4]);
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x07:
case 0x08:
+ case 0x07:
{
- // 3 is common point
- if (isoVal3 >= isoVal0)
+ pts.append
+ (
+ generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)
+ );
+ if (triIndex == 0x07)
{
- pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
- pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
- }
- else
- {
- pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
- pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
@@ -341,22 +374,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
snappedPoints,
- pVals_[f0[1]],
pVals[f0[1]],
pCoords[f0[1]],
snappedPoint[f0[1]],
- pVals_[f0[0]],
pVals[f0[0]],
pCoords[f0[0]],
snappedPoint[f0[0]],
- pVals_[f0[2]],
pVals[f0[2]],
pCoords[f0[2]],
snappedPoint[f0[2]],
- pVals_[oppositeI],
pVals[oppositeI],
pCoords[oppositeI],
snappedPoint[oppositeI],
@@ -370,22 +399,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
snappedPoints,
- pVals_[f0[0]],
pVals[f0[0]],
pCoords[f0[0]],
snappedPoint[f0[0]],
- pVals_[f0[1]],
pVals[f0[1]],
pCoords[f0[1]],
snappedPoint[f0[1]],
- pVals_[f0[2]],
pVals[f0[2]],
pCoords[f0[2]],
snappedPoint[f0[2]],
- pVals_[oppositeI],
pVals[oppositeI],
pCoords[oppositeI],
snappedPoint[oppositeI],
@@ -411,7 +436,6 @@ void Foam::isoSurfaceCell::generateTriPoints
}
label fp = f.fcIndex(fp0);
-
for (label i = 2; i < f.size(); i++)
{
label nextFp = f.fcIndex(fp);
@@ -425,22 +449,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
snappedPoints,
- pVals_[tri[1]],
pVals[tri[1]],
pCoords[tri[1]],
snappedPoint[tri[1]],
- pVals_[tri[0]],
pVals[tri[0]],
pCoords[tri[0]],
snappedPoint[tri[0]],
- pVals_[tri[2]],
pVals[tri[2]],
pCoords[tri[2]],
snappedPoint[tri[2]],
- cVals_[cellI],
cVals[cellI],
cCoords[cellI],
snappedCc[cellI],
@@ -454,22 +474,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
snappedPoints,
- pVals_[tri[0]],
pVals[tri[0]],
pCoords[tri[0]],
snappedPoint[tri[0]],
- pVals_[tri[1]],
pVals[tri[1]],
pCoords[tri[1]],
snappedPoint[tri[1]],
- pVals_[tri[2]],
pVals[tri[2]],
pCoords[tri[2]],
snappedPoint[tri[2]],
- cVals_[cellI],
cVals[cellI],
cCoords[cellI],
snappedCc[cellI],
@@ -495,7 +511,7 @@ void Foam::isoSurfaceCell::generateTriPoints
if (countNotFoundTets > 0)
{
- WarningInFunction
+ WarningIn("Foam::isoSurfaceCell::generateTriPoints(..)")
<< "Could not find " << countNotFoundTets
<< " tet base points, which may lead to inverted triangles."
<< endl;
@@ -507,68 +523,14 @@ void Foam::isoSurfaceCell::generateTriPoints
template<class Type>
-Foam::tmp<Foam::Field<Type>>
-Foam::isoSurfaceCell::interpolate
-(
- const scalarField& cVals,
- const scalarField& pVals,
- const Field<Type>& cCoords,
- const Field<Type>& pCoords
-) const
-{
- DynamicList<Type> triPoints(nCutCells_);
- DynamicList<label> triMeshCells(nCutCells_);
-
- // Dummy snap data
- DynamicList<Type> snappedPoints;
- labelList snappedCc(mesh_.nCells(), -1);
- labelList snappedPoint(mesh_.nPoints(), -1);
-
-
- generateTriPoints
- (
- cVals,
- pVals,
-
- cCoords,
- pCoords,
-
- snappedPoints,
- snappedCc,
- snappedPoint,
-
- triPoints,
- triMeshCells
- );
-
-
- // One value per point
- tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
- Field<Type>& values = tvalues();
-
- forAll(triPoints, i)
- {
- label mergedPointI = triPointMergeMap_[i];
-
- if (mergedPointI >= 0)
- {
- values[mergedPointI] = triPoints[i];
- }
- }
-
- return tvalues;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type>>
+Foam::tmp<Foam::Field<Type> >
Foam::isoSurfaceCell::interpolate
(
const Field<Type>& cCoords,
const Field<Type>& pCoords
) const
{
- DynamicList<Type> triPoints(nCutCells_);
+ DynamicList<Type> triPoints(3*nCutCells_);
DynamicList<label> triMeshCells(nCutCells_);
// Dummy snap data
@@ -576,7 +538,6 @@ Foam::isoSurfaceCell::interpolate
labelList snappedCc(mesh_.nCells(), -1);
labelList snappedPoint(mesh_.nPoints(), -1);
-
generateTriPoints
(
cVals_,
@@ -593,22 +554,12 @@ Foam::isoSurfaceCell::interpolate
triMeshCells
);
-
- // One value per point
- tmp<Field<Type>> tvalues(new Field<Type>(points().size()));
- Field<Type>& values = tvalues.ref();
-
- forAll(triPoints, i)
- {
- label mergedPointI = triPointMergeMap_[i];
-
- if (mergedPointI >= 0)
- {
- values[mergedPointI] = triPoints[i];
- }
- }
-
- return tvalues;
+ return isoSurface::interpolate
+ (
+ points().size(),
+ triPointMergeMap_,
+ triPoints
+ );
}
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index ea92b65ef74..c0444401db1 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -120,7 +120,7 @@ Foam::isoSurface::adaptPatchFields
{
fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
(
- fld.boundaryField()[patchI]
+ sliceFld.boundaryField()[patchI]
);
const scalarField& w = mesh.weights().boundaryField()[patchI];
@@ -240,8 +240,8 @@ void Foam::isoSurface::generateTriPoints
case 0x0F:
break;
- case 0x0E:
case 0x01:
+ case 0x0E:
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
@@ -254,10 +254,16 @@ void Foam::isoSurface::generateTriPoints
(
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
);
+ if (triIndex == 0x0E)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x0D:
case 0x02:
+ case 0x0D:
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
@@ -270,97 +276,133 @@ void Foam::isoSurface::generateTriPoints
(
generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
);
+ if (triIndex == 0x0D)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x0C:
case 0x03:
- {
- Type tp1 =
- generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
- Type tp2 =
- generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
+ case 0x0C:
+ {
+ Type p0p2 =
+ generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
+ Type p1p3 =
+ generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
- points.append
- (
- generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
- );
- points.append(tp1);
- points.append(tp2);
- points.append(tp2);
- points.append
- (
- generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
- );
- points.append(tp1);
- }
+ points.append
+ (
+ generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+ );
+ points.append(p1p3);
+ points.append(p0p2);
+
+ points.append(p1p3);
+ points.append
+ (
+ generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+ );
+ points.append(p0p2);
+ }
+ if (triIndex == 0x0C)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-5], points[sz-4]);
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x0B:
case 0x04:
- {
- points.append
- (
- generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
- );
- points.append
- (
- generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
- );
- points.append
- (
- generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
- );
- }
+ case 0x0B:
+ {
+ points.append
+ (
+ generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
+ );
+ points.append
+ (
+ generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
+ );
+ points.append
+ (
+ generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
+ );
+ }
+ if (triIndex == 0x0B)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x0A:
case 0x05:
- {
- Type tp0 =
- generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
- Type tp1 =
- generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+ case 0x0A:
+ {
+ Type p0p1 =
+ generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+ Type p2p3 =
+ generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+
+ points.append(p0p1);
+ points.append(p2p3);
+ points.append
+ (
+ generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+ );
- points.append(tp0);
- points.append(tp1);
- points.append
- (
- generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
- );
- points.append(tp0);
- points.append
- (
- generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
- );
- points.append(tp1);
- }
+ points.append(p0p1);
+ points.append
+ (
+ generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+ );
+ points.append(p2p3);
+ }
+ if (triIndex == 0x0A)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-5], points[sz-4]);
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x09:
case 0x06:
- {
- Type tp0 =
- generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
- Type tp1 =
- generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+ case 0x09:
+ {
+ Type p0p1 =
+ generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+ Type p2p3 =
+ generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
- points.append(tp0);
- points.append
- (
- generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
- );
- points.append(tp1);
- points.append(tp0);
- points.append
- (
- generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
- );
- points.append(tp1);
- }
+ points.append(p0p1);
+ points.append
+ (
+ generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
+ );
+ points.append(p2p3);
+
+ points.append(p0p1);
+ points.append(p2p3);
+ points.append
+ (
+ generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
+ );
+ }
+ if (triIndex == 0x09)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-5], points[sz-4]);
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x07:
case 0x08:
+ case 0x07:
points.append
(
generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
@@ -373,6 +415,12 @@ void Foam::isoSurface::generateTriPoints
(
generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
);
+ if (triIndex == 0x07)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
}
}
@@ -681,16 +729,45 @@ void Foam::isoSurface::generateTriPoints
}
-//template<class Type>
-//Foam::tmp<Foam::Field<Type>>
-//Foam::isoSurface::sample(const Field<Type>& vField) const
-//{
-// return tmp<Field<Type>>(new Field<Type>(vField, meshCells()));
-//}
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::isoSurface::interpolate
+(
+ const label nPoints,
+ const labelList& triPointMergeMap,
+ const DynamicList<Type>& unmergedValues
+)
+{
+ // One value per point
+ tmp<Field<Type> > tvalues(new Field<Type>(nPoints, pTraits<Type>::zero));
+ Field<Type>& values = tvalues.ref();
+ labelList nValues(values.size(), 0);
+
+ forAll(unmergedValues, i)
+ {
+ label mergedPointI = triPointMergeMap[i];
+
+ if (mergedPointI >= 0)
+ {
+ values[mergedPointI] += unmergedValues[i];
+ nValues[mergedPointI]++;
+ }
+ }
+
+ forAll(values, i)
+ {
+ if (nValues[i] > 0)
+ {
+ values[i] /= scalar(nValues[i]);
+ }
+ }
+
+ return tvalues;
+}
template<class Type>
-Foam::tmp<Foam::Field<Type>>
+Foam::tmp<Foam::Field<Type> >
Foam::isoSurface::interpolate
(
const GeometricField<Type, fvPatchField, volMesh>& cCoords,
@@ -707,7 +784,7 @@ Foam::isoSurface::interpolate
>> c2(adaptPatchFields(cCoords));
- DynamicList<Type> triPoints(nCutCells_);
+ DynamicList<Type> triPoints(3*nCutCells_);
DynamicList<label> triMeshCells(nCutCells_);
// Dummy snap data
@@ -731,52 +808,12 @@ Foam::isoSurface::interpolate
triMeshCells
);
-
- // One value per point
- tmp<Field<Type>> tvalues
+ return interpolate
(
- new Field<Type>(points().size(), pTraits<Type>::zero)
+ points().size(),
+ triPointMergeMap_,
+ triPoints
);
- Field<Type>& values = tvalues.ref();
- labelList nValues(values.size(), 0);
-
- forAll(triPoints, i)
- {
- label mergedPointI = triPointMergeMap_[i];
-
- if (mergedPointI >= 0)
- {
- values[mergedPointI] += triPoints[i];
- nValues[mergedPointI]++;
- }
- }
-
- if (debug)
- {
- Pout<< "nValues:" << values.size() << endl;
- label nMult = 0;
- forAll(nValues, i)
- {
- if (nValues[i] == 0)
- {
- FatalErrorInFunction
- << "point:" << i << " nValues:" << nValues[i]
- << abort(FatalError);
- }
- else if (nValues[i] > 1)
- {
- nMult++;
- }
- }
- Pout<< "Of which mult:" << nMult << endl;
- }
-
- forAll(values, i)
- {
- values[i] /= scalar(nValues[i]);
- }
-
- return tvalues;
}
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 3f8486efa4a..d789abed83f 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -124,18 +124,52 @@ void Foam::sampledIsoSurface::getIsoFields() const
// Get pointField
// ~~~~~~~~~~~~~~
+ // In case of multiple iso values we don't want to calculate multiple e.g.
+ // "volPointInterpolate(p)" so register it and re-use it. This is the
+ // same as the 'cache' functionality from volPointInterpolate but
+ // unfortunately that one does not guarantee that the field pointer
+ // remain: e.g. some other functionObject might delete the cached version.
+ // (volPointInterpolation::interpolate with cache=false deletes any
+ // registered one or if mesh.changing())
+
if (!subMeshPtr_.valid())
{
- word pointFldName = "volPointInterpolate(" + isoField_ + ')';
+ const word pointFldName =
+ "volPointInterpolate_"
+ + type()
+ + "("
+ + isoField_
+ + ')';
if (fvm.foundObject<pointScalarField>(pointFldName))
{
if (debug)
{
InfoInFunction
- << "Lookup pointField " << pointFldName << endl;
+ << "lookup pointField " << pointFldName << endl;
}
- pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
+ const pointScalarField& pfld = fvm.lookupObject<pointScalarField>
+ (
+ pointFldName
+ );
+
+ if (!pfld.upToDate(*volFieldPtr_))
+ {
+ if (debug)
+ {
+ InfoInFunction
+ << "updating pointField "
+ << pointFldName << endl;
+ }
+ // Update the interpolated value
+ volPointInterpolation::New(fvm).interpolate
+ (
+ *volFieldPtr_,
+ const_cast<pointScalarField&>(pfld)
+ );
+ }
+
+ pointFieldPtr_ = &pfld;
}
else
{
@@ -149,27 +183,20 @@ void Foam::sampledIsoSurface::getIsoFields() const
<< endl;
}
- if
+ // Interpolate without cache. Note that we're registering it
+ // below so next time round it goes into the condition
+ // above.
+ tmp<pointScalarField> tpfld
(
- storedPointFieldPtr_.empty()
- || (fvm.time().timeName() != storedPointFieldPtr_().instance())
- )
- {
- if (debug)
- {
- InfoInFunction
- << "Interpolating volField " << volFieldPtr_->name()
- << " to get pointField " << pointFldName << endl;
- }
-
- storedPointFieldPtr_.reset
+ volPointInterpolation::New(fvm).interpolate
(
- volPointInterpolation::New(fvm)
- .interpolate(*volFieldPtr_).ptr()
- );
- storedPointFieldPtr_->checkOut();
- pointFieldPtr_ = storedPointFieldPtr_.operator->();
- }
+ *volFieldPtr_,
+ pointFldName,
+ false
+ )
+ );
+ pointFieldPtr_ = tpfld.ptr();
+ const_cast<pointScalarField*>(pointFieldPtr_)->store();
}
@@ -427,7 +454,6 @@ Foam::sampledIsoSurface::sampledIsoSurface
prevTimeIndex_(-1),
storedVolFieldPtr_(NULL),
volFieldPtr_(NULL),
- storedPointFieldPtr_(NULL),
pointFieldPtr_(NULL)
{
if (!sampledSurface::interpolate())
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index d0a0b0ba986..e128560284d 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -94,7 +94,6 @@ class sampledIsoSurface
mutable const volScalarField* volFieldPtr_;
//- Cached pointfield
- mutable autoPtr<pointScalarField> storedPointFieldPtr_;
mutable const pointScalarField* pointFieldPtr_;
// And on subsetted mesh
--
GitLab