0001487: Slice Sampling (interpolation)on ANSAgenerated mesh fails (glibc free) for large parallel meshes (>100 mln cells)

ID	Project	Category	View Status	Date Submitted	Last Update

0001487	OpenFOAM	Bug	public	2015-01-13 15:41	2016-03-12 11:55

Reporter	matteoL	Assigned To	henry
Priority	normal	Severity	major	Reproducibility	sometimes
Status	resolved	Resolution	fixed
Platform	unix	OS	centos	OS Version	(please specify)

Summary	0001487: Slice Sampling (interpolation)on ANSAgenerated mesh fails (glibc free) for large parallel meshes (>100 mln cells)
Description	Hello, I am using OF2.2 and I am having trouble running a sampling slice when working with a volume mesh generated with ANSA. Here is my functionIncludeDict: slices { type surfaces; storeFilter false; functionObjectLibs ( "libsampling.so" ); outputControl timeStep; outputInterval 100; setFormat raw; surfaceFormat vtk; interpolationScheme cellPoint; fields U; surfaces ( X_Movie_0 { type cuttingPlane; planeType pointAndNormal; pointAndNormalDict { normalVector (1.000000 0.000000 0.000000); basePoint (-1.000000 -0.000000 -0.000000); } interpolate true; mergeTol 1e-8; } ); } When I run it (with execFlowFunctionObjects -latestTime ) I get the following error: Surfaces slices: Reading surface description: X_Movie_0 * glibc detected * execFlowFunctionObjects: double free or corruption (!prev): 0x00000000112848d0 * * glibc detected * execFlowFunctionObjects: free(): invalid next size (normal): 0x0000000010fe17e0 * * glibc detected * execFlowFunctionObjects: double free or corruption (!prev): 0x000000001087acd0 * * glibc detected * execFlowFunctionObjects: double free or corruption (!prev): 0x0000000013098d00 * * glibc detected * execFlowFunctionObjects: free(): invalid next size (normal): 0x00000000136d0880 * * glibc detected * execFlowFunctionObjects: double free or corruption (!prev): 0x0000000012675d60 * * glibc detected * execFlowFunctionObjects: double free or corruption (!prev): 0x0000000011f10520 * * glibc detected * execFlowFunctionObjects: malloc(): memory corruption: 0x0000000010bb2400 * * glibc detected * execFlowFunctionObjects: free(): invalid next size (normal): 0x00000000120b4550 * * glibc detected * execFlowFunctionObjects: free(): invalid next size (normal): 0x0000000013f0ea50 * ======= Backtrace: ========= /lib64/libc.so.6(+0x750c6)[0x2b08c8f590c6] Notice that: 1)The same dictionary works perfectly on a SnappyHexMesh generated mesh. 2)This dictionary works if I change the flag interpolate to false. 3)On a small test case, it works fine 4)On all my big mesh it fails I have tried the other interpolation schemes options but those fail as well..
Steps To Reproduce	It Always happens on my big meshes. Unfortunately I cannot share a mesh with the problem because it is too big to upload and for confidentiality reasons.
Additional Information	I would be glad to give you any further information required or perform any test you need. Thanks
Tags	No tags attached.

~~user4~~ 2015-01-14 08:55 ~0003537	- is your mesh ok? checkMesh -allGeometry - does it work if you slightly perturb the basePoint - maybe you're cutting along edges or faces. - set the mergeTol smaller - switch on the debug flag for isoSurface

matteoL 2015-01-14 10:33 reporter	StlFilesDebug.PNG (24,183 bytes) StlFilesDebug.PNG (24,183 bytes)

