diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C index f3821e1d91..888b51b9c6 100644 --- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C +++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C @@ -45,6 +45,11 @@ namespace Foam void Foam::distanceSurface::createGeometry() { + if (debug) + { + Pout<< "distanceSurface::createGeometry :updating geometry." << endl; + } + // Clear any stored topologies facesPtr_.clear(); @@ -67,7 +72,7 @@ void Foam::distanceSurface::createGeometry() false ), fvm, - dimensionedScalar("zero", dimless/dimTime, 0) + dimensionedScalar("zero", dimLength, 0) ) ); volScalarField& cellDistance = cellDistancePtr_(); @@ -157,6 +162,7 @@ void Foam::distanceSurface::createGeometry() } } + // On processor patches the mesh.C() will already be the cell centre // on the opposite side so no need to swap cellDistance. @@ -164,23 +170,70 @@ void Foam::distanceSurface::createGeometry() // Distance to points pointDistance_.setSize(fvm.nPoints()); { + const pointField& pts = fvm.points(); + List nearest; surfPtr_().findNearest ( - fvm.points(), - scalarField(fvm.nPoints(), GREAT), + pts, + scalarField(pts.size(), GREAT), nearest ); - forAll(pointDistance_, pointI) + + if (signed_) { - pointDistance_[pointI] = Foam::mag - ( - nearest[pointI].hitPoint() - - fvm.points()[pointI] - ); + vectorField normal; + surfPtr_().getNormal(nearest, normal); + + forAll(nearest, i) + { + vector d(pts[i]-nearest[i].hitPoint()); + + if ((d&normal[i]) > 0) + { + pointDistance_[i] = Foam::mag(d); + } + else + { + pointDistance_[i] = -Foam::mag(d); + } + } + } + else + { + forAll(nearest, i) + { + pointDistance_[i] = Foam::mag(pts[i]-nearest[i].hitPoint()); + } } } + + if (debug) + { + Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl; + cellDistance.write(); + pointScalarField pDist + ( + IOobject + ( + "pointDistance", + fvm.time().timeName(), + fvm.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + pointMesh::New(fvm), + dimensionedScalar("zero", dimLength, 0) + ); + pDist.internalField() = pointDistance_; + + Pout<< "Writing point distance:" << pDist.objectPath() << endl; + pDist.write(); + } + + //- Direct from cell field and point field. isoSurfPtr_.reset ( @@ -196,6 +249,7 @@ void Foam::distanceSurface::createGeometry() if (debug) { print(Pout); + Pout<< endl; } } @@ -264,6 +318,13 @@ bool Foam::distanceSurface::needsUpdate() const bool Foam::distanceSurface::expire() { + if (debug) + { + Pout<< "distanceSurface::expire :" + << " have-facesPtr_:" << facesPtr_.valid() + << " needsUpdate_:" << needsUpdate_ << endl; + } + // Clear any stored topologies facesPtr_.clear(); @@ -280,6 +341,13 @@ bool Foam::distanceSurface::expire() bool Foam::distanceSurface::update() { + if (debug) + { + Pout<< "distanceSurface::update :" + << " have-facesPtr_:" << facesPtr_.valid() + << " needsUpdate_:" << needsUpdate_ << endl; + } + if (!needsUpdate_) { return false; diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C index 8522dd8442..3571df4a49 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurface.C +++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C @@ -58,6 +58,31 @@ Foam::scalar Foam::isoSurface::isoFraction } +bool Foam::isoSurface::isEdgeOfFaceCut +( + const scalarField& pVals, + const face& f, + const bool ownLower, + const bool neiLower +) const +{ + forAll(f, fp) + { + bool fpLower = (pVals[f[fp]] < iso_); + if + ( + (fpLower != ownLower) + || (fpLower != neiLower) + || (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)) + ) + { + return true; + } + } + return false; +} + + // Determine for every face/cell whether it (possibly) generates triangles. void Foam::isoSurface::calcCutTypes ( @@ -84,22 +109,13 @@ void Foam::isoSurface::calcCutTypes } else { - // Mesh edge. + // See if any mesh edge is cut by looping over all the edges of the + // face. const face f = mesh_.faces()[faceI]; - forAll(f, fp) + if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower)) { - bool fpLower = (pVals[f[fp]] < iso_); - if - ( - (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)) - || (fpLower != ownLower) - || (fpLower != neiLower) - ) - { - faceCutType_[faceI] = CUT; - break; - } + faceCutType_[faceI] = CUT; } } } @@ -117,22 +133,13 @@ void Foam::isoSurface::calcCutTypes { bool ownLower = (cVals[own[faceI]] < iso_); - // Mesh edge. const face f = mesh_.faces()[faceI]; - forAll(f, fp) + if (isEdgeOfFaceCut(pVals, f, ownLower, ownLower)) { - bool fpLower = (pVals[f[fp]] < iso_); - if - ( - (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)) - || (fpLower != ownLower) - ) - { - faceCutType_[faceI] = CUT; - break; - } + faceCutType_[faceI] = CUT; } + faceI++; } } @@ -152,19 +159,9 @@ void Foam::isoSurface::calcCutTypes // Mesh edge. const face f = mesh_.faces()[faceI]; - forAll(f, fp) + if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower)) { - bool fpLower = (pVals[f[fp]] < iso_); - if - ( - (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)) - || (fpLower != ownLower) - || (fpLower != neiLower) - ) - { - faceCutType_[faceI] = CUT; - break; - } + faceCutType_[faceI] = CUT; } } faceI++; @@ -1355,6 +1352,30 @@ Foam::isoSurface::isoSurface iso_(iso), mergeDistance_(mergeTol*mesh_.bounds().mag()) { + const polyBoundaryMesh& patches = mesh_.boundaryMesh(); + const labelList& own = mesh_.faceOwner(); + + // Check + forAll(patches, patchI) + { + if (isA(patches[patchI])) + { + FatalErrorIn + ( + "isoSurface::isoSurface\n" + "(\n" + " const volScalarField& cVals,\n" + " const scalarField& pVals,\n" + " const scalar iso,\n" + " const bool regularise,\n" + " const scalar mergeTol\n" + ")\n" + ) << "Iso surfaces not supported on case with empty patches." + << exit(FatalError); + } + } + + // Determine if any cut through face/cell calcCutTypes(cVals, pVals); @@ -1364,8 +1385,6 @@ Foam::isoSurface::isoSurface PackedList<1> isBoundaryPoint(mesh_.nPoints()); labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces()); - const polyBoundaryMesh& patches = mesh_.boundaryMesh(); - const labelList& own = mesh_.faceOwner(); forAll(patches, patchI) { diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H index 832773018e..abe3d51b56 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurface.H +++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H @@ -31,6 +31,9 @@ Description G.M. Treece, R.W. Prager and A.H. Gee. Note: + - not possible on patches of type 'empty'. There are no values on + 'empty' patch fields so even the api would have to change + (no volScalarField as argument). Too messy. - in parallel the regularisation (coarsening) always takes place and slightly different surfaces will be created compared to non-parallel. The surface will still be continuous though! @@ -123,6 +126,15 @@ class isoSurface const scalar s1 ) const; + //- Check if any edge of a face is cut + bool isEdgeOfFaceCut + ( + const scalarField& pVals, + const face& f, + const bool ownLower, + const bool neiLower + ) const; + //- Set faceCutType,cellCutType. void calcCutTypes ( @@ -217,7 +229,7 @@ class isoSurface ) const; template - label generateTriPoints + label generateFaceTriPoints ( const volScalarField& cVals, const scalarField& pVals, diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C index ffff99f104..67d05a5a1c 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C +++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C @@ -200,7 +200,7 @@ void Foam::isoSurface::generateTriPoints template -Foam::label Foam::isoSurface::generateTriPoints +Foam::label Foam::isoSurface::generateFaceTriPoints ( const volScalarField& cVals, const scalarField& pVals, @@ -288,6 +288,19 @@ void Foam::isoSurface::generateTriPoints const labelList& own = mesh_.faceOwner(); const labelList& nei = mesh_.faceNeighbour(); + if + ( + (cVals.size() != mesh_.nCells()) + || (pVals.size() != mesh_.nPoints()) + || (cCoords.size() != mesh_.nCells()) + || (pCoords.size() != mesh_.nPoints()) + || (snappedCc.size() != mesh_.nCells()) + || (snappedPoint.size() != mesh_.nPoints()) + ) + { + FatalErrorIn("isoSurface::generateTriPoints(..)") + << "Incorrect size." << abort(FatalError); + } // Determine neighbouring snap status labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1); @@ -319,7 +332,7 @@ void Foam::isoSurface::generateTriPoints { if (faceCutType_[faceI] != NOTCUT) { - generateTriPoints + generateFaceTriPoints ( cVals, pVals, @@ -357,7 +370,7 @@ void Foam::isoSurface::generateTriPoints { if (faceCutType_[faceI] != NOTCUT) { - generateTriPoints + generateFaceTriPoints ( cVals, pVals, @@ -384,14 +397,14 @@ void Foam::isoSurface::generateTriPoints } else if (isA(pp)) { - // Assume zero-gradient. + // Assume zero-gradient. But what about coordinates? label faceI = pp.start(); forAll(pp, i) { if (faceCutType_[faceI] != NOTCUT) { - generateTriPoints + generateFaceTriPoints ( cVals, pVals, @@ -423,7 +436,7 @@ void Foam::isoSurface::generateTriPoints { if (faceCutType_[faceI] != NOTCUT) { - generateTriPoints + generateFaceTriPoints ( cVals, pVals,