Skip to content
GitLab
Commit 36f19c88 authored by Henry Weller's avatar Henry Weller
Browse files

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
parent 3ba1d176
Expand all Hide whitespace changes
Inline Side-by-side
src/sampling/sampledSurface/isoSurface/isoSurface.H
View file @ 36f19c88
......@@ -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>
......
src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
View file @ 36f19c88
......@@ -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);
}
}
......
src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
View file @ 36f19c88
......@@ -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
......
src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
View file @ 36f19c88
......@@ -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));
}