matteoL 2015-01-14 10:41 reporter ~0003543	Hello mattijs, Thanks for your help. here is the log of chekMesh: Checking geometry... Overall domain bounding box (-28.796305 -42.500068 -1.027792) (54.959113 3.7873575e-05 20.510296) Mesh (non-empty, non-wedge) directions (1 1 1) Mesh (non-empty) directions (1 1 1) Boundary openness (1.0127749e-15 -1.7754966e-16 -2.3397573e-15) OK. Max cell openness = 8.9573582e-16 OK. Max aspect ratio = 70.509793 OK. Minimum face area = 2.4086684e-10. Maximum face area = 1.7096958. Face area magnitudes OK. Min volume = 1.9372266e-14. Max volume = 2.057266. Total volume = 59463.08. Cell volumes OK. Mesh non-orthogonality Max: 85.118003 average: 16.411001 Number of severely non-orthogonal (> 70 degrees) faces: 278. Non-orthogonality check OK. <<Writing 278 non-orthogonal faces to set nonOrthoFaces Face pyramids OK. Max skewness = 5.105728, 24 highly skew faces detected which may impair the quality of the results <<Writing 24 skew faces to set skewFaces Coupled point location match (average 0) OK. Error in face tets: 1 faces with low quality or negative volume decomposition tets. <<Writing 1 faces with low quality or negative volume decomposition tets to set lowQualityTetFaces Edges too small, min/max edge length = 3.0707654e-06 2.3819928, number too small: 2011612 <<Writing 2011612 points on short edges to set shortEdges <<Writing 2600899 near (closer than 9.6359319e-05 apart) points to set nearPoints There are 51431 faces with concave angles between consecutive edges. Max concave angle = 56.166339 degrees. <<Writing 51431 faces with concave angles to set concaveFaces Face flatness (1 = flat, 0 = butterfly) : average = 0.99977716 min = 0.59449512 There are 2 faces with ratio between projected and actual area < 0.8 Minimum ratio (minimum flatness, maximum warpage) = 0.59449512 <<Writing 2 warped faces to set warpedFaces Cell determinant (wellposedness) : minimum: 0 average: 3.9231018 *Cells with small determinant (< 0.001) found, number of cells: 50913 <<Writing 50913 under-determined cells to set underdeterminedCells *Concave cells (using face planes) found, number of cells: 2231003 <<Writing 2231003 concave cells to set concaveCells Failed 4 mesh checks. I have reduced the mergetol to 1e-11, perturbed the basepoint and activated the debug flag. It is still crashing. I am uploading the log. With the debug it has created lots of stl files, most of them empty (I have uploaded a snapshot of my disk). Any idea? Thanks, again, matteo

matteoL 2015-01-14 10:45 reporter	WithDebug.txt.gz (37,034 bytes)

~~user4~~ 2015-01-19 20:50 ~0003567	From your traceback the problem is in the interpolation from volField to pointField (volPointInterpolation) and not related to the cutting plane. Do you see the same problem if e.g. you run foamToVTK? - do you see the same problem in 2.3.x? - what are your patch types? - try running a small problem under valgrind (using e.g. the mpirunDebug script)

matteoL 2015-01-20 10:36 reporter ~0003570	Hello, foamToVTK works fine. - yes. - my patch types are: wall, patch,cyclycAMI and processor - unfortunately I can't run it on valgrind because I have those issues only on big cases. I haven't been able to generate a small case with such problem. Indeed the problem is in the interpolation step, in fact if i run the sampling procedure without the interpolate option to TRUE everything works fine. Thanks.

wyldckat 2015-02-21 20:52 updater ~0003851	Greetings to all! @Matteo: Any chance you can test this case with the latest bug fix version of OpenFOAM 2.3.x? I ask this because with luck, this issue has been solved with the fix implemented for the bug report #1506: http://www.openfoam.org/mantisbt/view.php?id=1506 Best regards, Bruno

~~user4~~ 2015-03-11 12:21 ~0004082	We can replicate your problem. It is due to a dangling pointer in the iso surface algorithm which is used by the cuttingPlane sampling method. We expect to have a fix for this soon. As a workaround sometimes it helps to change the order of your functionObjects

matteoL 2015-08-20 15:09 reporter ~0005277	Hello Mattijs, Any news on this? Thanks, Matteo

wyldckat 2016-02-21 20:21 updater ~0005976	This seems to have been fixed in the following commit on the OpenFOAM-history repository by Mattijs: https://github.com/OpenCFD/OpenFOAM-history/commit/a11a6cd044f4d5ce4ac34926be6f777bbc2bf5ec @Matteo: Since you could not share the case with which this issue could be reproduced, can you please let us know if the commit I pointed out solves the problem? If not, any chance you can provide a dummy test case that reproduces the same error?

MattijsJ 2016-03-11 17:14 reporter	0001-BUG-isosurface-stores-reference-to-volPointInterpola.patch (67,596 bytes) From 13d3ace2604b97361771b32b8d3e82f99cf0957f Mon Sep 17 00:00:00 2001 From: mattijs <mattijs> Date: Fri, 11 Mar 2016 17:08:52 +0000 Subject: [PATCH] BUG: isosurface: stores reference to volPointInterpolation-cached field. If volPointInterpolation decides to remove field it invalidates the reference Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1487 --- .../sampledSurface/isoSurface/isoSurface.C \| 680 ++------------------- .../sampledSurface/isoSurface/isoSurface.H \| 80 +-- .../sampledSurface/isoSurface/isoSurfaceCell.C \| 218 ++----- .../sampledSurface/isoSurface/isoSurfaceCell.H \| 45 -- .../isoSurface/isoSurfaceCellTemplates.C \| 361 +++++------ .../isoSurface/isoSurfaceTemplates.C \| 289 +++++---- .../sampledSurface/isoSurface/sampledIsoSurface.C \| 72 ++- .../sampledSurface/isoSurface/sampledIsoSurface.H \| 1 - 8 files changed, 470 insertions(+), 1276 deletions(-) diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C index a9d57ad..1ca392f 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurface.C +++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C @@ -471,43 +471,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,73 +484,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 @@ -1202,435 +1098,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(3surf.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 +1420,58 @@ Foam::isoSurface::isoSurface } + { + DynamicList<point> triPoints(3nCutCells_); + 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 +1487,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 967536b..6937bfb 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurface.H +++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H @@ -193,20 +193,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 @@ -323,6 +311,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 +330,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 +340,10 @@ class isoSurface public: + //- Declare friendship with isoSurfaceCell to share some functionality + friend class isoSurfaceCell; + + //- Runtime type information TypeName("isoSurface"); @@ -413,18 +365,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 c2fec19..21359a8 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C +++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C @@ -1222,142 +1222,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 ( @@ -1605,57 +1469,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() diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H index 034fbd3..21b54f2 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H +++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H @@ -214,19 +214,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 +267,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 +292,6 @@ class isoSurfaceCell labelList& newToOldPoints ); - //- Combine all triangles inside a cell into a minimal triangulation - void combineCellTriangles(); public: @@ -348,24 +321,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 1eba206..ea26461 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(3nCutCells_); 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 ea92b65..295733c 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]; @@ -189,6 +189,8 @@ Type Foam::isoSurface::generatePoint } +// Note: cannot use simpler isoSurfaceCell::generateTriPoints since +// the need here to sometimes pass in remote 'snappedPoints' template<class Type> void Foam::isoSurface::generateTriPoints ( @@ -240,8 +242,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 +256,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 +278,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 +417,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; } } @@ -693,6 +743,43 @@ 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::isoSurface::interpolate +( const GeometricField<Type, fvPatchField, volMesh>& cCoords, const Field<Type>& pCoords ) const @@ -707,7 +794,7 @@ Foam::isoSurface::interpolate >> c2(adaptPatchFields(cCoords)); - DynamicList<Type> triPoints(nCutCells_); + DynamicList<Type> triPoints(3nCutCells_); DynamicList<label> triMeshCells(nCutCells_); // Dummy snap data @@ -731,52 +818,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 3f8486e..d789abe 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 d0a0b0b..e128560 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 -- 2.7.2 0001-BUG-isosurface-stores-reference-to-volPointInterpola.patch (67,596 bytes)

MattijsJ 2016-03-11 17:18 reporter ~0006031	Attached basically that OpenFOAM-history patch. Should also fix #905 (orientation of triangles) What was happening is that - isoSurface stores a reference to the pointField - this pointField was cached by volPointInterpolation - so if for any reason volPointInterpolation decided to create a new pointField the reference would become invalid. I tried to replicate the problem but failed.

henry 2016-03-12 11:54 manager ~0006032	Resolved in OpenFOAM-dev by commit 36f19c88c4362a28cfc7b13c70ab4f0431416997

Date Modified	Username	Field	Change
2015-01-13 15:41	matteoL	New Issue
2015-01-14 08:55	~~user4~~	Note Added: 0003537
2015-01-14 10:33	matteoL	File Added: StlFilesDebug.PNG
2015-01-14 10:41	matteoL	Note Added: 0003543
2015-01-14 10:45	matteoL	File Added: WithDebug.txt.gz
2015-01-19 20:50	~~user4~~	Note Added: 0003567
2015-01-20 10:36	matteoL	Note Added: 0003570
2015-02-21 20:52	wyldckat	Note Added: 0003851
2015-03-11 12:21	~~user4~~	Note Added: 0004082
2015-03-24 00:17	liuhuafei	Issue cloned: 0001621
2015-08-20 15:09	matteoL	Note Added: 0005277
2016-02-21 20:21	wyldckat	Note Added: 0005976
2016-03-11 17:14	MattijsJ	File Added: 0001-BUG-isosurface-stores-reference-to-volPointInterpola.patch
2016-03-11 17:18	MattijsJ	Note Added: 0006031
2016-03-12 11:54	henry	Note Added: 0006032
2016-03-12 11:54	henry	Status	new => resolved
2016-03-12 11:54	henry	Resolution	open => fixed
2016-03-12 11:54	henry	Assigned To	=> henry

View Issue Details

Activities

Issue History