diff --git a/applications/utilities/mesh/generation/blockMesh/blockPoints.C b/applications/utilities/mesh/generation/blockMesh/blockPoints.C index 0cb024215dbcf57f76a754fb2f7c9900fe802bb8..2fecb21bf54b2340e3e56afd9b5acc408d1d56a2 100644 --- a/applications/utilities/mesh/generation/blockMesh/blockPoints.C +++ b/applications/utilities/mesh/generation/blockMesh/blockPoints.C @@ -41,23 +41,23 @@ namespace Foam void block::blockPoints() { // set local variables for mesh specification - label ni = blockDef_.n().x(); - label nj = blockDef_.n().y(); - label nk = blockDef_.n().z(); + const label ni = blockDef_.n().x(); + const label nj = blockDef_.n().y(); + const label nk = blockDef_.n().z(); - point start = blockDef_.points()[blockDef_.blockShape()[0]]; - point xEnd = blockDef_.points()[blockDef_.blockShape()[1]]; - point xyEnd = blockDef_.points()[blockDef_.blockShape()[2]]; - point yEnd = blockDef_.points()[blockDef_.blockShape()[3]]; + const point p000 = blockDef_.points()[blockDef_.blockShape()[0]]; + const point p100 = blockDef_.points()[blockDef_.blockShape()[1]]; + const point p110 = blockDef_.points()[blockDef_.blockShape()[2]]; + const point p010 = blockDef_.points()[blockDef_.blockShape()[3]]; - point zEnd = blockDef_.points()[blockDef_.blockShape()[4]]; - point xzEnd = blockDef_.points()[blockDef_.blockShape()[5]]; - point xyzEnd = blockDef_.points()[blockDef_.blockShape()[6]]; - point yzEnd = blockDef_.points()[blockDef_.blockShape()[7]]; + const point p001 = blockDef_.points()[blockDef_.blockShape()[4]]; + const point p101 = blockDef_.points()[blockDef_.blockShape()[5]]; + const point p111 = blockDef_.points()[blockDef_.blockShape()[6]]; + const point p011 = blockDef_.points()[blockDef_.blockShape()[7]]; - // set reference to the list of edge point and weighting factors - const List<List<point> >& edgePoints = blockDef_.blockEdgePoints(); - const scalarListList& edgeWeights = blockDef_.blockEdgeWeights(); + // list of edge point and weighting factors + const List<List<point> >& p = blockDef_.blockEdgePoints(); + const scalarListList& w = blockDef_.blockEdgeWeights(); // generate vertices @@ -69,193 +69,138 @@ void block::blockPoints() { label vertexNo = vtxLabel(i, j, k); - vector edgex1 = start*(1.0 - edgeWeights[0][i]) - + xEnd*edgeWeights[0][i]; + // points on edges + vector edgex1 = p000 + (p100 - p000)*w[0][i]; + vector edgex2 = p010 + (p110 - p010)*w[1][i]; + vector edgex3 = p011 + (p111 - p011)*w[2][i]; + vector edgex4 = p001 + (p101 - p001)*w[3][i]; - vector edgex2 = yEnd*(1.0 - edgeWeights[1][i]) - + xyEnd*edgeWeights[1][i]; + vector edgey1 = p000 + (p010 - p000)*w[4][j]; + vector edgey2 = p100 + (p110 - p100)*w[5][j]; + vector edgey3 = p101 + (p111 - p101)*w[6][j]; + vector edgey4 = p001 + (p011 - p001)*w[7][j]; - vector edgex3 = yzEnd*(1.0 - edgeWeights[2][i]) - + xyzEnd*edgeWeights[2][i]; - - vector edgex4 = zEnd*(1.0 - edgeWeights[3][i]) - + xzEnd*edgeWeights[3][i]; - - - - vector edgey1 = start*(1.0 - edgeWeights[4][j]) - + yEnd*edgeWeights[4][j]; - - vector edgey2 = xEnd*(1.0 - edgeWeights[5][j]) - + xyEnd*edgeWeights[5][j]; - - vector edgey3 = xzEnd*(1.0 - edgeWeights[6][j]) - + xyzEnd*edgeWeights[6][j]; - - vector edgey4 = zEnd*(1.0 - edgeWeights[7][j]) - + yzEnd*edgeWeights[7][j]; - - - - vector edgez1 = start*(1.0 - edgeWeights[8][k]) - + zEnd*edgeWeights[8][k]; - - vector edgez2 = xEnd*(1.0 - edgeWeights[9][k]) - + xzEnd*edgeWeights[9][k]; - - vector edgez3 = xyEnd*(1.0 - edgeWeights[10][k]) - + xyzEnd*edgeWeights[10][k]; - - vector edgez4 = yEnd*(1.0 - edgeWeights[11][k]) - + yzEnd*edgeWeights[11][k]; + vector edgez1 = p000 + (p001 - p000)*w[8][k]; + vector edgez2 = p100 + (p101 - p100)*w[9][k]; + vector edgez3 = p110 + (p111 - p110)*w[10][k]; + vector edgez4 = p010 + (p011 - p010)*w[11][k]; // calculate the importance factors for all edges // x - direction scalar impx1 = ( - (1.0 - edgeWeights[0][i]) - *(1.0 - edgeWeights[4][j]) - *(1.0 - edgeWeights[8][k]) - + edgeWeights[0][i] - *(1.0 - edgeWeights[5][j]) - *(1.0 - edgeWeights[9][k]) + (1.0 - w[0][i])*(1.0 - w[4][j])*(1.0 - w[8][k]) + + w[0][i]*(1.0 - w[5][j])*(1.0 - w[9][k]) ); scalar impx2 = ( - (1.0 - edgeWeights[1][i]) - *edgeWeights[4][j] - *(1.0 - edgeWeights[11][k]) - + edgeWeights[1][i] - *edgeWeights[5][j] - *(1.0 - edgeWeights[10][k]) + (1.0 - w[1][i])*w[4][j]*(1.0 - w[11][k]) + + w[1][i]*w[5][j]*(1.0 - w[10][k]) ); - scalar impx3 = - ( - (1.0 - edgeWeights[2][i]) - *edgeWeights[7][j] - *edgeWeights[11][k] - + edgeWeights[2][i] - *edgeWeights[6][j] - *edgeWeights[10][k] - ); + scalar impx3 = + ( + (1.0 - w[2][i])*w[7][j]*w[11][k] + + w[2][i]*w[6][j]*w[10][k] + ); - scalar impx4 = - ( - (1.0 - edgeWeights[3][i]) - *(1.0 - edgeWeights[7][j]) - *edgeWeights[8][k] - + edgeWeights[3][i] - *(1.0 - edgeWeights[6][j]) - *edgeWeights[9][k] - ); + scalar impx4 = + ( + (1.0 - w[3][i])*(1.0 - w[7][j])*w[8][k] + + w[3][i]*(1.0 - w[6][j])*w[9][k] + ); + scalar magImpx = impx1 + impx2 + impx3 + impx4; + impx1 /= magImpx; + impx2 /= magImpx; + impx3 /= magImpx; + impx4 /= magImpx; // y - direction scalar impy1 = ( - (1.0 - edgeWeights[4][j]) - *(1.0 - edgeWeights[0][i]) - *(1.0 - edgeWeights[8][k]) - + edgeWeights[4][j] - *(1.0 - edgeWeights[1][i]) - *(1.0 - edgeWeights[11][k]) + (1.0 - w[4][j])*(1.0 - w[0][i])*(1.0 - w[8][k]) + + w[4][j]*(1.0 - w[1][i])*(1.0 - w[11][k]) ); scalar impy2 = ( - (1.0 - edgeWeights[5][j]) - *edgeWeights[0][i] - *(1.0 - edgeWeights[9][k]) - + edgeWeights[5][j] - *edgeWeights[1][i] - *(1.0 - edgeWeights[10][k]) + (1.0 - w[5][j])*w[0][i]*(1.0 - w[9][k]) + + w[5][j]*w[1][i]*(1.0 - w[10][k]) ); scalar impy3 = ( - (1.0 - edgeWeights[6][j]) - *edgeWeights[3][i] - *edgeWeights[9][k] - + edgeWeights[6][j] - *edgeWeights[2][i] - *edgeWeights[10][k] + (1.0 - w[6][j])*w[3][i]*w[9][k] + + w[6][j]*w[2][i]*w[10][k] ); scalar impy4 = ( - (1.0 - edgeWeights[7][j]) - *(1.0 - edgeWeights[3][i]) - *edgeWeights[8][k] - + edgeWeights[7][j] - *(1.0 - edgeWeights[2][i]) - *edgeWeights[11][k] + (1.0 - w[7][j])*(1.0 - w[3][i])*w[8][k] + + w[7][j]*(1.0 - w[2][i])*w[11][k] ); + scalar magImpy = impy1 + impy2 + impy3 + impy4; + impy1 /= magImpy; + impy2 /= magImpy; + impy3 /= magImpy; + impy4 /= magImpy; // z - direction scalar impz1 = ( - (1.0 - edgeWeights[8][k]) - *(1.0 - edgeWeights[0][i]) - *(1.0 - edgeWeights[4][j]) - + edgeWeights[8][k] - *(1.0 - edgeWeights[3][i]) - *(1.0 - edgeWeights[7][j]) + (1.0 - w[8][k])*(1.0 - w[0][i])*(1.0 - w[4][j]) + + w[8][k]*(1.0 - w[3][i])*(1.0 - w[7][j]) ); scalar impz2 = ( - (1.0 - edgeWeights[9][k]) - *edgeWeights[0][i] - *(1.0 - edgeWeights[5][j]) - + edgeWeights[9][k] - *edgeWeights[3][i] - *(1.0 - edgeWeights[6][j]) + (1.0 - w[9][k])*w[0][i]*(1.0 - w[5][j]) + + w[9][k]*w[3][i]*(1.0 - w[6][j]) ); scalar impz3 = ( - (1.0 - edgeWeights[10][k]) - *edgeWeights[1][i] - *edgeWeights[5][j] - + edgeWeights[10][k] - *edgeWeights[2][i] - *edgeWeights[6][j] + (1.0 - w[10][k])*w[1][i]*w[5][j] + + w[10][k]*w[2][i]*w[6][j] ); scalar impz4 = ( - (1.0 - edgeWeights[11][k]) - *(1.0 - edgeWeights[1][i]) - *edgeWeights[4][j] - + edgeWeights[11][k] - *(1.0 - edgeWeights[2][i]) - *edgeWeights[7][j] + (1.0 - w[11][k])*(1.0 - w[1][i])*w[4][j] + + w[11][k]*(1.0 - w[2][i])*w[7][j] ); + scalar magImpz = impz1 + impz2 + impz3 + impz4; + impz1 /= magImpz; + impz2 /= magImpz; + impz3 /= magImpz; + impz4 /= magImpz; + // calculate the correction vectors - vector corx1 = impx1*(edgePoints[0][i] - edgex1); - vector corx2 = impx2*(edgePoints[1][i] - edgex2); - vector corx3 = impx3*(edgePoints[2][i] - edgex3); - vector corx4 = impx4*(edgePoints[3][i] - edgex4); + vector corx1 = impx1*(p[0][i] - edgex1); + vector corx2 = impx2*(p[1][i] - edgex2); + vector corx3 = impx3*(p[2][i] - edgex3); + vector corx4 = impx4*(p[3][i] - edgex4); - vector cory1 = impy1*(edgePoints[4][j] - edgey1); - vector cory2 = impy2*(edgePoints[5][j] - edgey2); - vector cory3 = impy3*(edgePoints[6][j] - edgey3); - vector cory4 = impy4*(edgePoints[7][j] - edgey4); + vector cory1 = impy1*(p[4][j] - edgey1); + vector cory2 = impy2*(p[5][j] - edgey2); + vector cory3 = impy3*(p[6][j] - edgey3); + vector cory4 = impy4*(p[7][j] - edgey4); - vector corz1 = impz1*(edgePoints[8][k] - edgez1); - vector corz2 = impz2*(edgePoints[9][k] - edgez2); - vector corz3 = impz3*(edgePoints[10][k] - edgez3); - vector corz4 = impz4*(edgePoints[11][k] - edgez4); + vector corz1 = impz1*(p[8][k] - edgez1); + vector corz2 = impz2*(p[9][k] - edgez2); + vector corz3 = impz3*(p[10][k] - edgez3); + vector corz4 = impz4*(p[11][k] - edgez4); // multiply by the importance factor + // x - direction edgex1 *= impx1; edgex2 *= impx2; @@ -285,7 +230,6 @@ void block::blockPoints() vertices_[vertexNo] += corx1 + corx2 + corx3 + corx4; vertices_[vertexNo] += cory1 + cory2 + cory3 + cory4; vertices_[vertexNo] += corz1 + corz2 + corz3 + corz4; - } } } diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict index b3e489a1763decc5f17f62569d386ccf38dafd52..db42299d1eba33304fecd603edfaf95f03880477 100644 --- a/applications/utilities/postProcessing/sampling/sample/sampleDict +++ b/applications/utilities/postProcessing/sampling/sample/sampleDict @@ -40,7 +40,7 @@ surfaceFormat vtk; // 1] vertex values determined from neighbouring cell-centre values // 2] face values determined using the current face interpolation scheme // for the field (linear, gamma, etc.) -interpolationScheme cellPointFace; +interpolationScheme cellPoint; // Fields to sample. fields @@ -155,21 +155,25 @@ surfaces // Optional: whether to leave as faces (=default) or triangulate } -/* not yet (re)implemented -- - constantIso + interpolatedIso { - name iso; - field rho; - value 0.5; + // Iso surface for interpolated values only + type isoSurface; + isoField rho; + isoValue 0.5; + interpolate true; + //regularise false; //optional: do not simplify } - someIso + constantIso { - type iso; - field rho; - value 0.5; - interpolate true; + // Iso surface for constant values. Guarantees triangles to not + // cross cells. + type isoSurfaceCell; + isoField rho; + isoValue 0.5; + interpolate false; + //regularise false; //optional: do not simplify } - */ ); diff --git a/applications/utilities/surface/surfaceCheck/surfaceCheck.C b/applications/utilities/surface/surfaceCheck/surfaceCheck.C index 3506b05de7942769a2a425868ba0ddb2f7a51f59..2cd26f757e6ed75f3e9a00516f15b354eacb3a03 100644 --- a/applications/utilities/surface/surfaceCheck/surfaceCheck.C +++ b/applications/utilities/surface/surfaceCheck/surfaceCheck.C @@ -36,7 +36,7 @@ License using namespace Foam; // Does face use valid vertices? -bool validTri(const triSurface& surf, const label faceI) +bool validTri(const bool verbose, const triSurface& surf, const label faceI) { // Simple check on indices ok. @@ -49,20 +49,21 @@ bool validTri(const triSurface& surf, const label faceI) || (f[2] < 0) || (f[2] >= surf.points().size()) ) { - //WarningIn("validTri(const triSurface&, const label)") - // << "triangle " << faceI << " vertices " << f - // << " uses point indices outside point range 0.." - // << surf.points().size()-1 << endl; + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI << " vertices " << f + << " uses point indices outside point range 0.." + << surf.points().size()-1 << endl; return false; } if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2])) { - //WarningIn("validTri(const triSurface&, const label)") - // << "triangle " << faceI - // << " uses non-unique vertices " << f - // << endl; + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI + << " uses non-unique vertices " << f + << " coords:" << f.points(surf.points()) + << endl; return false; } @@ -91,11 +92,12 @@ bool validTri(const triSurface& surf, const label faceI) && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2])) ) { - //WarningIn("validTri(const triSurface&, const label)") - // << "triangle " << faceI << " vertices " << f - // << " has the same vertices as triangle " << nbrFaceI - // << " vertices " << nbrF - // << endl; + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI << " vertices " << f + << " has the same vertices as triangle " << nbrFaceI + << " vertices " << nbrF + << " coords:" << f.points(surf.points()) + << endl; return false; } @@ -170,9 +172,11 @@ int main(int argc, char *argv[]) argList::validArgs.clear(); argList::validArgs.append("surface file"); argList::validOptions.insert("noSelfIntersection", ""); + argList::validOptions.insert("verbose", ""); argList args(argc, argv); bool checkSelfIntersection = !args.options().found("noSelfIntersection"); + bool verbose = args.options().found("verbose"); fileName surfFileName(args.additionalArgs()[0]); Pout<< "Reading surface from " << surfFileName << " ..." << nl << endl; @@ -232,7 +236,7 @@ int main(int argc, char *argv[]) forAll(surf, faceI) { - if (!validTri(surf, faceI)) + if (!validTri(verbose, surf, faceI)) { illegalFaces.append(faceI); } diff --git a/applications/utilities/surface/surfaceSubset/surfaceSubsetDict b/applications/utilities/surface/surfaceSubset/surfaceSubsetDict index 70879af6f428700fe0b8d8f9dcd7cbd325732327..83f5399b4431f23bb58557ed8eab60fba36f3888 100644 --- a/applications/utilities/surface/surfaceSubset/surfaceSubsetDict +++ b/applications/utilities/surface/surfaceSubset/surfaceSubsetDict @@ -42,4 +42,7 @@ surface // Extend selection with edge neighbours addFaceNeighbours no; +// Invert selection +invertSelection false; + // ************************************************************************* // diff --git a/src/OpenFOAM/containers/Lists/SortableList/SortableList.C b/src/OpenFOAM/containers/Lists/SortableList/SortableList.C index 7ebfa9de9ef185d82ee77f4e1e0c7c07d0d3921d..9eaf99c801b994d8f2a93f31c3ca2b8b8a83ca03 100644 --- a/src/OpenFOAM/containers/Lists/SortableList/SortableList.C +++ b/src/OpenFOAM/containers/Lists/SortableList/SortableList.C @@ -32,7 +32,7 @@ License // Construct from List template <class Type> -Foam::SortableList<Type>::SortableList(const List<Type>& values) +Foam::SortableList<Type>::SortableList(const UList<Type>& values) : List<Type>(values), indices_(values.size()) diff --git a/src/OpenFOAM/containers/Lists/SortableList/SortableList.H b/src/OpenFOAM/containers/Lists/SortableList/SortableList.H index 1a5272e5b7823098dec984927c1ede3d23066286..305da6e32eed03ef9194d42a1247925e992b203f 100644 --- a/src/OpenFOAM/containers/Lists/SortableList/SortableList.H +++ b/src/OpenFOAM/containers/Lists/SortableList/SortableList.H @@ -88,7 +88,7 @@ public: //- Construct from List, sorting the elements. // Starts with indices set to index in argument - explicit SortableList(const List<Type>&); + explicit SortableList(const UList<Type>&); //- Construct from tranferred List, sorting the elements. // Starts with indices set to index in argument diff --git a/src/OpenFOAM/fields/Fields/Field/Field.C b/src/OpenFOAM/fields/Fields/Field/Field.C index 994f8ee6301b63f66c26e5532b5c755c2f86dece..c906bcbf5559b9e601109e01f4f10281af24b96c 100644 --- a/src/OpenFOAM/fields/Fields/Field/Field.C +++ b/src/OpenFOAM/fields/Fields/Field/Field.C @@ -575,20 +575,6 @@ void Field<Type>::replace } -template<class Type> -void Field<Type>::transfer(Field<Type>& f) -{ - List<Type>::transfer(f); -} - - -template<class Type> -void Field<Type>::transfer(List<Type>& lst) -{ - List<Type>::transfer(lst); -} - - template<class Type> tmp<Field<Type> > Field<Type>::T() const { diff --git a/src/OpenFOAM/fields/Fields/Field/Field.H b/src/OpenFOAM/fields/Fields/Field/Field.H index 8e168cb0ac9b5911a8b52da259615bf751365cca..1b372925e340ea410302b1162ad6598c7b65ccee 100644 --- a/src/OpenFOAM/fields/Fields/Field/Field.H +++ b/src/OpenFOAM/fields/Fields/Field/Field.H @@ -297,12 +297,6 @@ public: //- Replace a component field of the field void replace(const direction, const cmptType&); - //- Transfer the contents of the argument Field into this Field - void transfer(Field<Type>&); - - //- Transfer the contents of the argument List into this Field - void transfer(List<Type>&); - //- Return the field transpose (only defined for second rank tensors) tmp<Field<Type> > T() const; diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointFaces.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointFaces.C index dc8d865f59cf8a48175efd0e7e4b3f9946a19302..b7a90375b7ce1cfba075922e47529dce67fd1ea7 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointFaces.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointFaces.C @@ -38,6 +38,11 @@ const labelListList& primitiveMesh::pointFaces() const { if (!pfPtr_) { + if (debug) + { + Pout<< "primitiveMesh::pointFaces() : " + << "calculating pointFaces" << endl; + } // Invert faces() pfPtr_ = new labelListList(nPoints()); invertManyToMany(nPoints(), faces(), *pfPtr_); diff --git a/src/Pstream/mpi/IPread.C b/src/Pstream/mpi/IPread.C index 487749b24b23bdf94787e72c8bfda3caf871f89a..ac7c2e787a958138a0ffb2f772e3f68ae4a111a0 100644 --- a/src/Pstream/mpi/IPread.C +++ b/src/Pstream/mpi/IPread.C @@ -193,15 +193,13 @@ void Foam::IPstream::waitRequests() { if (IPstream_outstandingRequests_.size() > 0) { - List<MPI_Status> status(IPstream_outstandingRequests_.size()); - if ( MPI_Waitall ( IPstream_outstandingRequests_.size(), IPstream_outstandingRequests_.begin(), - status.begin() + MPI_STATUSES_IGNORE ) ) { @@ -231,9 +229,7 @@ bool Foam::IPstream::finishedRequest(const label i) } int flag; - MPI_Status status; - - MPI_Test(&IPstream_outstandingRequests_[i], &flag, &status); + MPI_Test(&IPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE); return flag != 0; } diff --git a/src/Pstream/mpi/OPwrite.C b/src/Pstream/mpi/OPwrite.C index 86d9c998ab6360bb37f5ec08ea2ecdc7ced1ecdf..d0cdcbd0bcc28041fc258c9884bd9e815d61b549 100644 --- a/src/Pstream/mpi/OPwrite.C +++ b/src/Pstream/mpi/OPwrite.C @@ -131,15 +131,13 @@ void Foam::OPstream::waitRequests() { if (OPstream_outstandingRequests_.size() > 0) { - List<MPI_Status> status(OPstream_outstandingRequests_.size()); - if ( MPI_Waitall ( OPstream_outstandingRequests_.size(), OPstream_outstandingRequests_.begin(), - status.begin() + MPI_STATUSES_IGNORE ) ) { @@ -169,9 +167,7 @@ bool Foam::OPstream::finishedRequest(const label i) } int flag; - MPI_Status status; - - MPI_Test(&OPstream_outstandingRequests_[i], &flag, &status); + MPI_Test(&OPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE); return flag != 0; } diff --git a/src/Pstream/mpi/Pstream.C b/src/Pstream/mpi/Pstream.C index a6361ab102077eb21922a224d99dbc4ed86b7ba2..caa21b6dcc000bc7cc6b026afd346ee35b00e674 100644 --- a/src/Pstream/mpi/Pstream.C +++ b/src/Pstream/mpi/Pstream.C @@ -157,8 +157,6 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop) if (Pstream::nProcs() <= Pstream::nProcsSimpleSum) { - MPI_Status status; - if (Pstream::master()) { for @@ -180,7 +178,7 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop) Pstream::procID(slave), Pstream::msgType(), MPI_COMM_WORLD, - &status + MPI_STATUS_IGNORE ) ) { @@ -260,7 +258,7 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop) Pstream::procID(Pstream::masterNo()), Pstream::msgType(), MPI_COMM_WORLD, - &status + MPI_STATUS_IGNORE ) ) { @@ -279,8 +277,6 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop) Value = sum; /* - MPI_Status status; - int myProcNo = Pstream::myProcNo(); int nProcs = Pstream::nProcs(); @@ -314,7 +310,7 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop) Pstream::procID(childProcId), Pstream::msgType(), MPI_COMM_WORLD, - &status + MPI_STATUS_IGNORE ) ) { @@ -370,7 +366,7 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop) Pstream::procID(parentId), Pstream::msgType(), MPI_COMM_WORLD, - &status + MPI_STATUS_IGNORE ) ) { diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C index f117968fe11e11823e73274b8da827f8cea2db25..4cfba30501e21a3676290aa06deba01e5f505cc6 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C @@ -2875,6 +2875,8 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2 } // 2. Extend to 2:1. I don't understand yet why this is not done + // 2. Extend to 2:1. For non-cube cells the scalar distance does not work + // so make sure it at least provides 2:1. PackedList<1> refineCell(mesh_.nCells(), 0); forAll(allCellInfo, cellI) { @@ -2887,6 +2889,25 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2 } faceConsistentRefinement(true, refineCell); + while (true) + { + label nChanged = faceConsistentRefinement(true, refineCell); + + reduce(nChanged, sumOp<label>()); + + if (debug) + { + Pout<< "hexRef8::consistentSlowRefinement2 : Changed " << nChanged + << " refinement levels due to 2:1 conflicts." + << endl; + } + + if (nChanged == 0) + { + break; + } + } + // 3. Convert back to labelList. label nRefined = 0; diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 2f301c040687b8da341e9e392dcbaa202a1c2c5f..199daec5544e3c8a1a0d8f66ac34be4fb8144f38 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -297,15 +297,4 @@ cfdTools/general/SRF/SRFModel/SRFModel/newSRFModel.C cfdTools/general/SRF/SRFModel/rpm/rpm.C cfdTools/general/SRF/derivedFvPatchFields/SRFVelocityFvPatchVectorField/SRFVelocityFvPatchVectorField.C -fvMeshCutSurface = fvMesh/fvMeshCutSurface - -meshCut = $(fvMeshCutSurface)/meshCut -$(meshCut)/meshCutSurface.C -$(meshCut)/cellAddressing.C - -edgeCuts = $(fvMeshCutSurface)/edgeCuts -$(edgeCuts)/meshEdgeCuts.C -$(edgeCuts)/faceDecompIsoSurfaceCuts.C -$(edgeCuts)/cellDecompIsoSurfaceCuts.C - LIB = $(FOAM_LIBBIN)/libfiniteVolume diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C index 00ce11355a3ceebabe00dffd59f1605453f76ee2..321eeba6a87be806089b4cafd53ea2195286ca56 100644 --- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C +++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C @@ -33,15 +33,68 @@ void Foam::setRefCell const volScalarField& field, const dictionary& dict, label& refCelli, - scalar& refValue + scalar& refValue, + bool forceReference ) { - if (field.needReference()) + if (field.needReference() || forceReference) { word refCellName = field.name() + "RefCell"; + word refPointName = field.name() + "RefPoint"; + word refValueName = field.name() + "RefValue"; - refCelli = readLabel(dict.lookup(refCellName)); + if (dict.found(refCellName)) + { + if (Pstream::master()) + { + refCelli = readLabel(dict.lookup(refCellName)); + } + else + { + refCelli = -1; + } + } + else if (dict.found(refPointName)) + { + point refPointi(dict.lookup(refPointName)); + refCelli = field.mesh().findCell(refPointi); + label hasRef = (refCelli >= 0 ? 1 : 0); + label sumHasRef = returnReduce<label>(hasRef, sumOp<label>()); + if (sumHasRef != 1) + { + FatalErrorIn + ( + "void Foam::setRefCell" + "(" + " const volScalarField&," + " const dictionary&," + " label& scalar&," + " bool" + ")" + ) + << "Unable to set reference cell for field " << field.name() + << nl << " Reference point " << refPointName + << " found on multiple domains" << nl << abort(FatalError); + } + } + else + { + FatalErrorIn + ( + "void Foam::setRefCell" + "(" + " const volScalarField&," + " const dictionary&," + " label& scalar&," + " bool" + ")" + ) + << "Unable to set reference cell for field" << field.name() << nl + << " Please supply either " << refCellName + << " or " << refPointName << nl << abort(FatalError); + } + refValue = readScalar(dict.lookup(refValueName)); } } diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H index 724f44ed175186e6c765dc4b3fd7472cace9b41c..db9858ce27aed5b35cd2518281721d71c6c524ad 100644 --- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H +++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.H @@ -52,7 +52,8 @@ void setRefCell const volScalarField& field, const dictionary& dict, label& refCelli, - scalar& refValue + scalar& refValue, + bool forceReference = false ); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C index 279e72a0c55620bc270bd21ba3d65c022529434c..ab3af1717e006ce47754d5affaa9b7d93de5bc3a 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/directMappedFixedValue/directMappedFixedValueFvPatchField.C @@ -234,6 +234,16 @@ void directMappedFixedValueFvPatchField<Type>::updateCoeffs() ( mpp.samplePatch() ); + if (patchID < 0) + { + FatalErrorIn + ( + "void directMappedFixedValueFvPatchField<Type>::" + "updateCoeffs()" + )<< "Unable to find sample patch " << mpp.samplePatch() + << " for patch " << this->patch().name() << nl + << abort(FatalError); + } typedef GeometricField<Type, fvPatchField, volMesh> fieldType; const word& fieldName = this->dimensionedInternalField().name(); const fieldType& sendField = diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C index 9dab83d168a4e57c7041331ca349b95dd1fdcf44..b8698eb17738d2e5e39adfd6ac0eaa3568356a15 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C @@ -483,7 +483,7 @@ void Foam::fvMatrix<Type>::setReference { if (psi_.needReference() || forceReference) { - if (Pstream::master()) + if (cell >= 0) { source()[cell] += diag()[cell]*value; diag()[cell] += diag()[cell]; diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompCuts.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompCuts.H deleted file mode 100644 index a1cc320e8fe6dfd8aaf8052428564f7546ec6e2d..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompCuts.H +++ /dev/null @@ -1,161 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::cellDecompCuts - -Description - Container for cuts of edges of (implicit) tet decomposition. - Used to collect data for meshCut. - - As much as possible, cuts are defined using mesh information: - - cut (exactly) through mesh vertex - - cut (exactly) through cell centre - - - cut through mesh edge. Both edge label and position on edge given. - - - cut through tet pyramidEdge (edge between vertex and cell centre). - Edge and position on edge given. - - - cut through diagonalEdge (edge between vertices of a face) - Edge and position on edge given. - -SourceFiles - -\*---------------------------------------------------------------------------*/ - -#ifndef cellDecompCuts_H -#define cellDecompCuts_H - -#include "meshEdgeCuts.H" -#include "pyramidEdge.H" -#include "diagonalEdge.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes - -/*---------------------------------------------------------------------------*\ - Class cellDecompCuts Declaration -\*---------------------------------------------------------------------------*/ - -class cellDecompCuts -: - public meshEdgeCuts -{ - -protected: - - labelList meshCellCentres_; - - List<pyramidEdge> pyrEdges_; - scalarField pyrEdgeWeights_; - - List<diagonalEdge> diagEdges_; - scalarField diagEdgeWeights_; - - // Private Member Functions - - -public: - - // Constructors - - //- Construct from components - cellDecompCuts - ( - const primitiveMesh& mesh, - const labelList& cells, - - const labelList& meshVerts, - const labelList& meshCellCentres, - - const labelList& meshEdges, - const scalarField& meshEdgeWeights, - - const List<pyramidEdge>& pyrEdges, - const scalarField& pyrEdgeWeights, - - const List<diagonalEdge>& diagEdges, - const scalarField& diagEdgeWeights - ) - : - meshEdgeCuts(mesh, cells, meshVerts, meshEdges, meshEdgeWeights), - meshCellCentres_(meshCellCentres), - pyrEdges_(pyrEdges), - pyrEdgeWeights_(pyrEdgeWeights), - diagEdges_(diagEdges), - diagEdgeWeights_(diagEdgeWeights) - {} - - - // Member Functions - - const labelList& meshCellCentres() const - { - return meshCellCentres_; - } - - const List<pyramidEdge>& pyrEdges() const - { - return pyrEdges_; - } - - const scalarField& pyrEdgeWeights() const - { - return pyrEdgeWeights_; - } - - const List<diagonalEdge>& diagEdges() const - { - return diagEdges_; - } - - const scalarField& diagEdgeWeights() const - { - return diagEdgeWeights_; - } - - label size() const - { - return - meshEdgeCuts::size() + meshCellCentres_.size() - + pyrEdges_.size() + diagEdges_.size(); - } - -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompIsoSurfaceCuts.C b/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompIsoSurfaceCuts.C deleted file mode 100644 index 1f8f7d6d1fc516fca8c4bb30dfc81d6159e322a1..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompIsoSurfaceCuts.C +++ /dev/null @@ -1,267 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - -\*---------------------------------------------------------------------------*/ - -#include "cellDecompIsoSurfaceCuts.H" -#include "volPointInterpolation.H" - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::cellDecompIsoSurfaceCuts::constructEdgeCuts -( - const volScalarField& volField, - const pointScalarField& ptField, - const scalar isoVal, - const scalar tol -) -{ - const primitiveMesh& mesh = volField.mesh(); - - // Intermediate storage for labels of cut edges and cut points - label nBFaces = mesh.nFaces() - mesh.nInternalFaces(); - - labelHashSet cutCells(nBFaces); - labelHashSet cutVerts(nBFaces); - labelHashSet cutCellCentres(nBFaces); - - DynamicList<label> meshCutEdges(nBFaces); - DynamicList<scalar> meshEdgeWeights(nBFaces); - - DynamicList<pyramidEdge> cutPyrEdges(nBFaces); - DynamicList<scalar> pyrEdgeWeights(nBFaces); - - DynamicList<diagonalEdge> cutDiagEdges(nBFaces); - DynamicList<scalar> diagEdgeWeights(nBFaces); - - - // Walk over faces - forAll(mesh.faces(), faceI) - { - const face& f = mesh.faces()[faceI]; - - label ownerI = mesh.faceOwner()[faceI]; - label neighbourI = -1; - if (mesh.isInternalFace(faceI)) - { - neighbourI = mesh.faceNeighbour()[faceI]; - } - - scalar weight; - - // - // Check face diagonal. Simple triangulation from f[0] - // - - for(label fp = 2; fp < f.size() - 1; fp++) - { - if (crosses(isoVal, ptField[f[0]], ptField[f[fp]], weight)) - { - mark(ownerI, cutCells); - if (neighbourI != -1) - { - mark(neighbourI, cutCells); - } - - if (weight < tol) - { - mark(f[0], cutVerts); - } - else if (weight > 1-tol) - { - mark(f[fp], cutVerts); - } - else - { - cutDiagEdges.append(diagonalEdge(faceI, 0, fp)); - diagEdgeWeights.append(weight); - } - } - } - } - - - // Walk over all mesh points - forAll(mesh.points(), pointI) - { - const labelList& myCells = mesh.pointCells()[pointI]; - - forAll(myCells, myCellI) - { - label cellI = myCells[myCellI]; - - // - // Check pyramid edge crossing - // - - scalar weight; - - if (crosses(isoVal, ptField[pointI], volField[cellI], weight)) - { - mark(cellI, cutCells); - - if (weight < tol) - { - mark(pointI, cutVerts); - } - else if (weight > 1-tol) - { - mark(cellI, cutCellCentres); - } - else - { - cutPyrEdges.append(pyramidEdge(pointI, cellI)); - pyrEdgeWeights.append(weight); - } - } - } - } - - - // Walk over all mesh edges - forAll(mesh.edges(), edgeI) - { - const edge& e = mesh.edges()[edgeI]; - - scalar weight; - - if (crosses(isoVal, ptField[e.start()], ptField[e.end()], weight)) - { - mark(mesh.edgeCells()[edgeI], cutCells); - - if (weight < tol) - { - mark(e.start(), cutVerts); - } - else if (weight > 1-tol) - { - mark(e.end(), cutVerts); - } - else - { - meshCutEdges.append(edgeI); - meshEdgeWeights.append(weight); - } - } - } - - meshCutEdges.shrink(); - meshEdgeWeights.shrink(); - - cutPyrEdges.shrink(); - pyrEdgeWeights.shrink(); - - cutDiagEdges.shrink(); - diagEdgeWeights.shrink(); - - - // Tranfer lists to cellDecompCuts - - cells_ = cutCells.toc(); - - meshVerts_ = cutVerts.toc(); - meshCellCentres_ = cutCellCentres.toc(); - - meshEdges_.transfer(meshCutEdges); - meshEdgeWeights_.transfer(meshEdgeWeights); - - pyrEdges_.transfer(cutPyrEdges); - pyrEdgeWeights_.transfer(pyrEdgeWeights); - - diagEdges_.transfer(cutDiagEdges); - diagEdgeWeights_.transfer(diagEdgeWeights); -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::cellDecompIsoSurfaceCuts::cellDecompIsoSurfaceCuts -( - const volScalarField& volField, - const pointScalarField& ptField, - const scalar isoVal, - const scalar tol -) -: - cellDecompCuts - ( - volField.mesh(), - labelList(), - - labelList(), // mesh vertices - labelList(), // mesh cell centres - - labelList(), // mesh edges - scalarField(), - - List<pyramidEdge>(), // pyramid edges - scalarField(), - - List<diagonalEdge>(), // face diagonal edges - scalarField() - ) -{ - constructEdgeCuts(volField, ptField, isoVal, tol); -} - - -// Construct from interpolator and field (does interpolation) -Foam::cellDecompIsoSurfaceCuts::cellDecompIsoSurfaceCuts -( - const volScalarField& volField, - const volPointInterpolation& pInterp, - const scalar isoVal, - const scalar tol -) -: - cellDecompCuts - ( - volField.mesh(), - labelList(), - - labelList(), // mesh vertices - labelList(), // mesh cell centres - - labelList(), // mesh edges - scalarField(), - - List<pyramidEdge>(), // pyramid edges - scalarField(), - - List<diagonalEdge>(), // face diagonal edges - scalarField() - ) -{ - // Get field on vertices - tmp<pointScalarField> tptField = pInterp.interpolate(volField); - const pointScalarField& ptField = tptField(); - - - constructEdgeCuts(volField, ptField, isoVal, tol); -} - - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompIsoSurfaceCuts.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompIsoSurfaceCuts.H deleted file mode 100644 index f7ee4d2b1159e69150f4b50067d0ffc3c20224dd..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/cellDecompIsoSurfaceCuts.H +++ /dev/null @@ -1,101 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::cellDecompIsoSurfaceCuts - -Description - cellCuts resulting from isoSurface reconstruction on cellDecomp tet mesh - -SourceFiles - cellDecompIsoSurfaceCuts.C - -\*---------------------------------------------------------------------------*/ - -#ifndef cellDecompIsoSurfaceCuts_H -#define cellDecompIsoSurfaceCuts_H - -#include "cellDecompCuts.H" -#include "volFields.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes -class volPointInterpolation; - -/*---------------------------------------------------------------------------*\ - Class cellDecompIsoSurfaceCuts Declaration -\*---------------------------------------------------------------------------*/ - -class cellDecompIsoSurfaceCuts -: - public cellDecompCuts -{ - // Private Member Functions - - //- Do all the hard work: determine edge crossings - void constructEdgeCuts - ( - const volScalarField&, - const pointScalarField&, - const scalar, - const scalar - ); - -public: - - // Constructors - - //- Construct from field on cells and field on vertices - cellDecompIsoSurfaceCuts - ( - const volScalarField&, - const pointScalarField&, - const scalar isoValue, - const scalar tol - ); - - //- Construct from interpolator and field (does interpolation) - cellDecompIsoSurfaceCuts - ( - const volScalarField&, - const volPointInterpolation&, - const scalar isoValue, - const scalar tol - ); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompCuts.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompCuts.H deleted file mode 100644 index c5f374e051266d817aa21ba62024e25c1878a363..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompCuts.H +++ /dev/null @@ -1,200 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::faceDecompCuts - -Description - Container for cuts of edges of (implicit) tet decomposition. Used to - collect data for meshCut. As much as possible cuts are defined using - mesh information: - - - cut (exactly) through mesh vertex - - ,, ,, cell centre - - ,, ,, face centre - - - cut through mesh edge. Both edge label and position on edge given. - - - cut through tet pyramidEdge (edge between vertex and cell centre). - Edge and position on edge given. - - cut through tet faceEdge (edge between vertex and face centre). - Edge and position on edge given. - - cut through tet centreEdge (edge between face centre and cell centre). - Edge and position on edge given. - -SourceFiles - -\*---------------------------------------------------------------------------*/ - -#ifndef faceDecompCuts_H -#define faceDecompCuts_H - -#include "meshEdgeCuts.H" -#include "pyramidEdge.H" -#include "faceEdge.H" -#include "centreEdge.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes - -/*---------------------------------------------------------------------------*\ - Class faceDecompCuts Declaration -\*---------------------------------------------------------------------------*/ - -class faceDecompCuts -: - public meshEdgeCuts -{ - -protected: - - //- List of faces whose centre is exactly cut - labelList meshFaceCentres_; - - //- List of cells whose centre is exactly cut - labelList meshCellCentres_; - - //- List of pyramid edge descriptions - List<pyramidEdge> pyrEdges_; - - //- For each pyramid edge description the position of the cut - scalarField pyrEdgeWeights_; - - //- List of face edge descriptions - List<faceEdge> faceEdges_; - - //- For each face edge description the position of the cut - scalarField faceEdgeWeights_; - - //- List of centre edge descriptions - List<centreEdge> centreEdges_; - - //- For each centre edge description the position of the cut - scalarField centreEdgeWeights_; - -public: - - // Constructors - - //- Construct from components - faceDecompCuts - ( - const primitiveMesh& mesh, - const labelList& cells, - - const labelList& meshVerts, - const labelList& meshFaceCentres, - const labelList& meshCellCentres, - - const labelList& meshEdges, - const scalarField& meshEdgeWeights, - - const List<pyramidEdge>& pyrEdges, - const scalarField& pyrEdgeWeights, - - const List<faceEdge>& faceEdges, - const scalarField& faceEdgeWeights, - - const List<centreEdge>& centreEdges, - const scalarField& centreEdgeWeights - ) - : - meshEdgeCuts(mesh, cells, meshVerts, meshEdges, meshEdgeWeights), - meshFaceCentres_(meshFaceCentres), - meshCellCentres_(meshCellCentres), - pyrEdges_(pyrEdges), - pyrEdgeWeights_(pyrEdgeWeights), - faceEdges_(faceEdges), - faceEdgeWeights_(faceEdgeWeights), - centreEdges_(centreEdges), - centreEdgeWeights_(centreEdgeWeights) - {} - - - // Member Functions - - const labelList& meshFaceCentres() const - { - return meshFaceCentres_; - } - - const labelList& meshCellCentres() const - { - return meshCellCentres_; - } - - const List<pyramidEdge>& pyrEdges() const - { - return pyrEdges_; - } - - const scalarField& pyrEdgeWeights() const - { - return pyrEdgeWeights_; - } - - const List<centreEdge>& centreEdges() const - { - return centreEdges_; - } - - const scalarField& centreEdgeWeights() const - { - return centreEdgeWeights_; - } - - const List<faceEdge>& faceEdges() const - { - return faceEdges_; - } - - const scalarField& faceEdgeWeights() const - { - return faceEdgeWeights_; - } - - label size() const - { - return - meshEdgeCuts::size() - + meshFaceCentres_.size() + meshCellCentres_.size() - + pyrEdges_.size() + centreEdges_.size() + faceEdges_.size(); - } - -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.C b/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.C deleted file mode 100644 index 6f2ad0a8f14a99148039821f16f07681d9a9f731..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.C +++ /dev/null @@ -1,342 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - -\*---------------------------------------------------------------------------*/ - -#include "faceDecompIsoSurfaceCuts.H" -#include "volPointInterpolation.H" - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::faceDecompIsoSurfaceCuts::constructEdgeCuts -( - const volScalarField& volField, - const pointScalarField& ptField, - const scalar isoVal, - const scalar tol -) -{ - const primitiveMesh& mesh = volField.mesh(); - - // Intermediate storage for labels of cut edges and cut points - label nBFaces = mesh.nFaces() - mesh.nInternalFaces(); - - labelHashSet cutCells(nBFaces); - labelHashSet cutVerts(nBFaces); - labelHashSet cutFaceCentres(nBFaces); - labelHashSet cutCellCentres(nBFaces); - - DynamicList<label> meshCutEdges(nBFaces); - DynamicList<scalar> meshEdgeWeights(nBFaces); - - DynamicList<pyramidEdge> cutPyrEdges(nBFaces); - DynamicList<scalar> pyrEdgeWeights(nBFaces); - - DynamicList<faceEdge> cutFaceEdges(nBFaces); - DynamicList<scalar> faceEdgeWeights(nBFaces); - - DynamicList<centreEdge> cutCentreEdges(nBFaces); - DynamicList<scalar> centreEdgeWeights(nBFaces); - - - // Walk over faces - forAll(mesh.faces(), faceI) - { - const face& f = mesh.faces()[faceI]; - - // Face centre value - scalar fcVal = 0; - - forAll(f, fp) - { - fcVal += ptField[f[fp]]; - } - fcVal /= f.size(); - - - label ownerI = mesh.faceOwner()[faceI]; - label neighbourI = -1; - if (mesh.isInternalFace(faceI)) - { - neighbourI = mesh.faceNeighbour()[faceI]; - } - - scalar weight; - - // - // Check faceCentre-cellCentre crossing for owner and neigbour - // - - if (crosses(isoVal, fcVal, volField[ownerI], weight)) - { - mark(ownerI, cutCells); - - if (weight < tol) - { - mark(faceI, cutFaceCentres); - } - else if (weight > 1-tol) - { - mark(ownerI, cutCellCentres); - } - else - { - cutCentreEdges.append(centreEdge(faceI, true)); - centreEdgeWeights.append(weight); - } - } - - if (neighbourI != -1) - { - if (crosses(isoVal, fcVal, volField[neighbourI], weight)) - { - mark(neighbourI, cutCells); - - if (weight < tol) - { - mark(faceI, cutFaceCentres); - } - else if (weight > 1-tol) - { - mark(neighbourI, cutCellCentres); - } - else - { - cutCentreEdges.append(centreEdge(faceI, false)); - centreEdgeWeights.append(weight); - } - } - } - - - forAll(f, fp) - { - // - // Check face internal (faceEdge) edge crossing - // - - if (crosses(isoVal, ptField[f[fp]], fcVal, weight)) - { - mark(ownerI, cutCells); - if (neighbourI != -1) - { - mark(neighbourI, cutCells); - } - - if (weight < tol) - { - mark(f[fp], cutVerts); - } - else if (weight > 1-tol) - { - mark(faceI, cutFaceCentres); - } - else - { - cutFaceEdges.append(faceEdge(faceI, fp)); - faceEdgeWeights.append(weight); - } - } - } - } - - - // Walk over all mesh points - forAll(mesh.points(), pointI) - { - const labelList& myCells = mesh.pointCells()[pointI]; - - forAll(myCells, myCellI) - { - label cellI = myCells[myCellI]; - - // - // Check pyramid edge crossing - // - - scalar weight; - - if (crosses(isoVal, ptField[pointI], volField[cellI], weight)) - { - mark(cellI, cutCells); - - if (weight < tol) - { - mark(pointI, cutVerts); - } - else if (weight > 1-tol) - { - mark(cellI, cutCellCentres); - } - else - { - cutPyrEdges.append(pyramidEdge(pointI, cellI)); - pyrEdgeWeights.append(weight); - } - } - } - } - - - // Walk over all mesh edges - forAll(mesh.edges(), edgeI) - { - const edge& e = mesh.edges()[edgeI]; - - scalar weight; - - if (crosses(isoVal, ptField[e.start()], ptField[e.end()], weight)) - { - mark(mesh.edgeCells()[edgeI], cutCells); - - if (weight < tol) - { - mark(e.start(), cutVerts); - } - else if (weight > 1-tol) - { - mark(e.end(), cutVerts); - } - else - { - meshCutEdges.append(edgeI); - meshEdgeWeights.append(weight); - } - } - } - - meshCutEdges.shrink(); - meshEdgeWeights.shrink(); - - cutPyrEdges.shrink(); - pyrEdgeWeights.shrink(); - - cutFaceEdges.shrink(); - faceEdgeWeights.shrink(); - - cutCentreEdges.shrink(); - centreEdgeWeights.shrink(); - - - - // Tranfer lists to faceDecompCuts - - cells_ = cutCells.toc(); - meshVerts_ = cutVerts.toc(); - meshFaceCentres_ = cutFaceCentres.toc(); - meshCellCentres_ = cutCellCentres.toc(); - - meshEdges_.transfer(meshCutEdges); - meshEdgeWeights_.transfer(meshEdgeWeights); - - pyrEdges_.transfer(cutPyrEdges); - pyrEdgeWeights_.transfer(pyrEdgeWeights); - - faceEdges_.transfer(cutFaceEdges); - faceEdgeWeights_.transfer(faceEdgeWeights); - - centreEdges_.transfer(cutCentreEdges); - centreEdgeWeights_.transfer(centreEdgeWeights); -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::faceDecompIsoSurfaceCuts::faceDecompIsoSurfaceCuts -( - const volScalarField& volField, - const pointScalarField& ptField, - const scalar isoVal, - const scalar tol -) -: - faceDecompCuts - ( - volField.mesh(), - labelList(), - - labelList(), - labelList(), - labelList(), - - labelList(), // mesh edges - scalarField(), - - List<pyramidEdge>(), // pyramid edges - scalarField(), - - List<faceEdge>(), // face edges - scalarField(), - - List<centreEdge>(), // cell centre edges - scalarField() - ) -{ - constructEdgeCuts(volField, ptField, isoVal, tol); -} - - -// Construct from interpolator and field (does interpolation) -Foam::faceDecompIsoSurfaceCuts::faceDecompIsoSurfaceCuts -( - const volScalarField& volField, - const volPointInterpolation& pInterp, - const scalar isoVal, - const scalar tol -) -: - faceDecompCuts - ( - volField.mesh(), - labelList(), - - labelList(), - labelList(), - labelList(), - - labelList(), // mesh edges - scalarField(), - - List<pyramidEdge>(), // pyramid edges - scalarField(), - - List<faceEdge>(), // face edges - scalarField(), - - List<centreEdge>(), // cell centre edges - scalarField() - ) -{ - // Get field on vertices - tmp<pointScalarField> tptField = pInterp.interpolate(volField); - const pointScalarField& ptField = tptField(); - - - constructEdgeCuts(volField, ptField, isoVal, tol); -} - - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.H deleted file mode 100644 index d845de2b4faa5896b6b7f2616100e2c426cc7c03..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/faceDecompIsoSurfaceCuts.H +++ /dev/null @@ -1,102 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::faceDecompIsoSurfaceCuts - -Description - List of triangles based on isoSurface reconstruction on faceDecomp mesh. - -SourceFiles - faceDecompIsoSurfaceCuts.C - -\*---------------------------------------------------------------------------*/ - -#ifndef faceDecompIsoSurfaceCuts_H -#define faceDecompIsoSurfaceCuts_H - -#include "faceDecompCuts.H" -#include "labelHashSet.H" -#include "volFields.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes -class volPointInterpolation; - -/*---------------------------------------------------------------------------*\ - Class faceDecompIsoSurfaceCuts Declaration -\*---------------------------------------------------------------------------*/ - -class faceDecompIsoSurfaceCuts -: - public faceDecompCuts -{ - // Private Member Functions - - //- Do all the hard work: determine edge crossings - void constructEdgeCuts - ( - const volScalarField&, - const pointScalarField&, - const scalar isoVal, - const scalar tol - ); - -public: - - // Constructors - - //- Construct from field on cells and field on vertices - faceDecompIsoSurfaceCuts - ( - const volScalarField&, - const pointScalarField&, - const scalar isoVal, - const scalar tol - ); - - //- Construct from interpolator and field (does interpolation) - faceDecompIsoSurfaceCuts - ( - const volScalarField&, - const volPointInterpolation&, - const scalar isoVal, - const scalar tol - ); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/meshEdgeCuts.C b/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/meshEdgeCuts.C deleted file mode 100644 index 752aee07ff810df68a2fc8e34d7e4707f0b9d671..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/meshEdgeCuts.C +++ /dev/null @@ -1,109 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - -\*---------------------------------------------------------------------------*/ - -#include "meshEdgeCuts.H" - -// * * * * * * * * * * * * * Private Static Functions * * * * * * * * * * * // - -void Foam::meshEdgeCuts::mark(const label elem, labelHashSet& markedElems) -{ - labelHashSet::const_iterator iter = markedElems.find(elem); - - if (iter == markedElems.end()) - { - markedElems.insert(elem); - } -} - - -void Foam::meshEdgeCuts::mark -( - const labelList& elems, - labelHashSet& markedElems -) -{ - forAll(elems, elemI) - { - mark(elems[elemI], markedElems); - } -} - - -bool Foam::meshEdgeCuts::crosses -( - const scalar isoVal, - const scalar val0, - const scalar val1, - scalar& weight -) -{ - if - ( - ((val0 <= isoVal) && (val1 >= isoVal)) - || ((val1 <= isoVal) && (val0 >= isoVal)) - ) - { - if (val0 == val1) - { - return false; - } - else - { - weight = (isoVal - val0)/((val1 - val0) + VSMALL); - - return true; - } - } - else - { - return false; - } -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -// Construct from components -Foam::meshEdgeCuts::meshEdgeCuts -( - const primitiveMesh& mesh, - const labelList& cells, - const labelList& meshVerts, - const labelList& meshEdges, - const scalarField& meshEdgeWeights -) -: - mesh_(mesh), - cells_(cells), - meshVerts_(meshVerts), - meshEdges_(meshEdges), - meshEdgeWeights_(meshEdgeWeights) -{} - - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/meshEdgeCuts.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/meshEdgeCuts.H deleted file mode 100644 index dab50b77e7d60385863fe6651418378a32c242bc..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/edgeCuts/meshEdgeCuts.H +++ /dev/null @@ -1,161 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::meshEdgeCuts - -Description - Container for cuts of edges of mesh. Cuts either exactly through existing - vertices or through edges. - -SourceFiles - -\*---------------------------------------------------------------------------*/ - -#ifndef meshEdgeCuts_H -#define meshEdgeCuts_H - -#include "labelList.H" -#include "labelHashSet.H" -#include "scalarField.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes -class primitiveMesh; - -/*---------------------------------------------------------------------------*\ - Class meshEdgeCuts Declaration -\*---------------------------------------------------------------------------*/ - -class meshEdgeCuts -{ - // Private data - - const primitiveMesh& mesh_; - - // Private Member Functions - - //- Disallow default bitwise copy construct - meshEdgeCuts(const meshEdgeCuts&); - - //- Disallow default bitwise assignment - void operator=(const meshEdgeCuts&); - - -protected: - - //- List of cells containing the cuts - labelList cells_; - - //- Points exactly cut by cuts - labelList meshVerts_; - - //- List of edge labels cut - labelList meshEdges_; - - //- Positions on edges - scalarField meshEdgeWeights_; - - - // Protected Static Functions - - //- Mark element in a hashSet - static void mark(const label elem, labelHashSet& markedElems); - - //- Mark list of elements in a hashSet - static void mark(const labelList& elems, labelHashSet& markedElems); - - //- Return true and set weight if linear interpolation between - // val0 and val1 crosses isoVal. weight=1 if isoVal==val1 - static bool crosses - ( - const scalar isoVal, - const scalar val0, - const scalar val1, - scalar& weight - ); - -public: - - // Constructors - - //- Construct from components - meshEdgeCuts - ( - const primitiveMesh& mesh, - const labelList& cells, - const labelList& meshVerts, - const labelList& meshEdges, - const scalarField& meshEdgeWeights - ); - - - // Member Functions - - const primitiveMesh& mesh() const - { - return mesh_; - } - - const labelList& cells() const - { - return cells_; - } - - const labelList& meshVerts() const - { - return meshVerts_; - } - - const labelList& meshEdges() const - { - return meshEdges_; - } - - const scalarField& meshEdgeWeights() const - { - return meshEdgeWeights_; - } - - label size() const - { - return meshVerts_.size() + meshEdges_.size(); - } - -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/cellAddressing.C b/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/cellAddressing.C deleted file mode 100644 index be1fbb99df3bf14272858280f3b50f44a033ad12..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/cellAddressing.C +++ /dev/null @@ -1,132 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - -\*---------------------------------------------------------------------------*/ - -#include "cellAddressing.H" -#include "labelHashSet.H" -#include "cellModel.H" - - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -Foam::labelList Foam::cellAddressing::makeIdentity(const label size) -{ - labelList result(size); - forAll(result, resultI) - { - result[resultI] = resultI; - } - return result; -} - - -Foam::label Foam::cellAddressing::findEdge(const edge& e) -{ - forAll(edges_, edgeI) - { - if (edges_[edgeI] == e) - { - return edgeI; - } - } - FatalErrorIn - ( - "cellAddressing::findEdge(const edge&)" - ) << "Problem: cannot find edge " << e << " in edges " << edges_ - << exit(FatalError); - - return -1; -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -// Construct from components -Foam::cellAddressing::cellAddressing(const Foam::cellModel& model) -: - modelFaces_(model.modelFaces()), - edges_(model.edges(makeIdentity(model.nPoints()))), - edgeFaces_(model.nEdges()), - faceEdges_(model.nFaces()), - pointEdges_(model.nPoints()) -{ - List<labelHashSet > dynEdgeFaces(edgeFaces_.size()); - - forAll(faceEdges_, faceI) - { - DynamicList<label> fEdges; - - const face& f = modelFaces_[faceI]; - - forAll(f, fp) - { - label edgeI = findEdge(edge(f[fp], f[(fp+1)%f.size()])); - - // Append to faceEdges - fEdges.append(edgeI); - - // Append to edgefaces - labelHashSet& eFaces = dynEdgeFaces[edgeI]; - - if (!eFaces.found(faceI)) - { - eFaces.insert(faceI); - } - } - fEdges.shrink(); - - // Copy into faceEdges_ - faceEdges_[faceI].transfer(fEdges.shrink()); - } - - // Convert dynEdgeFaces into edgeFaces_ - - forAll(dynEdgeFaces, edgeI) - { - edgeFaces_[edgeI] = dynEdgeFaces[edgeI].toc(); - } - - - // Convert edges_ into pointEdges - List<DynamicList<label> > dynPointEdges(model.nPoints()); - - forAll(edges_, edgeI) - { - const edge& e = edges_[edgeI]; - - dynPointEdges[e.start()].append(edgeI); - dynPointEdges[e.end()].append(edgeI); - } - - forAll(dynPointEdges, pointI) - { - pointEdges_[pointI].transfer(dynPointEdges[pointI].shrink()); - } -} - - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/cellAddressing.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/cellAddressing.H deleted file mode 100644 index d51937b4b41ce4a9479b175426dfbd1d891b8eab..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/cellAddressing.H +++ /dev/null @@ -1,145 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::cellAddressing - -Description - Additional addressing for a cellModel. (cellModel themselves only - contain cell-points). Used in cellFeatures class. - -SourceFiles - cellAddressing.C - -\*---------------------------------------------------------------------------*/ - -#ifndef cellAddressing_H -#define cellAddressing_H - -#include "cellModel.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes -class cellModel; - -/*---------------------------------------------------------------------------*\ - Class cellAddressing Declaration -\*---------------------------------------------------------------------------*/ - -class cellAddressing -{ - // Private data - - faceList modelFaces_; - - edgeList edges_; - - labelListList edgeFaces_; - - labelListList faceEdges_; - - labelListList pointEdges_; - - // Private static functions - - static labelList makeIdentity(const label size); - - // Private static functions - - label findEdge(const edge& e); - -public: - - // Constructors - - //- Construct from components - cellAddressing(const cellModel& model); - - // Destructor - - // Member Functions - - // Access - - inline const faceList& modelFaces() const - { - return modelFaces_; - } - - inline const edgeList& edges() const - { - return edges_; - } - - inline const labelListList& edgeFaces() const - { - return edgeFaces_; - } - - inline const labelListList& faceEdges() const - { - return faceEdges_; - } - - inline const labelListList& pointEdges() const - { - return pointEdges_; - } - - - // Check - - // Edit - - // Write - - - // Member Operators - - // Friend Functions - - // Friend Operators - - // IOstream Operators - -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -//#include "cellAddressingI.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurface.C b/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurface.C deleted file mode 100644 index 6d0017adcfaa690456629d5d166552379851dc36..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurface.C +++ /dev/null @@ -1,1068 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Description - - Cuts tets resulting from tet decomposition. See Userguide, - 'Cell shapes in the traditional mesh description' - - f = face - fp = vertex of face - f[fp] = mesh vertex - fc = face centre - cc = cell centre - - Caters for both types of decomposition: - - cellDecomp: introduction of cellCentre. Faces get triangulated by - introducing diagonals. - - faceDecomp: introduction of cellCentre and faceCentres. - - - Tet given by its 4 vertices: - - faceDecomp cellDecomp - ---------- ---------- - 0: f[fp] f[fp] - 1: fc f[0] - 2: f[fp+1] f[fp+1] - 3: cc cc - - Tet given by its 6 edges: - - faceDecomp cellDecomp - ---------- ---------- - 0: fp-fc (face edge) fp-f0 (diagonal edge) - 1: fp-fp+1 ('real', mesh edge) fp-fp+1 (mesh edge) - 2: fp-cc (pyramid edge) fp-cc (pyramid edge) - 3: fc-cc (centre edge) f0-cc (pyramid edge) - 4: fc-fp+1 (face edge) f0-fp+1 (diagonal edge) - 5: fp+1 - cc (pyramid edge) fp+1-cc (pyramid edge) - -\*---------------------------------------------------------------------------*/ - -#include "meshCutSurface.H" -#include "faceDecompCuts.H" -#include "cellDecompCuts.H" -#include "primitiveMesh.H" -#include "cellAddressing.H" -#include "Map.H" -#include "boolList.H" -#include "cellModeller.H" -#include "orientedSurface.H" -#include "boundBox.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -// Symbolic names of tet vertices in faceDecomposition (replace FC by f[0] for -// cellDecomposition) -Foam::label Foam::meshCutSurface::FP0 = 0; -Foam::label Foam::meshCutSurface::FC = 1; -Foam::label Foam::meshCutSurface::FP1 = 2; -Foam::label Foam::meshCutSurface::CC = 3; - -// Symbolic names of tet edges in faceDecomposition -Foam::label Foam::meshCutSurface::FP0_FC = 0; -Foam::label Foam::meshCutSurface::FP0_FP1 = 1; -Foam::label Foam::meshCutSurface::FP0_CC = 2; -Foam::label Foam::meshCutSurface::FC_CC = 3; -Foam::label Foam::meshCutSurface::FC_FP1 = 4; -Foam::label Foam::meshCutSurface::FP1_CC = 5; - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -template<class T, class Hash> -Foam::label Foam::meshCutSurface::find -( - const HashTable<label, T, Hash>& table, - const T& key -) -{ - typename HashTable<label, T, Hash>::const_iterator iter = table.find(key); - - if (iter == table.end()) - { - return -1; - } - else - { - return iter(); - } -} - - -Foam::label Foam::meshCutSurface::findEdge -( - const edgeList& edges, - const labelList& edgeLabels, - const edge& e -) -{ - forAll(edgeLabels, edgeLabelI) - { - label edgeI = edgeLabels[edgeLabelI]; - - if (edges[edgeI] == e) - { - return edgeI; - } - } - - forAll(edgeLabels, edgeLabelI) - { - label edgeI = edgeLabels[edgeLabelI]; - - Pout<< "edge:" << edgeI << " vertices:" << edges[edgeI] << endl; - } - - FatalErrorIn - ( - "meshCutSurface::findEdge(const edgeList&, const labelList&, const edge&)" - ) << "Can not find edge " << e << " among edge labels " << edgeLabels - << exit(FatalError); - - return -1; -} - - -Foam::label Foam::meshCutSurface::findCutEdge -( - const cellAddressing& model, - const labelList& tetEdgeCuts, - const label faceI, - const label cutEdgeI -) -{ - const labelList& fEdges = model.faceEdges()[faceI]; - - forAll(fEdges, fEdgeI) - { - label edgeI = fEdges[fEdgeI]; - - if ((edgeI != cutEdgeI) && (tetEdgeCuts[edgeI] != -1)) - { - return edgeI; - } - } - - FatalErrorIn - ( - "meshCutSurface::findCutEdge(const cellAddressing&, const labelList&" - ", const label, const label)" - ) << "Problem: face " << faceI << " consists of edges " << fEdges - << " of which only one is cut" - << exit(FatalError); - - return -1; -} - - -Foam::triSurface Foam::meshCutSurface::removeDuplicates(const triSurface& surf) -{ - boolList includeTri(surf.size(), true); - - forAll(surf, triI) - { - const labelledTri& t = surf[triI]; - - const labelList& nbs = surf.faceFaces()[triI]; - - forAll(nbs, nbI) - { - if (nbs[nbI] <= triI) - { - // lower numbered faces already checked - continue; - } - - const labelledTri& nb = surf[nbs[nbI]]; - - if - ( - ((t[0] == nb[0]) || (t[0] == nb[1]) || (t[0] == nb[2])) - && ((t[1] == nb[0]) || (t[1] == nb[1]) || (t[1] == nb[2])) - && ((t[2] == nb[0]) || (t[2] == nb[1]) || (t[2] == nb[2])) - ) - { - includeTri[nbs[nbI]] = false; - } - } - } - - labelList pointMap, faceMap; - - return surf.subsetMesh(includeTri, pointMap, faceMap); -} - - -void Foam::meshCutSurface::cutTet -( - const cellAddressing& model, - const labelList& tetVertCuts, - const labelList& tetEdgeCuts, - const label cellI, - - List<labelledTri>& tris, - label& triI -) -{ - label nVertCuts = 0; - forAll(tetVertCuts, tetVertI) - { - if (tetVertCuts[tetVertI] != -1) - { - nVertCuts++; - } - } - - label nEdgeCuts = 0; - forAll(tetEdgeCuts, tetEdgeI) - { - if (tetEdgeCuts[tetEdgeI] != -1) - { - nEdgeCuts++; - } - } - - if (nVertCuts != 0) - { - // Cuts through vertices as well as edges. Handle separately. - cutTetThroughVerts - ( - model, - tetVertCuts, - nVertCuts, - tetEdgeCuts, - nEdgeCuts, - cellI, - - tris, - triI - ); - } - else - { - // Cuts through edges. 'Simple' case. - cutTetThroughEdges - ( - model, - tetEdgeCuts, - nEdgeCuts, - cellI, - - tris, - triI - ); - } -} - - -void Foam::meshCutSurface::cutTetThroughEdges -( - const cellAddressing& model, - const labelList& tetEdgeCuts, - const label nEdgeCuts, - const label cellI, - - List<labelledTri>& tris, - label& triI -) -{ - if (nEdgeCuts < 3) - { - //Pout<< "Tet only cut by " << nEdgeCuts << endl; - } - else if (nEdgeCuts == 3) - { - label triVertI = 0; - - labelledTri& tri = tris[triI++]; - - forAll(tetEdgeCuts, tetEdgeI) - { - if (tetEdgeCuts[tetEdgeI] != -1) - { - tri[triVertI++] = tetEdgeCuts[tetEdgeI]; - } - } - - //Pout<< "Triangle formed by vertices:" << tri << endl; - tri.region() = cellI; - } - else if (nEdgeCuts == 4) - { - // Find ordering of vertices - - label edge0 = -1; - forAll(tetEdgeCuts, tetEdgeI) - { - if (tetEdgeCuts[tetEdgeI] != -1) - { - edge0 = tetEdgeI; - - break; - } - } - - const labelList& nbFaces = model.edgeFaces()[edge0]; - - // Find cut edges among faces using edge0 - label edge1 = findCutEdge(model, tetEdgeCuts, nbFaces[0], edge0); - label edge2 = findCutEdge(model, tetEdgeCuts, nbFaces[1], edge0); - - // Find the other edge - label edge3 = -1; - forAll(tetEdgeCuts, tetEdgeI) - { - if - ( - (tetEdgeCuts[tetEdgeI] != -1) - && (tetEdgeI != edge0) - && (tetEdgeI != edge1) - && (tetEdgeI != edge2) - ) - { - // Fourth cut edge - edge3 = tetEdgeI; - - break; - } - } - - //Pout<< "Triangle 1 formed by vertices:" << tri1 << endl; - tris[triI++] = - labelledTri - ( - tetEdgeCuts[edge0], - tetEdgeCuts[edge1], - tetEdgeCuts[edge2], - cellI - ); - - //Pout<< "Triangle 2 formed by vertices:" << tri2 << endl; - tris[triI++] = - labelledTri - ( - tetEdgeCuts[edge1], - tetEdgeCuts[edge2], - tetEdgeCuts[edge3], - cellI - ); - - } - else - { - //FatalErrorIn("meshCutSurface::cutTet") - // << "Cannot handle tets with " << nEdgeCuts << " edges cut" - // << exit(FatalError); - } -} - - -void Foam::meshCutSurface::cutTetThroughVerts -( - const cellAddressing& model, - const labelList& tetVertCuts, - const label nVertCuts, - const labelList& tetEdgeCuts, - const label nEdgeCuts, - const label cellI, - - List<labelledTri>& tris, - label& triI -) -{ - if (nVertCuts == 1) - { - if (nEdgeCuts == 1) - { - return; - } - else if (nEdgeCuts != 2) - { - // Illegal - return; - } - else // nEdgeCuts == 2 - { - // Find cut-edges not connected to cut-vertex. - - // Triangle to fill. - labelledTri tri; - label triVertI = 0; - - // Get label of vertex cut - label cutVertI = -1; - - forAll(tetVertCuts, tetVertI) - { - if (tetVertCuts[tetVertI] != -1) - { - tri[triVertI++] = tetVertCuts[tetVertI]; - - cutVertI = tetVertI; - - break; - } - } - - // Mark all edges connected to the cut vertex - boolList connectedEdge(6, false); - - const labelList& vEdges = model.pointEdges()[cutVertI]; - - forAll(vEdges, vEdgeI) - { - connectedEdge[vEdges[vEdgeI]] = true; - } - - // Go through all edges and collect vertices from cut edges. - forAll(tetEdgeCuts, tetEdgeI) - { - if (tetEdgeCuts[tetEdgeI] != -1) - { - // Cut edge. - - if (connectedEdge[tetEdgeI]) - { - // Illegal - return; - } - else - { - // Store edge cut vertex - tri[triVertI++] = tetEdgeCuts[tetEdgeI]; - } - } - } - - - if (triVertI == 3) - { - // Collected all triVerts. - tri.region() = cellI; - tris[triI++] = tri; - } - else - { - // Illegal. - } - } - } - else if (nVertCuts == 2) - { - if (nEdgeCuts == 1) - { - // Only interesting if the cut edge is the one opposite the - // edge connecting the two cut vertices. - - // Triangle to fill. - labelledTri tri; - label triVertI = 0; - - // Get label of edge cut - label edgeI = -1; - - forAll(tetEdgeCuts, tetEdgeI) - { - if (tetEdgeCuts[tetEdgeI] != -1) - { - tri[triVertI++] = tetEdgeCuts[tetEdgeI]; - - edgeI = tetEdgeI; - - break; - } - } - - - // Find cut vertices which are not on edge edgeI. - const edge& e = model.edges()[edgeI]; - - forAll(tetVertCuts, tetVertI) - { - if - ( - (tetVertCuts[tetVertI] != -1) - && (tetVertI != e.start()) - && (tetVertI != e.end()) - ) - { - tri[triVertI++] = tetVertCuts[tetVertI]; - } - } - - - if (triVertI == 3) - { - // Collected all triVerts. - //Pout<< "Triangle formed by vertices:" << tri << endl; - tri.region() = cellI; - tris[triI++] = tri; - } - else - { - // Illegal. - } - } - else // nEdgeCuts != 1 - { - // Illegal if 2 or more edges cut. - } - } - else if (nVertCuts == 3) - { - if (nEdgeCuts == 0) - { - // Cut coincides with face. Form triangle out of cut vertices. - - // Triangle to fill. - labelledTri tri; - label triVertI = 0; - - forAll(tetVertCuts, tetVertI) - { - if (tetVertCuts[tetVertI] != -1) - { - tri[triVertI++] = tetVertCuts[tetVertI]; - } - } - tri.region() = cellI; - tris[triI++] = tri; - } - else // nEdgeCuts != 0 - { - // Illegal. - } - } - else // nVertCuts <= 0 || > 3 - { - // Illegal. - } -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -// Construct from edge cuts. Is given edges and position on edge. -// - cuts.cells(): list of cells affected in any way by cutting. -// - cuts.meshVerts(): labels of mesh vertices that are cut -// (cut exactly through point -// - cuts.faceCentres(): face labels whose faceCentres are cut -// - cuts.cellCentres(): cell labels whose cellcentres are cut -// - cuts.meshEdges(): labels of mesh edges that are cut -// - cuts.pyrEdges(): pyramid edges cut -// - cuts.centreEdges(): face-cellCentre edges cut -// - cuts.faceEdges(): face-face edges cut -Foam::meshCutSurface::meshCutSurface(const faceDecompCuts& cuts) -: - triSurface() -{ - const primitiveMesh& mesh = cuts.mesh(); - - // Collect points and create reverse map from edge to new vertex. - pointField points(cuts.size()); - label pointI = 0; - - // - // Cuts through existing points - // - - // Mesh vertices cut - Map<label> meshVertToVert(cuts.meshVerts().size()*2); - - forAll(cuts.meshVerts(), cutVertI) - { - points[pointI++] = mesh.points()[cuts.meshVerts()[cutVertI]]; - - meshVertToVert.insert(cuts.meshVerts()[cutVertI], pointI - 1); - } - - // Face centres cut - Map<label> faceCentreToVert(cuts.meshFaceCentres().size()*2); - - forAll(cuts.meshFaceCentres(), cutFaceCentreI) - { - label faceI = cuts.meshFaceCentres()[cutFaceCentreI]; - - points[pointI++] = mesh.faceCentres()[faceI]; - - faceCentreToVert.insert(faceI, pointI - 1); - } - - // Cell centres cut - Map<label> cellCentreToVert(cuts.meshCellCentres().size()*2); - - forAll(cuts.meshCellCentres(), cutCellCentreI) - { - label cellI = cuts.meshCellCentres()[cutCellCentreI]; - - points[pointI++] = mesh.cellCentres()[cellI]; - - cellCentreToVert.insert(cellI, pointI - 1); - } - - - // - // Cuts through edges - // - - // Cuts through existing (mesh) edges - Map<label> meshEdgeToVert(cuts.meshEdges().size()*2); - - forAll(cuts.meshEdges(), meshCutSurfaceEdgeI) - { - label edgeI = cuts.meshEdges()[meshCutSurfaceEdgeI]; - - const edge& e = mesh.edges()[edgeI]; - - scalar weight = cuts.meshEdgeWeights()[meshCutSurfaceEdgeI]; - - points[pointI++] = - (1-weight)*mesh.points()[e.start()] - + weight*mesh.points()[e.end()]; - - meshEdgeToVert.insert(edgeI, pointI - 1); - } - - // Cuts on tetDecomposition: Pyramid edges - HashTable<label, pyramidEdge, pyramidEdge::pyramidEdgeHash> - pyrEdgeToVert(cuts.pyrEdges().size()*2); - - forAll(cuts.pyrEdges(), edgeI) - { - const pyramidEdge& e = cuts.pyrEdges()[edgeI]; - - points[pointI++] = e.coord(mesh, cuts.pyrEdgeWeights()[edgeI]); - - pyrEdgeToVert.insert(e, pointI - 1); - } - - // Cuts on tetDecomposition: Face-Cell centre edges - HashTable<label, centreEdge, centreEdge::centreEdgeHash> - centreEdgeToVert(cuts.centreEdges().size()*2); - - forAll(cuts.centreEdges(), edgeI) - { - const centreEdge& e = cuts.centreEdges()[edgeI]; - - points[pointI++] = e.coord(mesh, cuts.centreEdgeWeights()[edgeI]); - - centreEdgeToVert.insert(e, pointI - 1); - } - - // Cuts on tetDecomposition: Face edges - HashTable<label, faceEdge, faceEdge::faceEdgeHash> - faceEdgeToVert(cuts.faceEdges().size()*2); - - forAll(cuts.faceEdges(), edgeI) - { - const faceEdge& e = cuts.faceEdges()[edgeI]; - - points[pointI++] = e.coord(mesh, cuts.faceEdgeWeights()[edgeI]); - - faceEdgeToVert.insert(e, pointI - 1); - } - - if (pointI != points.size()) - { - FatalErrorIn - ( - "meshCutSurface::meshCutSurface(const faceDecompCuts& cuts)" - ) << "pointI:" << pointI << " points.size():" << points.size() - << abort(FatalError); - } - - - // Generate tet model + additional addressing - const cellModel& tet = *(cellModeller::lookup("tet")); - - cellAddressing tetPlus(tet); - - // Tet edges cut. -1 if edge not cut, vertex label otherwise. - labelList tetEdgeCuts(6); - - // Tet vertices cut. -1 if edge not cut, vertex label otherwise. - labelList tetVertCuts(4); - - List<labelledTri> tris(2*cuts.size()); - label triI = 0; - - forAll(cuts.cells(), cutCellI) - { - label cellI = cuts.cells()[cutCellI]; - - tetVertCuts[CC] = find(cellCentreToVert, cellI); - - const cell& cellFaces = mesh.cells()[cellI]; - - forAll(cellFaces, cellFaceI) - { - label faceI = cellFaces[cellFaceI]; - - tetVertCuts[FC] = find(faceCentreToVert, faceI); - - bool isOwner = (mesh.faceOwner()[faceI] == cellI); - - tetEdgeCuts[FC_CC] = find - ( - centreEdgeToVert, - centreEdge(faceI, isOwner) - ); - - const face& f = mesh.faces()[faceI]; - - forAll(f, fp) - { - label fp1 = (fp+1) % f.size(); - - tetVertCuts[FP0] = find(meshVertToVert, f[fp]); - tetVertCuts[FP1] = find(meshVertToVert, f[fp1]); - - - tetEdgeCuts[FP0_FC] = find - ( - faceEdgeToVert, - faceEdge(faceI, fp) - ); - - - label edgeI = findEdge - ( - mesh.edges(), - mesh.faceEdges()[faceI], - edge(f[fp], f[fp1]) - ); - tetEdgeCuts[FP0_FP1] = find - ( - meshEdgeToVert, - edgeI - ); - - tetEdgeCuts[FP0_CC] = find - ( - pyrEdgeToVert, - pyramidEdge(f[fp], cellI) - ); - - tetEdgeCuts[FC_FP1] = find - ( - faceEdgeToVert, - faceEdge(faceI, fp1) - ); - - tetEdgeCuts[FP1_CC] = find - ( - pyrEdgeToVert, - pyramidEdge(f[fp1], cellI) - ); - - - cutTet - ( - tetPlus, - tetVertCuts, - tetEdgeCuts, - cellI, - - tris, - triI - ); - } - } - } - - tris.setSize(triI); - - // Get a point outside all of the surface (over all processors) - boundBox bb(points, true); - - point outsidePt(1.5*bb.max() - 0.5*bb.min()); - - // Orient the surface after removing duplicate faces. - triSurface::operator= - ( - orientedSurface - ( - removeDuplicates - ( - triSurface(tris, points) - ), - outsidePt - ) - ); -} - - -// Construct from edge cuts. Is given edges and position on edge. -// - cuts.cells(): list of cells affected in any way by cutting. -// - cuts.meshVerts(): labels of mesh vertices that are cut (cut exactly through -// point -// - cuts.cellCentres(): cell labels whose cellcentres are cut -// - cuts.meshEdges(): labels of mesh edges that are cut -// - cuts.pyrEdges(): pyramid edges cut -// - cuts.diagEdges(): diagonal edges cut -Foam::meshCutSurface::meshCutSurface(const cellDecompCuts& cuts) -: - triSurface() -{ - const primitiveMesh& mesh = cuts.mesh(); - - // Collect points and create reverse map from edge to new vertex. - pointField points(cuts.size()); - label pointI = 0; - - // - // Cuts through existing points - // - - // Mesh vertices cut - Map<label> meshVertToVert(cuts.meshVerts().size()*2); - - forAll(cuts.meshVerts(), cutVertI) - { - points[pointI++] = mesh.points()[cuts.meshVerts()[cutVertI]]; - - meshVertToVert.insert(cuts.meshVerts()[cutVertI], pointI - 1); - } - - // Cell centres cut - Map<label> cellCentreToVert(cuts.meshCellCentres().size()*2); - - forAll(cuts.meshCellCentres(), cutCellCentreI) - { - label cellI = cuts.meshCellCentres()[cutCellCentreI]; - - points[pointI++] = mesh.cellCentres()[cellI]; - - cellCentreToVert.insert(cellI, pointI - 1); - } - - - // - // Cuts through edges - // - - // Cuts through existing (mesh) edges - Map<label> meshEdgeToVert(cuts.meshEdges().size()*2); - - forAll(cuts.meshEdges(), meshCutSurfaceEdgeI) - { - label edgeI = cuts.meshEdges()[meshCutSurfaceEdgeI]; - - const edge& e = mesh.edges()[edgeI]; - - scalar weight = cuts.meshEdgeWeights()[meshCutSurfaceEdgeI]; - - points[pointI++] = - (1-weight)*mesh.points()[e.start()] - + weight*mesh.points()[e.end()]; - - meshEdgeToVert.insert(edgeI, pointI - 1); - } - - // Cuts on tetDecomposition: Pyramid edges - HashTable<label, pyramidEdge, pyramidEdge::pyramidEdgeHash> - pyrEdgeToVert(cuts.pyrEdges().size()*2); - - forAll(cuts.pyrEdges(), edgeI) - { - const pyramidEdge& e = cuts.pyrEdges()[edgeI]; - - points[pointI++] = e.coord(mesh, cuts.pyrEdgeWeights()[edgeI]); - - pyrEdgeToVert.insert(e, pointI - 1); - } - - // Cuts on tetDecomposition: diagonal edges - HashTable<label, diagonalEdge, diagonalEdge::diagonalEdgeHash> - diagEdgeToVert(cuts.diagEdges().size()*2); - - forAll(cuts.diagEdges(), edgeI) - { - const diagonalEdge& e = cuts.diagEdges()[edgeI]; - - points[pointI++] = e.coord(mesh, cuts.diagEdgeWeights()[edgeI]); - - diagEdgeToVert.insert(e, pointI - 1); - } - - - - // Generate tet model + additional addressing - const cellModel& tet = *(cellModeller::lookup("tet")); - - cellAddressing tetPlus(tet); - - // Tet edges cut. -1 if edge not cut, vertex label otherwise. - labelList tetEdgeCuts(6); - - // Tet vertices cut. -1 if edge not cut, vertex label otherwise. - labelList tetVertCuts(4); - - List<labelledTri> tris(2*cuts.size()); - label triI = 0; - - forAll(cuts.cells(), cutCellI) - { - label cellI = cuts.cells()[cutCellI]; - - tetVertCuts[CC] = find(cellCentreToVert, cellI); - - const cell& cellFaces = mesh.cells()[cellI]; - - forAll(cellFaces, cellFaceI) - { - label faceI = cellFaces[cellFaceI]; - - const face& f = mesh.faces()[faceI]; - - // Handle f0 separately - tetVertCuts[FC] = find(meshVertToVert, f[0]); - tetEdgeCuts[FC_CC] = find - ( - pyrEdgeToVert, - pyramidEdge(f[0], cellI) - ); - - for (label fp = 1; fp < f.size()-1; fp++) - { - label fp1 = fp+1; - - tetVertCuts[FP0] = find(meshVertToVert, f[fp]); - tetVertCuts[FP1] = find(meshVertToVert, f[fp1]); - - tetEdgeCuts[FP0_FP1] = find - ( - meshEdgeToVert, - findEdge - ( - mesh.edges(), - mesh.faceEdges()[faceI], - edge(f[fp], f[fp1]) - ) - ); - - if (fp == 1) - { - // f[1] to f[0] is normal mesh edge - tetEdgeCuts[FP0_FC] = find - ( - meshEdgeToVert, - findEdge - ( - mesh.edges(), - mesh.faceEdges()[faceI], - edge(f[0], f[fp]) - ) - ); - } - else - { - // Diagonal cut - tetEdgeCuts[FP0_FC] = find - ( - diagEdgeToVert, - diagonalEdge(faceI, 0, fp) - ); - } - - if (fp1 == f.size()-1) - { - // Last face vertex to f[0] -> normal mesh edge - tetEdgeCuts[FC_FP1] = find - ( - meshEdgeToVert, - findEdge - ( - mesh.edges(), - mesh.faceEdges()[faceI], - edge(f[0], f[fp1]) - ) - ); - } - else - { - tetEdgeCuts[FC_FP1] = find - ( - diagEdgeToVert, - diagonalEdge(faceI, 0, fp1) - ); - } - - tetEdgeCuts[FP0_CC] = find - ( - pyrEdgeToVert, - pyramidEdge(f[fp], cellI) - ); - - tetEdgeCuts[FP1_CC] = find - ( - pyrEdgeToVert, - pyramidEdge(f[fp1], cellI) - ); - - - cutTet - ( - tetPlus, - tetVertCuts, - tetEdgeCuts, - cellI, - - tris, - triI - ); - } - } - } - - tris.setSize(triI); - - // Get a point outside all of the surface (over all processors) - boundBox bb(points, true); - - point outsidePt(1.5*bb.max() - 0.5*bb.min()); - - // Orient the surface after removing duplicate faces. - triSurface::operator= - ( - orientedSurface - ( - removeDuplicates - ( - triSurface(tris, points) - ), - outsidePt - ) - ); -} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurface.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurface.H deleted file mode 100644 index c5d92529ff6f9e3e6d52170c6519ca6a0860496e..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurface.H +++ /dev/null @@ -1,209 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::meshCutSurface - -Description - Description of surface resulting from cuts through mesh. Region on - labelledTris is the cell the triangle is resulting from. - -SourceFiles - meshCutSurface.C - -\*---------------------------------------------------------------------------*/ - -#ifndef meshCutSurface_H -#define meshCutSurface_H - -#include "labelList.H" -#include "HashTable.H" -#include "pointField.H" -#include "triSurface.H" -#include "boolList.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes -class primitiveMesh; -class cellAddressing; -class faceDecompCuts; -class cellDecompCuts; - -/*---------------------------------------------------------------------------*\ - Class meshCutSurface Declaration -\*---------------------------------------------------------------------------*/ - -class meshCutSurface -: - public triSurface -{ - - // Static data members - - //- Symbolic names of tet vertices - static label FP0; - static label FC; - static label FP1; - static label CC; - - //- Symbolic names of tet edges - static label FP0_FC; - static label FP0_FP1; - static label FP0_CC; - static label FC_CC; - static label FC_FP1; - static label FP1_CC; - - - // Private Static Functions - - //- Find key in hashtable. Returns value or -1 if not found - template<class T, class Hash> - static label find - ( - const HashTable<label, T, Hash>& table, - const T& key - ); - - //- Find edge in subset edgeLabels of edges. - static label findEdge - ( - const edgeList& edges, - const labelList& edgeLabels, - const edge& e - ); - - //- Find the other (i.e. not cutEdgeI) cut-edge on faceI of cell model - static label findCutEdge - ( - const cellAddressing& model, - const labelList& tetEdgeCuts, - const label faceI, - const label cutEdgeI - ); - - - // Private Member Functions - - //- Remove duplicate faces (resulting from 3 meshVerts cut - is seen - // from both sides) - static triSurface removeDuplicates(const triSurface& surf); - - -public: - - //- Generic tet cutting routine. - // Generate triangles given the cut vertices and edges of the tet. - // (Triangles can be incorrectly oriented) - void cutTet - ( - const cellAddressing& model, - const labelList& tetVertCuts, - const labelList& tetEdgeCuts, - const label cellI, - List<labelledTri>& tris, - label& triI - ); - - //- Tet cutting when cuts are only through edges. - void cutTetThroughEdges - ( - const cellAddressing& model, - const labelList& tetEdgeCuts, - const label nEdgeCuts, - const label cellI, - List<labelledTri>& tris, - label& triI - ); - - //- Tet cutting when cuts are through corner vertices (and possibly - // edges) - void cutTetThroughVerts - ( - const cellAddressing& model, - const labelList& tetVertCuts, - const label nVertCuts, - const labelList& tetEdgeCuts, - const label nEdgeCuts, - const label cellI, - List<labelledTri>& tris, - label& triI - ); - - - // Static functions - - //- Interpolate onto points given edge weights, cell centre values, - // face centre values and vertex values - template <class T> - static tmp<Field<T> > interpolate - ( - const faceDecompCuts& cuts, - const Field<T>& vField, - const Field<T>& fField, - const Field<T>& pField - ); - - //- Interpolate onto points given edge weights, cell centre values, - // and vertex values - template <class T> - static tmp<Field<T> > interpolate - ( - const cellDecompCuts& cuts, - const Field<T>& vField, - const Field<T>& pField - ); - - - // Constructors - - //- Construct from cuts on faceDecomp tet decomposition - // (introduces cellcentres and facecentres) - meshCutSurface(const faceDecompCuts& cuts); - - //- Construct from cuts on cellDecomp tet decomposition - // (introduces cellcentres but not facecentres) - meshCutSurface(const cellDecompCuts& cuts); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#ifdef NoRepository -# include "meshCutSurfaceInterpolate.C" -#endif - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurfaceInterpolate.C b/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurfaceInterpolate.C deleted file mode 100644 index 37327e04670b737aebea0bfc8db1677a3e6afc88..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/meshCut/meshCutSurfaceInterpolate.C +++ /dev/null @@ -1,174 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -\*---------------------------------------------------------------------------*/ - -#include "meshCutSurface.H" -#include "pyramidEdge.H" -#include "faceEdge.H" -#include "centreEdge.H" -#include "diagonalEdge.H" -#include "faceDecompCuts.H" -#include "cellDecompCuts.H" - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -template <class T> -Foam::tmp<Foam::Field<T> > Foam::meshCutSurface::interpolate -( - const faceDecompCuts& cuts, - const Field<T>& vField, - const Field<T>& fField, - const Field<T>& pField -) -{ - const primitiveMesh& mesh = cuts.mesh(); - - tmp<Field<T> > tvals(new Field<T>(cuts.size())); - Field<T>& vals = tvals(); - - label pointI = 0; - - forAll(cuts.meshVerts(), cutVertI) - { - vals[pointI++] = pField[cuts.meshVerts()[cutVertI]]; - } - - forAll(cuts.meshFaceCentres(), cutFaceCentreI) - { - label faceI = cuts.meshFaceCentres()[cutFaceCentreI]; - - vals[pointI++] = fField[faceI]; - } - - forAll(cuts.meshCellCentres(), cutCellCentreI) - { - label cellI = cuts.meshCellCentres()[cutCellCentreI]; - - vals[pointI++] = vField[cellI]; - } - - forAll(cuts.meshEdges(), meshCutSurfaceEdgeI) - { - label edgeI = cuts.meshEdges()[meshCutSurfaceEdgeI]; - - const edge& e = mesh.edges()[edgeI]; - - scalar weight = cuts.meshEdgeWeights()[meshCutSurfaceEdgeI]; - - vals[pointI++] = weight*pField[e.start()] + (1-weight)*pField[e.end()]; - } - - forAll(cuts.pyrEdges(), edgeI) - { - const pyramidEdge& e = cuts.pyrEdges()[edgeI]; - - scalar weight = cuts.pyrEdgeWeights()[edgeI]; - - vals[pointI++] = e.interpolate(vField, pField, weight); - } - - forAll(cuts.centreEdges(), edgeI) - { - const centreEdge& e = cuts.centreEdges()[edgeI]; - - scalar weight = cuts.centreEdgeWeights()[edgeI]; - - vals[pointI++] = e.interpolate(mesh, vField, fField, weight); - } - - forAll(cuts.faceEdges(), edgeI) - { - const faceEdge& e = cuts.faceEdges()[edgeI]; - - scalar weight = cuts.faceEdgeWeights()[edgeI]; - - vals[pointI++] = e.interpolate(mesh, fField, pField, weight); - } - - return tvals; -} - - -template <class T> -Foam::tmp<Foam::Field<T> > Foam::meshCutSurface::interpolate -( - const cellDecompCuts& cuts, - const Field<T>& vField, - const Field<T>& pField -) -{ - const primitiveMesh& mesh = cuts.mesh(); - - tmp<Field<T> > tvals(new Field<T>(cuts.size())); - Field<T>& vals = tvals(); - - label pointI = 0; - - forAll(cuts.meshVerts(), cutVertI) - { - vals[pointI++] = pField[cuts.meshVerts()[cutVertI]]; - } - - forAll(cuts.meshCellCentres(), cutCellCentreI) - { - label cellI = cuts.meshCellCentres()[cutCellCentreI]; - - vals[pointI++] = vField[cellI]; - } - - forAll(cuts.meshEdges(), meshCutSurfaceEdgeI) - { - label edgeI = cuts.meshEdges()[meshCutSurfaceEdgeI]; - - const edge& e = mesh.edges()[edgeI]; - - scalar weight = cuts.meshEdgeWeights()[meshCutSurfaceEdgeI]; - - vals[pointI++] = weight*pField[e.start()] + (1-weight)*pField[e.end()]; - } - - forAll(cuts.pyrEdges(), edgeI) - { - const pyramidEdge& e = cuts.pyrEdges()[edgeI]; - - scalar weight = cuts.pyrEdgeWeights()[edgeI]; - - vals[pointI++] = e.interpolate(vField, pField, weight); - } - - forAll(cuts.diagEdges(), edgeI) - { - const diagonalEdge& e = cuts.diagEdges()[edgeI]; - - scalar weight = cuts.diagEdgeWeights()[edgeI]; - - vals[pointI++] = e.interpolate(mesh, pField, weight); - } - - return tvals; -} - - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/centreEdge.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/centreEdge.H deleted file mode 100644 index a7646089e91908c2e754c6ba6508c473039c3184..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/centreEdge.H +++ /dev/null @@ -1,187 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::centreEdge - -Description - Implicit description of faceCentre to cellCentre edge - (from tet decomposition). - - Stores - - face label - - whether in owner or neighbour of face - -SourceFiles - -\*---------------------------------------------------------------------------*/ - -#ifndef centreEdge_H -#define centreEdge_H - -#include "label.H" -#include "primitiveMesh.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes - -/*---------------------------------------------------------------------------*\ - Class centreEdge Declaration -\*---------------------------------------------------------------------------*/ - -class centreEdge -{ - // Private data - - //- face label - label faceLabel_; - - //- is in owner or neighbour of face - bool inOwner_; - -public: - - // Public classes - - //- Hash function - class centreEdgeHash - { - - public: - - centreEdgeHash() - {} - - //- simple hashing function of labels - label operator()(const centreEdge& ce, const label tableSize) const - { - return (ce.faceLabel() + ce.inOwner()) % tableSize; - } - }; - - - // Constructors - - //- Construct null - inline centreEdge() - : - faceLabel_(-1), - inOwner_(-1) - {} - - //- Construct from components - inline centreEdge(const label faceLabel, const bool inOwner) - : - faceLabel_(faceLabel), - inOwner_(inOwner) - {} - - // Member Functions - - label faceLabel() const - { - return faceLabel_; - } - - bool inOwner() const - { - return inOwner_; - } - - template <class T> - T interpolate - ( - const primitiveMesh& mesh, - const Field<T>& cellField, - const Field<T>& faceField, - const scalar weight - ) const - { - label cellI; - if (inOwner_) - { - cellI = mesh.faceOwner()[faceLabel_]; - } - else - { - cellI = mesh.faceNeighbour()[faceLabel_]; - } - - return - (1-weight)*faceField[faceLabel_] - + weight*cellField[cellI]; - } - - point coord(const primitiveMesh& mesh, const scalar weight) const - { - return interpolate - ( - mesh, - mesh.cellCentres(), - mesh.faceCentres(), - weight - ); - } - - - // Member Operators - - bool operator==(const centreEdge& ce) const - { - return - (faceLabel() == ce.faceLabel()) - && (inOwner() == ce.inOwner()); - } - - - // IOstream Operators - - inline friend Ostream& operator<<(Ostream& os, const centreEdge& ce) - { - os << token::BEGIN_LIST - << ce.faceLabel_ << token::SPACE - << ce.inOwner_ - << token::END_LIST; - - // Check state of Ostream - os.check("Ostream& operator<<(Ostream&, const centreEdge&)"); - - return os; - } -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/diagonalEdge.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/diagonalEdge.H deleted file mode 100644 index fed2e1d9f5fee243a0e3ae6296fb47025e33f8b5..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/diagonalEdge.H +++ /dev/null @@ -1,199 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::diagonalEdge - -Description - Implicit description of diagonal edge (from tet decomposition). - Diangonal edge is edge between two vertices of same face. - - Stores - - face label - - indices in face - -SourceFiles - -\*---------------------------------------------------------------------------*/ - -#ifndef diagonalEdge_H -#define diagonalEdge_H - -#include "label.H" -#include "primitiveMesh.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class diagonalEdge Declaration -\*---------------------------------------------------------------------------*/ - -class diagonalEdge -{ - // Private data - - //- face label - label faceLabel_; - - //- index in face - label faceIndex0_; - - //- index in face - label faceIndex1_; - -public: - - // Public classes - - //- Hash function - class diagonalEdgeHash - { - - public: - - diagonalEdgeHash() - {} - - //- simple hashing function of labels. Commutative on face indices - // since edge is same. - label operator()(const diagonalEdge& de, const label tableSize) - const - { - // Note: mag since multiply might overflow and produce - // negative numbers - return - mag - ( - de.faceLabel() - + (de.faceLabel()*de.faceLabel()) - + de.faceIndex0() * de.faceIndex1() - + (de.faceIndex0() + de.faceIndex1()) - ) - % tableSize; - } - }; - - - // Constructors - - //- Construct null - inline diagonalEdge() - : - faceLabel_(-1), - faceIndex0_(-1), - faceIndex1_(-1) - {} - - //- Construct from components - inline diagonalEdge - ( - const label faceLabel, - const label faceIndex0, - const label faceIndex1 - ) - : - faceLabel_(faceLabel), - faceIndex0_(faceIndex0), - faceIndex1_(faceIndex1) - {} - - // Member Functions - - label faceLabel() const - { - return faceLabel_; - } - - label faceIndex0() const - { - return faceIndex0_; - } - - label faceIndex1() const - { - return faceIndex1_; - } - - template <class T> - T interpolate - ( - const primitiveMesh& mesh, - const Field<T>& vertField, - const scalar weight - ) const - { - const face& f = mesh.faces()[faceLabel_]; - - return - (1-weight)*vertField[f[faceIndex0_]] - + weight*vertField[f[faceIndex1_]]; - } - - point coord(const primitiveMesh& mesh, const scalar weight) const - { - return interpolate(mesh, mesh.points(), weight); - } - - - // Member Operators - - bool operator==(const diagonalEdge& de) const - { - return - (faceLabel() == de.faceLabel()) - && (faceIndex0() == de.faceIndex0()) - && (faceIndex1() == de.faceIndex1()); - } - - - // IOstream Operators - - inline friend Ostream& operator<<(Ostream& os, const diagonalEdge& de) - { - os << token::BEGIN_LIST - << de.faceLabel_ << token::SPACE - << de.faceIndex0_ << token::SPACE - << de.faceIndex1_ - << token::END_LIST; - - // Check state of Ostream - os.check("Ostream& operator<<(Ostream&, const diagonalEdge&)"); - - return os; - } -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/faceEdge.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/faceEdge.H deleted file mode 100644 index 6676862c7cb3832afb8c4a5bb77a83300c697584..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/faceEdge.H +++ /dev/null @@ -1,186 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::faceEdge - -Description - Implicit description of face edge (from tet decomposition). - Face edge is edge between vertex of face and face centre. - - Stores - - face label - - index in face (gives the mesh vertex) - - -SourceFiles - -\*---------------------------------------------------------------------------*/ - -#ifndef faceEdge_H -#define faceEdge_H - -#include "label.H" -#include "primitiveMesh.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes - -/*---------------------------------------------------------------------------*\ - Class faceEdge Declaration -\*---------------------------------------------------------------------------*/ - -class faceEdge -{ - // Private data - - //- face label - label faceLabel_; - - //- index in face - label faceIndex_; - -public: - - // Public classes - - //- Hash function - class faceEdgeHash - { - - public: - - faceEdgeHash() - {} - - //- simple hashing function of labels - label operator()(const faceEdge& fe, const label tableSize) - const - { - // Note: mag since multiply might overflow and produce - // negative numbers - return - mag - ( - fe.faceLabel() - + (fe.faceLabel() * fe.faceLabel()) - + fe.faceIndex() - + (fe.faceIndex() * fe.faceIndex()) - ) - % tableSize; - } - }; - - - // Constructors - - //- Construct null - inline faceEdge() - : - faceLabel_(-1), - faceIndex_(-1) - {} - - //- Construct from components - inline faceEdge(const label faceLabel, const label faceIndex) - : - faceLabel_(faceLabel), - faceIndex_(faceIndex) - {} - - // Member Functions - - label faceLabel() const - { - return faceLabel_; - } - - label faceIndex() const - { - return faceIndex_; - } - - template <class T> - T interpolate - ( - const primitiveMesh& mesh, - const Field<T>& faceField, - const Field<T>& vertField, - const scalar weight - ) const - { - const face& f = mesh.faces()[faceLabel_]; - - return - (1-weight)*vertField[f[faceIndex_]] - + weight*faceField[faceLabel_]; - } - - - point coord(const primitiveMesh& mesh, const scalar weight) const - { - return interpolate(mesh, mesh.faceCentres(), mesh.points(), weight); - } - - - // Member Operators - - bool operator==(const faceEdge& fe) const - { - return - (faceLabel() == fe.faceLabel()) - && (faceIndex() == fe.faceIndex()); - } - - - // IOstream Operators - - inline friend Ostream& operator<<(Ostream& os, const faceEdge& fe) - { - os << token::BEGIN_LIST - << fe.faceLabel_ << token::SPACE - << fe.faceIndex_ - << token::END_LIST; - - // Check state of Ostream - os.check("Ostream& operator<<(Ostream&, const faceEdge&)"); - - return os; - } -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/pyramidEdge.H b/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/pyramidEdge.H deleted file mode 100644 index b05fc6d2ef1756d93e39735caa0c3aa3ae4145ff..0000000000000000000000000000000000000000 --- a/src/finiteVolume/fvMesh/fvMeshCutSurface/tetEdges/pyramidEdge.H +++ /dev/null @@ -1,185 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software; you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by the - Free Software Foundation; either version 2 of the License, or (at your - option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Class - Foam::pyramidEdge - -Description - Implicit description of pyramid edge (from tet decomposition). - Pyramid edge is edge between mesh vertex and cell centre. - - Stores - - vertex label - - cell label - -SourceFiles - -\*---------------------------------------------------------------------------*/ - -#ifndef pyramidEdge_H -#define pyramidEdge_H - -#include "label.H" -#include "primitiveMesh.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// Forward declaration of classes - -/*---------------------------------------------------------------------------*\ - Class pyramidEdge Declaration -\*---------------------------------------------------------------------------*/ - -class pyramidEdge -{ - // Private data - - //- vertex label - label vertexLabel_; - - //- cell label - label cellLabel_; - -public: - - // Public classes - - //- Hash function - class pyramidEdgeHash - { - - public: - - pyramidEdgeHash() - {} - - //- simple hashing function of labels - label operator()(const pyramidEdge& pe, const label tableSize) const - { - // Note: mag since multiply might overflow and produce - // negative numbers - return - mag - ( - pe.vertexLabel() - + (pe.vertexLabel() * pe.vertexLabel()) - + pe.cellLabel() - + (pe.cellLabel() * pe.cellLabel()) - ) - % tableSize; - } - }; - - - // Constructors - - //- Construct null - inline pyramidEdge() - : - vertexLabel_(-1), - cellLabel_(-1) - {} - - //- Construct from components - inline pyramidEdge - ( - const label vertexLabel, - const label cellLabel - ) - : - vertexLabel_(vertexLabel), - cellLabel_(cellLabel) - {} - - // Member Functions - - label vertexLabel() const - { - return vertexLabel_; - } - - label cellLabel() const - { - return cellLabel_; - } - - template <class T> - T interpolate - ( - const Field<T>& cellField, - const Field<T>& vertField, - const scalar weight - ) const - { - return - (1-weight)*vertField[vertexLabel_] - + weight*cellField[cellLabel_]; - } - - - point coord(const primitiveMesh& mesh, const scalar weight) const - { - return interpolate(mesh.cellCentres(), mesh.points(), weight); - } - - - // Member Operators - - bool operator==(const pyramidEdge& pe) const - { - return - (vertexLabel() == pe.vertexLabel()) - && (cellLabel() == pe.cellLabel()); - } - - - // IOstream Operators - - inline friend Ostream& operator<<(Ostream& os, const pyramidEdge& pe) - { - os << token::BEGIN_LIST - << pe.vertexLabel_ << token::SPACE - << pe.cellLabel_ - << token::END_LIST; - - // Check state of Ostream - os.check("Ostream& operator<<(Ostream&, const pyramidEdge&)"); - - return os; - } -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/sampling/Make/files b/src/sampling/Make/files index f17e7222bde3226b94b774ab82ef2df25c7fad68..991057a157d640b8a0070c58e4e12967a74da0fe 100644 --- a/src/sampling/Make/files +++ b/src/sampling/Make/files @@ -24,6 +24,9 @@ sampledSurface/patch/sampledPatch.C sampledSurface/plane/sampledPlane.C sampledSurface/isoSurface/isoSurface.C sampledSurface/isoSurface/sampledIsoSurface.C +sampledSurface/isoSurface/isoSurfaceCell.C +sampledSurface/isoSurface/sampledIsoSurfaceCell.C +sampledSurface/distanceSurface/distanceSurface.C sampledSurface/sampledSurface/sampledSurface.C sampledSurface/sampledSurfaces/sampledSurfaces.C sampledSurface/sampledSurfacesFunctionObject/sampledSurfacesFunctionObject.C diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C new file mode 100644 index 0000000000000000000000000000000000000000..8c4acea9967b39942148f148308d5ae83c0af641 --- /dev/null +++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C @@ -0,0 +1,342 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "distanceSurface.H" +#include "dictionary.H" +#include "volFields.H" +#include "volPointInterpolation.H" +#include "addToRunTimeSelectionTable.H" +#include "fvMesh.H" +#include "isoSurfaceCell.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(distanceSurface, 0); + addToRunTimeSelectionTable(sampledSurface, distanceSurface, word); +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::distanceSurface::createGeometry() const +{ + // Clear any stored topo + facesPtr_.clear(); + + // Distance to cell centres + scalarField cellDistance(mesh().nCells()); + { + List<pointIndexHit> nearest; + surfPtr_().findNearest + ( + mesh().cellCentres(), + scalarField(mesh().nCells(), GREAT), + nearest + ); + + if (signed_) + { + vectorField normal; + surfPtr_().getNormal(nearest, normal); + + forAll(cellDistance, cellI) + { + vector d(mesh().cellCentres()[cellI]-nearest[cellI].hitPoint()); + + if ((d&normal[cellI]) > 0) + { + cellDistance[cellI] = Foam::mag(d); + } + else + { + cellDistance[cellI] = -Foam::mag(d); + } + } + } + else + { + forAll(cellDistance, cellI) + { + cellDistance[cellI] = Foam::mag + ( + nearest[cellI].hitPoint() + - mesh().cellCentres()[cellI] + ); + } + } + } + + // Distance to points + scalarField pointDistance(mesh().nPoints()); + { + List<pointIndexHit> nearest; + surfPtr_().findNearest + ( + mesh().points(), + scalarField(mesh().nPoints(), GREAT), + nearest + ); + forAll(pointDistance, pointI) + { + pointDistance[pointI] = Foam::mag + ( + nearest[pointI].hitPoint() + - mesh().points()[pointI] + ); + } + } + + //- Direct from cell field and point field. + const isoSurfaceCell iso + ( + mesh(), + cellDistance, + pointDistance, + distance_, + regularise_ + ); + + ////- From point field and interpolated cell. + //scalarField cellAvg(mesh().nCells(), scalar(0.0)); + //labelField nPointCells(mesh().nCells(), 0); + //{ + // for (label pointI = 0; pointI < mesh().nPoints(); pointI++) + // { + // const labelList& pCells = mesh().pointCells(pointI); + // + // forAll(pCells, i) + // { + // label cellI = pCells[i]; + // + // cellAvg[cellI] += pointDistance[pointI]; + // nPointCells[cellI]++; + // } + // } + //} + //forAll(cellAvg, cellI) + //{ + // cellAvg[cellI] /= nPointCells[cellI]; + //} + // + //const isoSurface iso + //( + // mesh(), + // cellAvg, + // pointDistance, + // distance_, + // regularise_ + //); + + + const_cast<distanceSurface&>(*this).triSurface::operator=(iso); + meshCells_ = iso.meshCells(); + + if (debug) + { + print(Pout); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::distanceSurface::distanceSurface +( + const word& name, + const polyMesh& mesh, + const dictionary& dict +) +: + sampledSurface(name, mesh, dict), + surfPtr_ + ( + searchableSurface::New + ( + dict.lookup("surfaceType"), + IOobject + ( + dict.lookupOrDefault("surfaceName", name), // name + mesh.time().constant(), // directory + "triSurface", // instance + mesh.time(), // registry + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + dict + ) + ), + distance_(readScalar(dict.lookup("distance"))), + signed_(readBool(dict.lookup("signed"))), + regularise_(dict.lookupOrDefault("regularise", true)), + zoneName_(word::null), + facesPtr_(NULL), + storedTimeIndex_(-1), + meshCells_(0) +{ +// label zoneId = -1; +// if (dict.readIfPresent("zone", zoneName_)) +// { +// zoneId = mesh.cellZones().findZoneID(zoneName_); +// if (debug && zoneId < 0) +// { +// Info<< "cellZone \"" << zoneName_ +// << "\" not found - using entire mesh" +// << endl; +// } +// } + correct(true); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::distanceSurface::~distanceSurface() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::distanceSurface::correct(const bool meshChanged) +{ + // Only change of mesh changes plane - zone restriction gets lost + if (meshChanged) + { + createGeometry(); + } +} + + +Foam::tmp<Foam::scalarField> +Foam::distanceSurface::sample +( + const volScalarField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::vectorField> +Foam::distanceSurface::sample +( + const volVectorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::sphericalTensorField> +Foam::distanceSurface::sample +( + const volSphericalTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::symmTensorField> +Foam::distanceSurface::sample +( + const volSymmTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::tensorField> +Foam::distanceSurface::sample +( + const volTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::scalarField> +Foam::distanceSurface::interpolate +( + const interpolation<scalar>& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp<Foam::vectorField> +Foam::distanceSurface::interpolate +( + const interpolation<vector>& interpolator +) const +{ + return interpolateField(interpolator); +} + +Foam::tmp<Foam::sphericalTensorField> +Foam::distanceSurface::interpolate +( + const interpolation<sphericalTensor>& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp<Foam::symmTensorField> +Foam::distanceSurface::interpolate +( + const interpolation<symmTensor>& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp<Foam::tensorField> +Foam::distanceSurface::interpolate +( + const interpolation<tensor>& interpolator +) const +{ + return interpolateField(interpolator); +} + + +void Foam::distanceSurface::print(Ostream& os) const +{ + os << "distanceSurface: " << name() << " :" + << " surface:" << surfPtr_().name() + << " distance:" << distance_ + << " faces:" << faces().size() + << " points:" << points().size(); +} + + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H new file mode 100644 index 0000000000000000000000000000000000000000..82fba8a83f3ff696122ffb23102d4838cb6f6c4d --- /dev/null +++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H @@ -0,0 +1,237 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::distanceSurface + +Description + A sampledSurface defined by a distance to a surface. + +SourceFiles + distanceSurface.C + +\*---------------------------------------------------------------------------*/ + +#ifndef distanceSurface_H +#define distanceSurface_H + +#include "sampledSurface.H" +#include "triSurface.H" +#include "searchableSurface.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class distanceSurface Declaration +\*---------------------------------------------------------------------------*/ + +class distanceSurface +: + public sampledSurface, + public triSurface +{ + // Private data + + //- Surface + const autoPtr<searchableSurface> surfPtr_; + + //- distance value + const scalar distance_; + + //- signed distance + const bool signed_; + + //- Whether to coarsen + const Switch regularise_; + + //- zone name (if restricted to zones) + word zoneName_; + + //- triangles converted to faceList + mutable autoPtr<faceList> facesPtr_; + + + // Recreated for every isoSurface + + //- Time at last call + mutable label storedTimeIndex_; + + //- For every triangle the original cell in mesh + mutable labelList meshCells_; + + + // Private Member Functions + + //- Create iso surface (if time has changed) + void createGeometry() const; + + //- sample field on faces + template <class Type> + tmp<Field<Type> > sampleField + ( + const GeometricField<Type, fvPatchField, volMesh>& vField + ) const; + + + template <class Type> + tmp<Field<Type> > + interpolateField(const interpolation<Type>&) const; + + +public: + + //- Runtime type information + TypeName("distanceSurface"); + + + // Constructors + + //- Construct from dictionary + distanceSurface + ( + const word& name, + const polyMesh& mesh, + const dictionary& dict + ); + + + // Destructor + + virtual ~distanceSurface(); + + + // Member Functions + + //- Points of surface + virtual const pointField& points() const + { + return triSurface::points(); + } + + //- Faces of surface + virtual const faceList& faces() const + { + if (!facesPtr_.valid()) + { + const triSurface& s = *this; + + facesPtr_.reset(new faceList(s.size())); + + forAll(s, i) + { + facesPtr_()[i] = s[i].triFaceFace(); + } + } + return facesPtr_; + } + + + //- Correct for mesh movement and/or field changes + virtual void correct(const bool meshChanged); + + + //- sample field on surface + virtual tmp<scalarField> sample + ( + const volScalarField& + ) const; + + //- sample field on surface + virtual tmp<vectorField> sample + ( + const volVectorField& + ) const; + + //- sample field on surface + virtual tmp<sphericalTensorField> sample + ( + const volSphericalTensorField& + ) const; + + //- sample field on surface + virtual tmp<symmTensorField> sample + ( + const volSymmTensorField& + ) const; + + //- sample field on surface + virtual tmp<tensorField> sample + ( + const volTensorField& + ) const; + + + //- interpolate field on surface + virtual tmp<scalarField> interpolate + ( + const interpolation<scalar>& + ) const; + + //- interpolate field on surface + virtual tmp<vectorField> interpolate + ( + const interpolation<vector>& + ) const; + + //- interpolate field on surface + virtual tmp<sphericalTensorField> interpolate + ( + const interpolation<sphericalTensor>& + ) const; + + //- interpolate field on surface + virtual tmp<symmTensorField> interpolate + ( + const interpolation<symmTensor>& + ) const; + + //- interpolate field on surface + virtual tmp<tensorField> interpolate + ( + const interpolation<tensor>& + ) const; + + //- Write + virtual void print(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "distanceSurfaceTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..ce983c7f5f08bffe29eda422a0f47dbff83d08cb --- /dev/null +++ b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C @@ -0,0 +1,83 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "distanceSurface.H" +#include "isoSurface.H" +#include "volFieldsFwd.H" +#include "pointFields.H" +#include "volPointInterpolation.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template <class Type> +Foam::tmp<Foam::Field<Type> > +Foam::distanceSurface::sampleField +( + const GeometricField<Type, fvPatchField, volMesh>& vField +) const +{ + return tmp<Field<Type> >(new Field<Type>(vField, meshCells_)); +} + + +template <class Type> +Foam::tmp<Foam::Field<Type> > +Foam::distanceSurface::interpolateField +( + const interpolation<Type>& interpolator +) const +{ + // One value per point + tmp<Field<Type> > tvalues(new Field<Type>(points().size())); + Field<Type>& values = tvalues(); + + boolList pointDone(points().size(), false); + + forAll(faces(), cutFaceI) + { + const face& f = faces()[cutFaceI]; + + forAll(f, faceVertI) + { + label pointI = f[faceVertI]; + + if (!pointDone[pointI]) + { + values[pointI] = interpolator.interpolate + ( + points()[pointI], + meshCells_[cutFaceI] + ); + pointDone[pointI] = true; + } + } + } + + return tvalues; +} + + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C index 45be2c1334b3f6c3540bdff3a831eccdd0547c1b..13c210d395bec3eafa6c287728dd9239954abba8 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurface.C +++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C @@ -25,12 +25,9 @@ License \*---------------------------------------------------------------------------*/ #include "isoSurface.H" -#include "dictionary.H" -#include "polyMesh.H" -#include "volFields.H" +#include "fvMesh.H" #include "mergePoints.H" -#include "tetMatcher.H" - +#include "syncTools.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -42,250 +39,1517 @@ namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +Foam::scalar Foam::isoSurface::isoFraction +( + const scalar s0, + const scalar s1 +) const +{ + scalar d = s1-s0; -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + if (mag(d) > VSMALL) + { + return (iso_-s0)/d; + } + else + { + return -1.0; + } +} -Foam::isoSurface::isoSurface + +// Determine for every face/cell whether it (possibly) generates triangles. +void Foam::isoSurface::calcCutTypes ( - const polyMesh& mesh, - const scalarField& cellValues, - const scalarField& pointValues, - const scalar iso + const volScalarField& cVals, + const scalarField& pVals ) -: - mesh_(mesh), - cellValues_(cellValues), - pointValues_(pointValues), - iso_(iso) { - const pointField& cellCentres = mesh.cellCentres(); + const polyBoundaryMesh& patches = mesh_.boundaryMesh(); + const labelList& own = mesh_.faceOwner(); + const labelList& nei = mesh_.faceNeighbour(); - tetMatcher tet; + faceCutType_.setSize(mesh_.nFaces()); + faceCutType_ = NOTCUT; - DynamicList<point> triPoints; - DynamicList<label> triMeshCells; - - forAll(mesh.cells(), cellI) + for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) { - label oldNPoints = triPoints.size(); - - const cell& cFaces = mesh.cells()[cellI]; + // CC edge. + bool ownLower = (cVals[own[faceI]] < iso_); + bool neiLower = (cVals[nei[faceI]] < iso_); - if (tet.isA(mesh, cellI)) + if (ownLower != neiLower) + { + faceCutType_[faceI] = CUT; + } + else { - // For tets don't do cell-centre decomposition, just use the - // tet points and values + // Mesh edge. + const face f = mesh_.faces()[faceI]; + + forAll(f, fp) + { + bool fpLower = (pVals[f[fp]] < iso_); + if + ( + (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)) + || (fpLower != ownLower) + || (fpLower != neiLower) + ) + { + faceCutType_[faceI] = CUT; + break; + } + } + } + } + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; - const face& f0 = mesh.faces()[cFaces[0]]; + label faceI = pp.start(); + forAll(pp, i) + { + bool ownLower = (cVals[own[faceI]] < iso_); + bool neiLower = (cVals.boundaryField()[patchI][i] < iso_); - // Get the other point - const face& f1 = mesh.faces()[cFaces[1]]; - label oppositeI = -1; - forAll(f1, fp) + if (ownLower != neiLower) + { + faceCutType_[faceI] = CUT; + } + else { - oppositeI = f1[fp]; + // Mesh edge. + const face f = mesh_.faces()[faceI]; - if (findIndex(f0, oppositeI) == -1) + forAll(f, fp) { - break; + bool fpLower = (pVals[f[fp]] < iso_); + if + ( + (fpLower != (pVals[f[f.fcIndex(fp)]] < iso_)) + || (fpLower != ownLower) + || (fpLower != neiLower) + ) + { + faceCutType_[faceI] = CUT; + break; + } } } + faceI++; + } + } + + + + nCutCells_ = 0; + cellCutType_.setSize(mesh_.nCells()); + cellCutType_ = NOTCUT; + + for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) + { + if (faceCutType_[faceI] != NOTCUT) + { + if (cellCutType_[own[faceI]] == NOTCUT) + { + cellCutType_[own[faceI]] = CUT; + nCutCells_++; + } + if (cellCutType_[nei[faceI]] == NOTCUT) + { + cellCutType_[nei[faceI]] = CUT; + nCutCells_++; + } + } + } + for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++) + { + if (faceCutType_[faceI] != NOTCUT) + { + if (cellCutType_[own[faceI]] == NOTCUT) + { + cellCutType_[own[faceI]] = CUT; + nCutCells_++; + } + } + } + + if (debug) + { + Pout<< "isoSurface : detected " << nCutCells_ + << " candidate cut cells." << endl; + } +} + + +// 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]); + } - vertexInterp + 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) +{ + vector sum = vector::zero; + + forAll(s, i) + { + sum += s[i].centre(s.points()); + } + return sum/s.size(); +} + + +// 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() == 0) + {} + else 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 ( - iso, - pointValues[f0[0]], - pointValues[f0[1]], - pointValues[f0[2]], - pointValues[oppositeI], - - mesh.points()[f0[0]], - mesh.points()[f0[1]], - mesh.points()[f0[2]], - mesh.points()[oppositeI], - - triPoints + 0.5 + * ( + tri0.centre(localPoints) + + tri1.centre(localPoints) + ) ); + info.setHit(); } - else + } + else + { + // 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; +} + + +// Get neighbour value and position. +void Foam::isoSurface::getNeighbour +( + const labelList& boundaryRegion, + const volScalarField& cVals, + const label cellI, + const label faceI, + scalar& nbrValue, + point& nbrPoint +) const +{ + const labelList& own = mesh_.faceOwner(); + const labelList& nei = mesh_.faceNeighbour(); + + if (mesh_.isInternalFace(faceI)) + { + label nbr = (own[faceI] == cellI ? nei[faceI] : own[faceI]); + nbrValue = cVals[nbr]; + nbrPoint = mesh_.C()[nbr]; + } + else + { + label bFaceI = faceI-mesh_.nInternalFaces(); + label patchI = boundaryRegion[bFaceI]; + label patchFaceI = faceI-mesh_.boundaryMesh()[patchI].start(); + + nbrValue = cVals.boundaryField()[patchI][patchFaceI]; + nbrPoint = mesh_.C().boundaryField()[patchI][patchFaceI]; + } +} + + +// Determine per cell centre whether all the intersections get collapsed +// to a single point +void Foam::isoSurface::calcSnappedCc +( + const labelList& boundaryRegion, + const volScalarField& cVals, + const scalarField& pVals, + + DynamicList<point>& snappedPoints, + labelList& snappedCc +) const +{ + const pointField& pts = mesh_.points(); + + snappedCc.setSize(mesh_.nCells()); + snappedCc = -1; + + // Work arrays + DynamicList<point, 64> localTriPoints(64); + + forAll(mesh_.cells(), cellI) + { + if (cellCutType_[cellI] == CUT) { - const cell& cFaces = mesh.cells()[cellI]; + scalar cVal = cVals[cellI]; + + const cell& cFaces = mesh_.cells()[cellI]; + + localTriPoints.clear(); + label nOther = 0; + point otherPointSum = vector::zero; + + // Create points for all intersections close to cell centre + // (i.e. from pyramid edges) forAll(cFaces, cFaceI) { label faceI = cFaces[cFaceI]; - const face& f = mesh.faces()[faceI]; - // Do a tetrahedrisation. Each face to cc becomes pyr. - // Each pyr gets split into two tets by diagionalisation - // of face. So - // - f[0], f[1], f[2], cc - // - f[0], f[2], f[3], cc + scalar nbrValue; + point nbrPoint; + getNeighbour + ( + boundaryRegion, + cVals, + cellI, + faceI, + nbrValue, + nbrPoint + ); + + // Calculate intersection points of edges to cell centre + FixedList<scalar, 3> s; + FixedList<point, 3> pt; - for(label fp = 1; fp < f.size() - 1; fp++) + // From cc to neighbour cc. + s[2] = isoFraction(cVal, nbrValue); + pt[2] = (1.0-s[2])*mesh_.C()[cellI] + s[2]*nbrPoint; + + const face& f = mesh_.faces()[cFaces[cFaceI]]; + + forAll(f, fp) + { + // From cc to fp + label p0 = f[fp]; + s[0] = isoFraction(cVal, pVals[p0]); + pt[0] = (1.0-s[0])*mesh_.C()[cellI] + s[0]*pts[p0]; + + // From cc to fp+1 + label p1 = f[f.fcIndex(fp)]; + s[1] = isoFraction(cVal, pVals[p1]); + pt[1] = (1.0-s[1])*mesh_.C()[cellI] + s[1]*pts[p1]; + + if + ( + (s[0] >= 0.0 && s[0] <= 0.5) + && (s[1] >= 0.0 && s[1] <= 0.5) + && (s[2] >= 0.0 && s[2] <= 0.5) + ) + { + localTriPoints.append(pt[0]); + localTriPoints.append(pt[1]); + localTriPoints.append(pt[2]); + } + else + { + // Get average of all other points + forAll(s, i) + { + if (s[i] >= 0.0 && s[i] <= 0.5) + { + otherPointSum += pt[i]; + nOther++; + } + } + } + } + } + + if (localTriPoints.size() == 0) + { + // No complete triangles. Get average of all intersection + // points. + if (nOther > 0) { - vertexInterp + snappedCc[cellI] = snappedPoints.size(); + snappedPoints.append(otherPointSum/nOther); + + //Pout<< " point:" << pointI + // << " replacing coord:" << mesh_.points()[pointI] + // << " by average:" << collapsedPoint[pointI] << endl; + } + } + else if (localTriPoints.size() == 3) + { + // Single triangle. No need for any analysis. Average points. + pointField points; + points.transfer(localTriPoints); + snappedCc[cellI] = snappedPoints.size(); + snappedPoints.append(sum(points)/points.size()); + + //Pout<< " point:" << pointI + // << " replacing coord:" << mesh_.points()[pointI] + // << " by average:" << collapsedPoint[pointI] << endl; + } + else + { + // Convert points into triSurface. + + // Merge points and compact out non-valid triangles + labelList triMap; // merged to unmerged triangle + labelList triPointReverseMap; // unmerged to merged point + triSurface surf + ( + stitchTriPoints ( - iso, - pointValues[f[0]], - pointValues[f[fp]], - pointValues[f[f.fcIndex(fp)]], - cellValues[cellI], - - mesh.points()[f[0]], - mesh.points()[f[fp]], - mesh.points()[f[f.fcIndex(fp)]], - cellCentres[cellI], - - triPoints - ); + false, // do not check for duplicate tris + localTriPoints, + triPointReverseMap, + triMap + ) + ); + + labelList faceZone; + label nZones = surf.markZones + ( + boolList(surf.nEdges(), false), + faceZone + ); + + if (nZones == 1) + { + snappedCc[cellI] = snappedPoints.size(); + snappedPoints.append(calcCentre(surf)); + //Pout<< " point:" << pointI << " nZones:" << nZones + // << " replacing coord:" << mesh_.points()[pointI] + // << " by average:" << collapsedPoint[pointI] << endl; + } + } + } + } +} + + +// Determine per meshpoint whether all the intersections get collapsed +// to a single point +void Foam::isoSurface::calcSnappedPoint +( + const PackedList<1>& isBoundaryPoint, + const labelList& boundaryRegion, + const volScalarField& cVals, + const scalarField& pVals, + + DynamicList<point>& snappedPoints, + labelList& snappedPoint +) const +{ + const pointField& pts = mesh_.points(); + + + const point greatPoint(VGREAT, VGREAT, VGREAT); + pointField collapsedPoint(mesh_.nPoints(), greatPoint); + + + // Work arrays + DynamicList<point, 64> localTriPoints(100); + + forAll(mesh_.pointFaces(), pointI) + { + if (isBoundaryPoint.get(pointI) == 1) + { + continue; + } + + const labelList& pFaces = mesh_.pointFaces()[pointI]; + + bool anyCut = false; + + forAll(pFaces, i) + { + label faceI = pFaces[i]; + + if (faceCutType_[faceI] == CUT) + { + anyCut = true; + break; + } + } + + if (!anyCut) + { + continue; + } + + + localTriPoints.clear(); + label nOther = 0; + point otherPointSum = vector::zero; + + forAll(pFaces, pFaceI) + { + label faceI = pFaces[pFaceI]; + const face& f = mesh_.faces()[faceI]; + label own = mesh_.faceOwner()[faceI]; + + // Create points for all intersections close to point + // (i.e. from pyramid edges) + + scalar nbrValue; + point nbrPoint; + getNeighbour + ( + boundaryRegion, + cVals, + own, + faceI, + nbrValue, + nbrPoint + ); + + // Calculate intersection points of edges emanating from point + FixedList<scalar, 4> s; + FixedList<point, 4> pt; + + label fp = findIndex(f, pointI); + s[0] = isoFraction(pVals[pointI], cVals[own]); + pt[0] = (1.0-s[0])*pts[pointI] + s[0]*mesh_.C()[own]; + + s[1] = isoFraction(pVals[pointI], nbrValue); + pt[1] = (1.0-s[1])*pts[pointI] + s[1]*nbrPoint; + + label nextPointI = f[f.fcIndex(fp)]; + s[2] = isoFraction(pVals[pointI], pVals[nextPointI]); + pt[2] = (1.0-s[2])*pts[pointI] + s[2]*pts[nextPointI]; + + label prevPointI = f[f.rcIndex(fp)]; + s[3] = isoFraction(pVals[pointI], pVals[prevPointI]); + pt[3] = (1.0-s[3])*pts[pointI] + s[3]*pts[prevPointI]; + + if + ( + (s[0] >= 0.0 && s[0] <= 0.5) + && (s[1] >= 0.0 && s[1] <= 0.5) + && (s[2] >= 0.0 && s[2] <= 0.5) + ) + { + localTriPoints.append(pt[0]); + localTriPoints.append(pt[1]); + localTriPoints.append(pt[2]); + } + if + ( + (s[0] >= 0.0 && s[0] <= 0.5) + && (s[1] >= 0.0 && s[1] <= 0.5) + && (s[3] >= 0.0 && s[3] <= 0.5) + ) + { + localTriPoints.append(pt[3]); + localTriPoints.append(pt[0]); + localTriPoints.append(pt[1]); + } + + // Get average of points as fallback + forAll(s, i) + { + if (s[i] >= 0.0 && s[i] <= 0.5) + { + otherPointSum += pt[i]; + nOther++; } } } + if (localTriPoints.size() == 0) + { + // No complete triangles. Get average of all intersection + // points. + if (nOther > 0) + { + collapsedPoint[pointI] = otherPointSum/nOther; + } + } + else if (localTriPoints.size() == 3) + { + // Single triangle. No need for any analysis. Average points. + pointField points; + points.transfer(localTriPoints); + collapsedPoint[pointI] = sum(points)/points.size(); + } + else + { + // Convert points into triSurface. + + // Merge points and compact out non-valid triangles + labelList triMap; // merged to unmerged triangle + labelList triPointReverseMap; // unmerged to merged point + triSurface surf + ( + stitchTriPoints + ( + false, // do not check for duplicate tris + localTriPoints, + triPointReverseMap, + triMap + ) + ); + + labelList faceZone; + label nZones = surf.markZones + ( + boolList(surf.nEdges(), false), + faceZone + ); + + if (nZones == 1) + { + collapsedPoint[pointI] = calcCentre(surf); + } + } + } + + syncTools::syncPointList + ( + mesh_, + collapsedPoint, + minEqOp<point>(), + greatPoint, + true // are coordinates so separate + ); + + snappedPoint.setSize(mesh_.nPoints()); + snappedPoint = -1; - // Every three triPoints is a cell - label nCells = (triPoints.size()-oldNPoints)/3; - for (label i = 0; i < nCells; i++) + forAll(collapsedPoint, pointI) + { + if (collapsedPoint[pointI] != greatPoint) { - triMeshCells.append(cellI); + snappedPoint[pointI] = snappedPoints.size(); + snappedPoints.append(collapsedPoint[pointI]); } } +} + + +Foam::triSurface Foam::isoSurface::stitchTriPoints +( + const bool checkDuplicates, + const List<point>& triPoints, + labelList& triPointReverseMap, // unmerged to merged point + labelList& triMap // merged to unmerged triangle +) const +{ + label nTris = triPoints.size()/3; - triPoints.shrink(); - triMeshCells.shrink(); - meshCells_.transfer(triMeshCells); + if ((triPoints.size() % 3) != 0) + { + FatalErrorIn("isoSurface::stitchTriPoints(..)") + << "Problem: number of points " << triPoints.size() + << " not a multiple of 3." << abort(FatalError); + } pointField newPoints; - mergePoints(triPoints, SMALL, false, triPointMergeMap_, newPoints); + mergePoints + ( + triPoints, + mergeDistance_, + false, + triPointReverseMap, + newPoints + ); - DynamicList<labelledTri> tris(meshCells_.size()); - forAll(meshCells_, triI) + // Check that enough merged. + if (debug) { - tris.append + pointField newNewPoints; + labelList oldToNew; + bool hasMerged = mergePoints ( - labelledTri - ( - triPointMergeMap_[3*triI], - triPointMergeMap_[3*triI+1], - triPointMergeMap_[3*triI+2], - 0 - ) + newPoints, + mergeDistance_, + true, + oldToNew, + newNewPoints ); + + if (hasMerged) + { + FatalErrorIn("isoSurface::stitchTriPoints(..)") + << "Merged points contain duplicates" + << " when merging with distance " << mergeDistance_ << endl + << "merged:" << newPoints.size() << " re-merged:" + << newNewPoints.size() + << abort(FatalError); + } } - //- 1.5.x: - //tris.shrink(); - //triSurface::operator= - //( - // triSurface(tris, geometricSurfacePatchList(0), newPoints) - //); - triSurface::operator= - ( - triSurface(tris, geometricSurfacePatchList(0), newPoints, true) - ); -} + List<labelledTri> tris; + { + DynamicList<labelledTri> dynTris(nTris); + label rawPointI = 0; + DynamicList<label> newToOldTri(nTris); + + for (label oldTriI = 0; oldTriI < nTris; oldTriI++) + { + labelledTri tri + ( + triPointReverseMap[rawPointI], + triPointReverseMap[rawPointI+1], + triPointReverseMap[rawPointI+2], + 0 + ); + rawPointI += 3; + if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2])) + { + newToOldTri.append(oldTriI); + dynTris.append(tri); + } + } -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -//Foam::tmp<Foam::scalarField> -//Foam::isoSurface::sample -//( -// const volScalarField& vField -//) const -//{ -// return sampleField(vField); -//} -// -// -//Foam::tmp<Foam::vectorField> -//Foam::isoSurface::sample -//( -// const volVectorField& vField -//) const -//{ -// return sampleField(vField); -//} -// -//Foam::tmp<Foam::sphericalTensorField> -//Foam::isoSurface::sample -//( -// const volSphericalTensorField& vField -//) const -//{ -// return sampleField(vField); -//} -// -// -//Foam::tmp<Foam::symmTensorField> -//Foam::isoSurface::sample -//( -// const volSymmTensorField& vField -//) const -//{ -// return sampleField(vField); -//} -// -// -//Foam::tmp<Foam::tensorField> -//Foam::isoSurface::sample -//( -// const volTensorField& vField -//) const -//{ -// return sampleField(vField); -//} -// -// -//Foam::tmp<Foam::scalarField> -//Foam::isoSurface::interpolate -//( -// const interpolation<scalar>& interpolator -//) const -//{ -// return interpolateField(interpolator); -//} -// -// -//Foam::tmp<Foam::vectorField> -//Foam::isoSurface::interpolate -//( -// const interpolation<vector>& interpolator -//) const -//{ -// return interpolateField(interpolator); -//} -// -//Foam::tmp<Foam::sphericalTensorField> -//Foam::isoSurface::interpolate -//( -// const interpolation<sphericalTensor>& interpolator -//) const -//{ -// return interpolateField(interpolator); -//} -// -// -//Foam::tmp<Foam::symmTensorField> -//Foam::isoSurface::interpolate -//( -// const interpolation<symmTensor>& interpolator -//) const -//{ -// return interpolateField(interpolator); -//} -// -// -//Foam::tmp<Foam::tensorField> -//Foam::isoSurface::interpolate -//( -// const interpolation<tensor>& interpolator -//) const -//{ -// return interpolateField(interpolator); -//} + triMap.transfer(newToOldTri.shrink()); + tris.transfer(dynTris.shrink()); + } + + + + // Use face centres to determine 'flat hole' situation (see RMT paper). + // Two unconnected triangles get connected because (some of) the edges + // separating them get collapsed. Below only checks for duplicate triangles, + // not non-manifold edge connectivity. + if (checkDuplicates) + { + if (debug) + { + Pout<< "isoSurface : merged from " << nTris + << " down to " << tris.size() << " triangles." << endl; + } + + pointField centres(tris.size()); + forAll(tris, triI) + { + centres[triI] = tris[triI].centre(newPoints); + } + + pointField mergedCentres; + labelList oldToMerged; + bool hasMerged = mergePoints + ( + centres, + mergeDistance_, + false, + oldToMerged, + mergedCentres + ); + + if (debug) + { + Pout<< "isoSurface : detected " + << centres.size()-mergedCentres.size() + << " duplicate triangles." << endl; + } + + if (hasMerged) + { + // Filter out duplicates. + label newTriI = 0; + DynamicList<label> newToOldTri(tris.size()); + labelList newToMaster(mergedCentres.size(), -1); + forAll(tris, triI) + { + label mergedI = oldToMerged[triI]; + + if (newToMaster[mergedI] == -1) + { + newToMaster[mergedI] = triI; + newToOldTri.append(triMap[triI]); + tris[newTriI++] = tris[triI]; + } + } + + triMap.transfer(newToOldTri.shrink()); + tris.setSize(newTriI); + } + } + + return triSurface(tris, geometricSurfacePatchList(0), newPoints, true); +} + + +// Does face use valid vertices? +bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI) +{ + // Simple check on indices ok. + + const labelledTri& f = surf[faceI]; + + if + ( + (f[0] < 0) || (f[0] >= surf.points().size()) + || (f[1] < 0) || (f[1] >= surf.points().size()) + || (f[2] < 0) || (f[2] >= surf.points().size()) + ) + { + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI << " vertices " << f + << " uses point indices outside point range 0.." + << surf.points().size()-1 << endl; + + return false; + } + + if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2])) + { + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI + << " uses non-unique vertices " << f + << endl; + return false; + } + + // duplicate triangle check + + const labelList& fFaces = surf.faceFaces()[faceI]; + + // Check if faceNeighbours use same points as this face. + // Note: discards normal information - sides of baffle are merged. + forAll(fFaces, i) + { + label nbrFaceI = fFaces[i]; + + if (nbrFaceI <= faceI) + { + // lower numbered faces already checked + continue; + } + + const labelledTri& nbrF = surf[nbrFaceI]; + + if + ( + ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2])) + && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2])) + && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2])) + ) + { + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI << " vertices " << f + << " fc:" << f.centre(surf.points()) + << " has the same vertices as triangle " << nbrFaceI + << " vertices " << nbrF + << " fc:" << nbrF.centre(surf.points()) + << endl; + + return false; + } + } + return true; +} + + +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() + << " 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(); + + forAll(oldToMerged, oldEdgeI) + { + label triI = oldEdgeI / 3; + label edgeI = oldToMerged[oldEdgeI]; + + if (edgeFace0[edgeI] == -1) + { + edgeFace0[edgeI] = triI; + } + else if (edgeFace1[edgeI] == -1) + { + edgeFace1[edgeI] = triI; + } + else + { + //WarningIn("orientSurface(triSurface&)") + // << "Edge " << edgeI << " with centre " << mergedCentres[edgeI] + // << " used by more than two triangles: " << edgeFace0[edgeI] + // << ", " + // << edgeFace1[edgeI] << " and " << triI << endl; + Map<labelList>::iterator iter = edgeFacesRest.find(edgeI); + + if (iter != edgeFacesRest.end()) + { + labelList& eFaces = iter(); + label sz = eFaces.size(); + eFaces.setSize(sz+1); + eFaces[sz] = triI; + } + else + { + edgeFacesRest.insert(edgeI, labelList(1, triI)); + } + } + } +} + + +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() > 0) + { + 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::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) + { + FatalErrorIn + ( + "isoSurface::orientSurface(triSurface&, const label)" + ) << "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, + const labelList& newToOldFaces, + labelList& oldToNewPoints, + labelList& newToOldPoints +) +{ + const boolList include + ( + createWithValues<boolList> + ( + s.size(), + false, + newToOldFaces, + true + ) + ); + + newToOldPoints.setSize(s.points().size()); + oldToNewPoints.setSize(s.points().size()); + oldToNewPoints = -1; + { + label pointI = 0; + + forAll(include, oldFacei) + { + if (include[oldFacei]) + { + // Renumber labels for triangle + const labelledTri& tri = s[oldFacei]; + + forAll(tri, fp) + { + label oldPointI = tri[fp]; + + if (oldToNewPoints[oldPointI] == -1) + { + oldToNewPoints[oldPointI] = pointI; + newToOldPoints[pointI++] = oldPointI; + } + } + } + } + newToOldPoints.setSize(pointI); + } + + // Extract points + pointField newPoints(newToOldPoints.size()); + forAll(newToOldPoints, i) + { + newPoints[i] = s.points()[newToOldPoints[i]]; + } + // Extract faces + List<labelledTri> newTriangles(newToOldFaces.size()); + + forAll(newToOldFaces, i) + { + // Get old vertex labels + const labelledTri& tri = s[newToOldFaces[i]]; + + newTriangles[i][0] = oldToNewPoints[tri[0]]; + newTriangles[i][1] = oldToNewPoints[tri[1]]; + newTriangles[i][2] = oldToNewPoints[tri[2]]; + newTriangles[i].region() = tri.region(); + } + + // Reuse storage. + return triSurface(newTriangles, s.patches(), newPoints, true); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::isoSurface::isoSurface +( + const volScalarField& cVals, + const scalarField& pVals, + const scalar iso, + const bool regularise, + const scalar mergeTol +) +: + mesh_(cVals.mesh()), + iso_(iso), + mergeDistance_(mergeTol*mag(mesh_.bounds().max()-mesh_.bounds().min())) +{ + // Determine if any cut through face/cell + calcCutTypes(cVals, pVals); + + + // Determine if point is on boundary. Points on boundaries are never + // snapped. Coupled boundaries are handled explicitly so not marked here. + PackedList<1> isBoundaryPoint(mesh_.nPoints()); + + labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces()); + const polyBoundaryMesh& patches = mesh_.boundaryMesh(); + const labelList& own = mesh_.faceOwner(); + + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + + if (!pp.coupled()) + { + label faceI = pp.start(); + + forAll(pp, i) + { + boundaryRegion[faceI-mesh_.nInternalFaces()] = patchI; + + const face& f = mesh_.faces()[faceI++]; + + forAll(f, fp) + { + isBoundaryPoint.set(f[fp], 1); + } + } + } + } + + + DynamicList<point> snappedPoints(nCutCells_); + + // Per cc -1 or a point inside snappedPoints. + labelList snappedCc; + if (regularise) + { + calcSnappedCc + ( + boundaryRegion, + cVals, + pVals, + + snappedPoints, + snappedCc + ); + } + else + { + snappedCc.setSize(mesh_.nCells()); + snappedCc = -1; + } + + // Determine neighbouring snap status + labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1); + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + + if (pp.coupled()) + { + label faceI = pp.start(); + forAll(pp, i) + { + neiSnappedCc[faceI-mesh_.nInternalFaces()] = + snappedCc[own[faceI]]; + faceI++; + } + } + } + syncTools::swapBoundaryFaceList(mesh_, neiSnappedCc, false); + + + if (debug) + { + Pout<< "isoSurface : shifted " << snappedPoints.size() + << " cell centres to intersection." << endl; + } + + label nCellSnaps = snappedPoints.size(); + + // Per point -1 or a point inside snappedPoints. + labelList snappedPoint; + if (regularise) + { + calcSnappedPoint + ( + isBoundaryPoint, + boundaryRegion, + cVals, + pVals, + + snappedPoints, + snappedPoint + ); + } + else + { + snappedPoint.setSize(mesh_.nPoints()); + snappedPoint = -1; + } + + if (debug) + { + Pout<< "isoSurface : shifted " << snappedPoints.size()-nCellSnaps + << " vertices to intersection." << endl; + } + + + DynamicList<point> triPoints(nCutCells_); + DynamicList<label> triMeshCells(nCutCells_); + + generateTriPoints + ( + cVals, + pVals, + + mesh_.C(), + mesh_.points(), + + snappedPoints, + snappedCc, + snappedPoint, + + triPoints, + triMeshCells + ); + + if (debug) + { + Pout<< "isoSurface : generated " << triMeshCells.size() + << " unmerged triangles." << endl; + } + + + // Merge points and compact out non-valid triangles + labelList triMap; // merged to unmerged triangle + triSurface::operator= + ( + stitchTriPoints + ( + true, // check for duplicate tris + triPoints, + triPointMergeMap_, // unmerged to merged point + triMap + ) + ); + + if (debug) + { + Pout<< "isoSurface : generated " << triMap.size() + << " merged triangles." << endl; + } + + meshCells_.setSize(triMap.size()); + forAll(triMap, i) + { + meshCells_[i] = triMeshCells[triMap[i]]; + } + + if (debug) + { + Pout<< "isoSurface : checking " << size() + << " triangles for validity." << endl; + + forAll(*this, triI) + { + // Copied from surfaceCheck + validTri(*this, triI); + } + } + + + //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); + //} +} // ************************************************************************* // diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H index ddd6854909f6529c5986c598cf30b1fbf7c6d657..28fc8b7c4a6ab68f579b4d2d72dc1509889875d9 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurface.H +++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H @@ -27,7 +27,20 @@ Class Description A surface formed by the iso value. - After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke + After "Regularised Marching Tetrahedra: improved iso-surface extraction", + G.M. Treece, R.W. Prager and A.H. Gee. + + Note: + - 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! + - does tets without using cell centres/cell values. Not tested. + - regularisation can create duplicate triangles/non-manifold surfaces. + Current handling of those is bit ad-hoc for now and not perfect. + - uses geometric merge with fraction of bounding box as distance. + - triangles can be between two cell centres so constant sampling + does not make sense. + - does not do 2D correctly, creates non-flat iso surface. SourceFiles isoSurface.C @@ -38,13 +51,17 @@ SourceFiles #define isoSurface_H #include "triSurface.H" +#include "labelPair.H" +#include "pointIndexHit.H" +#include "PackedList.H" +#include "volFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { -class polyMesh; +class fvMesh; /*---------------------------------------------------------------------------*\ Class isoSurface Declaration @@ -56,18 +73,39 @@ class isoSurface { // Private data - //- Reference to mesh - const polyMesh& mesh_; + enum segmentCutType + { + NEARFIRST, // intersection close to e.first() + NEARSECOND, // ,, e.second() + NOTNEAR // no intersection + }; - //- Reference to cell values - const scalarField& cellValues_; + enum cellCutType + { + NOTCUT, // not cut + SPHERE, // all edges to cell centre cut + CUT // normal cut + }; - //- Reference to point values - const scalarField& pointValues_; + + //- Reference to mesh + const fvMesh& mesh_; //- Isosurface value const scalar iso_; + //- When to merge points + const scalar mergeDistance_; + + + //- Whether face might be cut + List<cellCutType> faceCutType_; + + //- Whether cell might be cut + List<cellCutType> cellCutType_; + + //- Estimated number of cut cells + label nCutCells_; //- For every triangle the original cell in mesh labelList meshCells_; @@ -78,31 +116,210 @@ class isoSurface // Private Member Functions - template<class T> - static T vertexInterp + //- Get location of iso value as fraction inbetween s0,s1 + scalar isoFraction ( - const scalar iso, - const T& p0, - const T& p1, const scalar s0, const scalar s1 + ) const; + + //- Set faceCutType,cellCutType. + void calcCutTypes + ( + const volScalarField& cVals, + const scalarField& pVals ); - template<class T> - static void vertexInterp + static labelPair findCommonPoints ( - const scalar iso, + const labelledTri&, + const labelledTri& + ); + + static point calcCentre(const triSurface&); + + static pointIndexHit collapseSurface + ( + pointField& localPoints, + DynamicList<labelledTri, 64>& localTris + ); + + void getNeighbour + ( + const labelList& boundaryRegion, + const volScalarField& cVals, + const label cellI, + const label faceI, + scalar& nbrValue, + point& nbrPoint + ) const; + + //- Determine per cc whether all near cuts can be snapped to single + // point. + void calcSnappedCc + ( + const labelList& boundaryRegion, + const volScalarField& cVals, + const scalarField& pVals, + DynamicList<point>& snappedPoints, + labelList& snappedCc + ) const; + + //- Determine per point whether all near cuts can be snapped to single + // point. + void calcSnappedPoint + ( + const PackedList<1>& isBoundaryPoint, + const labelList& boundaryRegion, + const volScalarField& cVals, + const scalarField& pVals, + DynamicList<point>& snappedPoints, + labelList& snappedPoint + ) const; + + //- Generate single point by interpolation or snapping + template<class Type> + Type generatePoint + ( + const DynamicList<Type>& snappedPoints, + + const scalar s0, + const Type& p0, + const label p0Index, + + const scalar s1, + const Type& p1, + const label p1Index + ) const; + + template<class Type> + void generateTriPoints + ( + const DynamicList<Type>& snapped, + const scalar s0, + const Type& p0, + const label p0Index, + const scalar s1, + const Type& p1, + const label p1Index, + const scalar s2, + const Type& p2, + const label p2Index, + const scalar s3, + const Type& p3, + const label p3Index, + + DynamicList<Type>& points + ) const; + + template<class Type> + label generateTriPoints + ( + const volScalarField& cVals, + const scalarField& pVals, + + const GeometricField<Type, fvPatchField, volMesh>& cCoords, + const Field<Type>& pCoords, + + const DynamicList<Type>& snappedPoints, + const labelList& snappedCc, + const labelList& snappedPoint, + const label faceI, + + const scalar neiVal, + const Type& neiPt, + const label neiSnap, + + DynamicList<Type>& triPoints, + DynamicList<label>& triMeshCells + ) const; + + template<class Type> + void generateTriPoints + ( + const volScalarField& cVals, + const scalarField& pVals, + + const GeometricField<Type, fvPatchField, volMesh>& cCoords, + const Field<Type>& pCoords, + + const DynamicList<Type>& snappedPoints, + const labelList& snappedCc, + const labelList& snappedPoint, + + DynamicList<Type>& triPoints, + DynamicList<label>& triMeshCells + ) const; + + triSurface stitchTriPoints + ( + const bool checkDuplicates, + const List<point>& triPoints, + labelList& triPointReverseMap, // unmerged to merged point + labelList& triMap // merged to unmerged triangle + ) const; + + //- Check single triangle for (topological) validity + static bool validTri(const triSurface&, const label); - const T& p0, - const T& p1, - const T& p2, - const T& p3, + //- 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 + ); - DynamicList<T>& points + //- 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, + const labelList& newToOldFaces, + labelList& oldToNewPoints, + labelList& newToOldPoints ); public: @@ -113,13 +330,16 @@ public: // Constructors - //- Construct from dictionary + //- Construct from cell values and point values. Uses boundaryField + // for boundary values. Requires on coupled patchfields to be set + // to the opposite cell value. isoSurface ( - const polyMesh& mesh, - const scalarField& cellValues, - const scalarField& pointValues, - const scalar iso + const volScalarField& cellIsoVals, + const scalarField& pointIsoVals, + const scalar iso, + const bool regularise, + const scalar mergeTol = 1E-6 // fraction of bounding box ); @@ -132,27 +352,22 @@ public: } //- For every unmerged triangle point the point in the triSurface - const labelList triPointMergeMap() const + const labelList& triPointMergeMap() const { return triPointMergeMap_; } - - //- sample field on faces + //- Interpolates cCoords,pCoords. Takes the original fields + // used to create the iso surface. template <class Type> - tmp<Field<Type> > sample + tmp<Field<Type> > interpolate ( - const Field<Type>& sampleCellValues + const volScalarField& cellIsoVals, + const scalarField& pointIsoVals, + const GeometricField<Type, fvPatchField, volMesh>& cCoords, + const Field<Type>& pCoords ) const; - //- interpolate field to points - template <class Type> - tmp<Field<Type> > - interpolate - ( - const Field<Type>& sampleCellValues, - const Field<Type>& samplePointValues - ) const; }; diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C new file mode 100644 index 0000000000000000000000000000000000000000..5cffcdf0ddf4b8f21c8f2f2f01363ad97ec43df9 --- /dev/null +++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C @@ -0,0 +1,1753 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "isoSurfaceCell.H" +#include "dictionary.H" +#include "polyMesh.H" +#include "mergePoints.H" +#include "tetMatcher.H" +#include "syncTools.H" +#include "addToRunTimeSelectionTable.H" +#include "faceTriangulation.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(isoSurfaceCell, 0); +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::scalar Foam::isoSurfaceCell::isoFraction +( + const scalar s0, + const scalar s1 +) const +{ + scalar d = s1-s0; + + if (mag(d) > VSMALL) + { + return (iso_-s0)/d; + } + else + { + return -1.0; + } +} + + +//Foam::List<Foam::triFace> Foam::isoSurfaceCell::triangulate(const face& f) +// const +//{ +// faceList triFaces(f.nTriangles(mesh_.points())); +// label triFaceI = 0; +// f.triangles(mesh_.points(), triFaceI, triFaces); +// +// List<triFace> tris(triFaces.size()); +// forAll(triFaces, i) +// { +// tris[i][0] = triFaces[i][0]; +// tris[i][1] = triFaces[i][1]; +// tris[i][2] = triFaces[i][2]; +// } +// return tris; +//} + + +Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType +( + const PackedList<1>& isTet, + const scalarField& cellValues, + const scalarField& pointValues, + const label cellI +) const +{ + const cell& cFaces = mesh_.cells()[cellI]; + + if (isTet.get(cellI) == 1) + { + forAll(cFaces, cFaceI) + { + const face& f = mesh_.faces()[cFaces[cFaceI]]; + + for (label fp = 1; fp < f.size() - 1; fp++) + { + triFace tri(f[0], f[fp], f[f.fcIndex(fp)]); + + bool aLower = (pointValues[tri[0]] < iso_); + bool bLower = (pointValues[tri[1]] < iso_); + bool cLower = (pointValues[tri[2]] < iso_); + + if (aLower == bLower && aLower == cLower) + {} + else + { + return CUT; + } + } + } + return NOTCUT; + } + else + { + bool cellLower = (cellValues[cellI] < iso_); + + // First check if there is any cut in cell + bool edgeCut = false; + + forAll(cFaces, cFaceI) + { + const face& f = mesh_.faces()[cFaces[cFaceI]]; + + // Check pyramids cut + forAll(f, fp) + { + if ((pointValues[f[fp]] < iso_) != cellLower) + { + edgeCut = true; + break; + } + } + + if (edgeCut) + { + break; + } + + for (label fp = 1; fp < f.size() - 1; fp++) + { + triFace tri(f[0], f[fp], f[f.fcIndex(fp)]); + //List<triFace> tris(triangulate(f)); + //forAll(tris, i) + //{ + // const triFace& tri = tris[i]; + + bool aLower = (pointValues[tri[0]] < iso_); + bool bLower = (pointValues[tri[1]] < iso_); + bool cLower = (pointValues[tri[2]] < iso_); + + if (aLower == bLower && aLower == cLower) + {} + else + { + edgeCut = true; + break; + } + } + + if (edgeCut) + { + break; + } + } + + if (edgeCut) + { + // Count actual cuts (expensive since addressing needed) + // Note: not needed if you don't want to preserve maxima/minima + // centred around cellcentre. In that case just always return CUT + + const labelList& cPoints = mesh_.cellPoints(cellI); + + label nPyrCuts = 0; + + forAll(cPoints, i) + { + if ((pointValues[cPoints[i]] < iso_) != cellLower) + { + nPyrCuts++; + } + } + + if (nPyrCuts == cPoints.size()) + { + return SPHERE; + } + else + { + return CUT; + } + } + else + { + return NOTCUT; + } + } +} + + +void Foam::isoSurfaceCell::calcCutTypes +( + const PackedList<1>& isTet, + const scalarField& cVals, + const scalarField& pVals +) +{ + cellCutType_.setSize(mesh_.nCells()); + nCutCells_ = 0; + forAll(mesh_.cells(), cellI) + { + cellCutType_[cellI] = calcCutType(isTet, cVals, pVals, cellI); + + if (cellCutType_[cellI] == CUT) + { + nCutCells_++; + } + } + + if (debug) + { + Pout<< "isoSurfaceCell : detected " << nCutCells_ + << " candidate cut cells." << endl; + } +} + + + +// Return the two common points between two triangles +Foam::labelPair Foam::isoSurfaceCell::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::isoSurfaceCell::calcCentre(const triSurface& s) +{ + vector sum = vector::zero; + + forAll(s, i) + { + sum += s[i].centre(s.points()); + } + return sum/s.size(); +} + + +// Replace surface (localPoints, localTris) with single point. Returns +// point. Destructs arguments. +Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface +( + pointField& localPoints, + DynamicList<labelledTri, 64>& localTris +) +{ + pointIndexHit info(false, vector::zero, localTris.size()); + + if (localTris.size() == 0) + {} + else 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 + { + // 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; +} + + +void Foam::isoSurfaceCell::calcSnappedCc +( + const PackedList<1>& isTet, + const scalarField& cVals, + const scalarField& pVals, + + DynamicList<point>& snappedPoints, + labelList& snappedCc +) const +{ + const pointField& cc = mesh_.cellCentres(); + const pointField& pts = mesh_.points(); + + snappedCc.setSize(mesh_.nCells()); + snappedCc = -1; + + // Work arrays + DynamicList<point, 64> localPoints(64); + DynamicList<labelledTri, 64> localTris(64); + Map<label> pointToLocal(64); + + forAll(mesh_.cells(), cellI) + { + if (cellCutType_[cellI] == CUT && isTet.get(cellI) == 0) + { + scalar cVal = cVals[cellI]; + + const cell& cFaces = mesh_.cells()[cellI]; + + localPoints.clear(); + localTris.clear(); + pointToLocal.clear(); + + // Create points for all intersections close to cell centre + // (i.e. from pyramid edges) + + forAll(cFaces, cFaceI) + { + const face& f = mesh_.faces()[cFaces[cFaceI]]; + + forAll(f, fp) + { + label pointI = f[fp]; + + scalar s = isoFraction(cVal, pVals[pointI]); + + if (s >= 0.0 && s <= 0.5) + { + if (pointToLocal.insert(pointI, localPoints.size())) + { + localPoints.append((1.0-s)*cc[cellI]+s*pts[pointI]); + } + } + } + } + + if (localPoints.size() == 0) + { + // No near intersections + } + else if (localPoints.size() == 1) + { + // No need for any analysis. + snappedCc[cellI] = snappedPoints.size(); + snappedPoints.append(localPoints[0]); + + //Pout<< "cell:" << cellI + // << " at " << mesh_.cellCentres()[cellI] + // << " collapsing " << localPoints + // << " intersections down to " + // << snappedPoints[snappedCc[cellI]] << endl; + } + else if (localPoints.size() == 2) + { + //? No need for any analysis.??? + snappedCc[cellI] = snappedPoints.size(); + snappedPoints.append(0.5*(localPoints[0]+localPoints[1])); + + //Pout<< "cell:" << cellI + // << " at " << mesh_.cellCentres()[cellI] + // << " collapsing " << localPoints + // << " intersections down to " + // << snappedPoints[snappedCc[cellI]] << endl; + } + else + { + // Need to analyse + forAll(cFaces, cFaceI) + { + const face& f = mesh_.faces()[cFaces[cFaceI]]; + + // Do a tetrahedrisation. Each face to cc becomes pyr. + // Each pyr gets split into tets by diagonalisation + // of face. + + for (label fp = 1; fp < f.size() - 1; fp++) + { + triFace tri(f[0], f[fp], f[f.fcIndex(fp)]); + //List<triFace> tris(triangulate(f)); + //forAll(tris, i) + //{ + // const triFace& tri = tris[i]; + + // Get fractions for the three edges to cell centre + FixedList<scalar, 3> s(3); + s[0] = isoFraction(cVal, pVals[tri[0]]); + s[1] = isoFraction(cVal, pVals[tri[1]]); + s[2] = isoFraction(cVal, pVals[tri[2]]); + + if + ( + (s[0] >= 0.0 && s[0] <= 0.5) + && (s[1] >= 0.0 && s[1] <= 0.5) + && (s[2] >= 0.0 && s[2] <= 0.5) + ) + { + localTris.append + ( + labelledTri + ( + pointToLocal[tri[0]], + pointToLocal[tri[1]], + pointToLocal[tri[2]], + 0 + ) + ); + } + } + } + + pointField surfPoints; + surfPoints.transfer(localPoints); + pointIndexHit info = collapseSurface(surfPoints, localTris); + + if (info.hit()) + { + snappedCc[cellI] = snappedPoints.size(); + snappedPoints.append(info.hitPoint()); + + //Pout<< "cell:" << cellI + // << " at " << mesh_.cellCentres()[cellI] + // << " collapsing " << surfPoints + // << " intersections down to " + // << snappedPoints[snappedCc[cellI]] << endl; + } + } + } + } +} + + +// Generate triangles for face connected to pointI +void Foam::isoSurfaceCell::genPointTris +( + const scalarField& cellValues, + const scalarField& pointValues, + const label pointI, + const face& f, + const label cellI, + DynamicList<point, 64>& localTriPoints +) const +{ + const pointField& cc = mesh_.cellCentres(); + const pointField& pts = mesh_.points(); + + for (label fp = 1; fp < f.size() - 1; fp++) + { + triFace tri(f[0], f[fp], f[f.fcIndex(fp)]); + //List<triFace> tris(triangulate(f)); + //forAll(tris, i) + //{ + // const triFace& tri = tris[i]; + + label index = findIndex(tri, pointI); + + if (index == -1) + { + continue; + } + + // Tet between index..index-1, index..index+1, index..cc + label b = tri[tri.fcIndex(index)]; + label c = tri[tri.rcIndex(index)]; + + // Get fractions for the three edges emanating from point + FixedList<scalar, 3> s(3); + s[0] = isoFraction(pointValues[pointI], pointValues[b]); + s[1] = isoFraction(pointValues[pointI], pointValues[c]); + s[2] = isoFraction(pointValues[pointI], cellValues[cellI]); + + if + ( + (s[0] >= 0.0 && s[0] <= 0.5) + && (s[1] >= 0.0 && s[1] <= 0.5) + && (s[2] >= 0.0 && s[2] <= 0.5) + ) + { + localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[b]); + localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[c]); + localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*cc[cellI]); + } + } +} + + +// Generate triangle for tet connected to pointI +void Foam::isoSurfaceCell::genPointTris +( + const scalarField& pointValues, + const label pointI, + const label cellI, + DynamicList<point, 64>& localTriPoints +) const +{ + const pointField& pts = mesh_.points(); + const cell& cFaces = mesh_.cells()[cellI]; + + FixedList<label, 4> tet; + + label face0 = cFaces[0]; + const face& f0 = mesh_.faces()[face0]; + + if (mesh_.faceOwner()[face0] != cellI) + { + tet[0] = f0[0]; + tet[1] = f0[1]; + tet[2] = f0[2]; + } + else + { + tet[0] = f0[0]; + tet[1] = f0[2]; + tet[2] = f0[1]; + } + + // Find the point on the next face that is not on f0 + const face& f1 = mesh_.faces()[cFaces[1]]; + + forAll(f1, fp) + { + label p1 = f1[fp]; + + if (p1 != tet[0] && p1 != tet[1] && p1 != tet[2]) + { + tet[3] = p1; + break; + } + } + + + // Get the index of pointI + + label i0 = findIndex(tet, pointI); + label i1 = tet.fcIndex(i0); + label i2 = tet.fcIndex(i1); + label i3 = tet.fcIndex(i2); + + // Get fractions for the three edges emanating from point + FixedList<scalar, 3> s(3); + s[0] = isoFraction(pointValues[pointI], pointValues[tet[i1]]); + s[1] = isoFraction(pointValues[pointI], pointValues[tet[i2]]); + s[2] = isoFraction(pointValues[pointI], pointValues[tet[i3]]); + + if + ( + (s[0] >= 0.0 && s[0] <= 0.5) + && (s[1] >= 0.0 && s[1] <= 0.5) + && (s[2] >= 0.0 && s[2] <= 0.5) + ) + { + localTriPoints.append((1.0-s[0])*pts[pointI] + s[0]*pts[tet[i1]]); + localTriPoints.append((1.0-s[1])*pts[pointI] + s[1]*pts[tet[i2]]); + localTriPoints.append((1.0-s[2])*pts[pointI] + s[2]*pts[tet[i3]]); + } +} + + +void Foam::isoSurfaceCell::calcSnappedPoint +( + const PackedList<1>& isBoundaryPoint, + const PackedList<1>& isTet, + const scalarField& cVals, + const scalarField& pVals, + + DynamicList<point>& snappedPoints, + labelList& snappedPoint +) const +{ + const point greatPoint(VGREAT, VGREAT, VGREAT); + pointField collapsedPoint(mesh_.nPoints(), greatPoint); + + + // Work arrays + DynamicList<point, 64> localTriPoints(100); + labelHashSet localPointCells(100); + + forAll(mesh_.pointFaces(), pointI) + { + if (isBoundaryPoint.get(pointI) == 1) + { + continue; + } + + const labelList& pFaces = mesh_.pointFaces()[pointI]; + + bool anyCut = false; + + forAll(pFaces, i) + { + label faceI = pFaces[i]; + + if + ( + cellCutType_[mesh_.faceOwner()[faceI]] == CUT + || ( + mesh_.isInternalFace(faceI) + && cellCutType_[mesh_.faceNeighbour()[faceI]] == CUT + ) + ) + { + anyCut = true; + break; + } + } + + if (!anyCut) + { + continue; + } + + + localPointCells.clear(); + localTriPoints.clear(); + + forAll(pFaces, pFaceI) + { + label faceI = pFaces[pFaceI]; + const face& f = mesh_.faces()[faceI]; + label own = mesh_.faceOwner()[faceI]; + + // Triangulate around f[0] on owner side + if (isTet.get(own) == 1) + { + if (localPointCells.insert(own)) + { + genPointTris(pVals, pointI, own, localTriPoints); + } + } + else + { + genPointTris(cVals, pVals, pointI, f, own, localTriPoints); + } + + if (mesh_.isInternalFace(faceI)) + { + label nei = mesh_.faceNeighbour()[faceI]; + + if (isTet.get(nei) == 1) + { + if (localPointCells.insert(nei)) + { + genPointTris(pVals, pointI, nei, localTriPoints); + } + } + else + { + genPointTris(cVals, pVals, pointI, f, nei, localTriPoints); + } + } + } + + if (localTriPoints.size() == 0) + { + // No near intersections + } + else if (localTriPoints.size() == 3) + { + // Single triangle. No need for any analysis. Average points. + pointField points; + points.transfer(localTriPoints); + collapsedPoint[pointI] = sum(points)/points.size(); + + //Pout<< " point:" << pointI + // << " replacing coord:" << mesh_.points()[pointI] + // << " by average:" << collapsedPoint[pointI] << endl; + } + else + { + // Convert points into triSurface. + + // Merge points and compact out non-valid triangles + labelList triMap; // merged to unmerged triangle + labelList triPointReverseMap; // unmerged to merged point + triSurface surf + ( + stitchTriPoints + ( + false, // do not check for duplicate tris + localTriPoints, + triPointReverseMap, + triMap + ) + ); + + labelList faceZone; + label nZones = surf.markZones + ( + boolList(surf.nEdges(), false), + faceZone + ); + + if (nZones == 1) + { + collapsedPoint[pointI] = calcCentre(surf); + //Pout<< " point:" << pointI << " nZones:" << nZones + // << " replacing coord:" << mesh_.points()[pointI] + // << " by average:" << collapsedPoint[pointI] << endl; + } + } + } + + syncTools::syncPointList + ( + mesh_, + collapsedPoint, + minEqOp<point>(), + greatPoint, + true // are coordinates so separate + ); + + snappedPoint.setSize(mesh_.nPoints()); + snappedPoint = -1; + + forAll(collapsedPoint, pointI) + { + if (collapsedPoint[pointI] != greatPoint) + { + snappedPoint[pointI] = snappedPoints.size(); + snappedPoints.append(collapsedPoint[pointI]); + } + } +} + + + + +Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints +( + const bool checkDuplicates, + const List<point>& triPoints, + labelList& triPointReverseMap, // unmerged to merged point + labelList& triMap // merged to unmerged triangle +) const +{ + label nTris = triPoints.size()/3; + + if ((triPoints.size() % 3) != 0) + { + FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)") + << "Problem: number of points " << triPoints.size() + << " not a multiple of 3." << abort(FatalError); + } + + pointField newPoints; + mergePoints + ( + triPoints, + mergeDistance_, + false, + triPointReverseMap, + newPoints + ); + + // Check that enough merged. + if (debug) + { + pointField newNewPoints; + labelList oldToNew; + bool hasMerged = mergePoints + ( + newPoints, + mergeDistance_, + true, + oldToNew, + newNewPoints + ); + + if (hasMerged) + { + FatalErrorIn("isoSurfaceCell::stitchTriPoints(..)") + << "Merged points contain duplicates" + << " when merging with distance " << mergeDistance_ << endl + << "merged:" << newPoints.size() << " re-merged:" + << newNewPoints.size() + << abort(FatalError); + } + } + + + List<labelledTri> tris; + { + DynamicList<labelledTri> dynTris(nTris); + label rawPointI = 0; + DynamicList<label> newToOldTri(nTris); + + for (label oldTriI = 0; oldTriI < nTris; oldTriI++) + { + labelledTri tri + ( + triPointReverseMap[rawPointI], + triPointReverseMap[rawPointI+1], + triPointReverseMap[rawPointI+2], + 0 + ); + rawPointI += 3; + + if ((tri[0] != tri[1]) && (tri[0] != tri[2]) && (tri[1] != tri[2])) + { + newToOldTri.append(oldTriI); + dynTris.append(tri); + } + } + + triMap.transfer(newToOldTri.shrink()); + tris.transfer(dynTris.shrink()); + } + + + // Use face centres to determine 'flat hole' situation (see RMT paper). + // Two unconnected triangles get connected because (some of) the edges + // separating them get collapsed. Below only checks for duplicate triangles, + // not non-manifold edge connectivity. + if (checkDuplicates) + { + if (debug) + { + Pout<< "isoSurfaceCell : merged from " << nTris + << " down to " << tris.size() << " triangles." << endl; + } + + pointField centres(tris.size()); + forAll(tris, triI) + { + centres[triI] = tris[triI].centre(newPoints); + } + + pointField mergedCentres; + labelList oldToMerged; + bool hasMerged = mergePoints + ( + centres, + mergeDistance_, + false, + oldToMerged, + mergedCentres + ); + + if (debug) + { + Pout<< "isoSurfaceCell : detected " + << centres.size()-mergedCentres.size() + << " duplicate triangles." << endl; + } + + if (hasMerged) + { + // Filter out duplicates. + label newTriI = 0; + DynamicList<label> newToOldTri(tris.size()); + labelList newToMaster(mergedCentres.size(), -1); + forAll(tris, triI) + { + label mergedI = oldToMerged[triI]; + + if (newToMaster[mergedI] == -1) + { + newToMaster[mergedI] = triI; + newToOldTri.append(triMap[triI]); + tris[newTriI++] = tris[triI]; + } + } + + triMap.transfer(newToOldTri.shrink()); + tris.setSize(newTriI); + } + } + + return triSurface(tris, geometricSurfacePatchList(0), newPoints, true); +} + + +// Does face use valid vertices? +bool Foam::isoSurfaceCell::validTri(const triSurface& surf, const label faceI) +{ + // Simple check on indices ok. + + const labelledTri& f = surf[faceI]; + + if + ( + (f[0] < 0) || (f[0] >= surf.points().size()) + || (f[1] < 0) || (f[1] >= surf.points().size()) + || (f[2] < 0) || (f[2] >= surf.points().size()) + ) + { + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI << " vertices " << f + << " uses point indices outside point range 0.." + << surf.points().size()-1 << endl; + + return false; + } + + if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2])) + { + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI + << " uses non-unique vertices " << f + << endl; + return false; + } + + // duplicate triangle check + + const labelList& fFaces = surf.faceFaces()[faceI]; + + // Check if faceNeighbours use same points as this face. + // Note: discards normal information - sides of baffle are merged. + forAll(fFaces, i) + { + label nbrFaceI = fFaces[i]; + + if (nbrFaceI <= faceI) + { + // lower numbered faces already checked + continue; + } + + const labelledTri& nbrF = surf[nbrFaceI]; + + if + ( + ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2])) + && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2])) + && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2])) + ) + { + WarningIn("validTri(const triSurface&, const label)") + << "triangle " << faceI << " vertices " << f + << " has the same vertices as triangle " << nbrFaceI + << " vertices " << nbrF + << endl; + + return false; + } + } + return true; +} + + +void Foam::isoSurfaceCell::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<< "isoSurfaceCell : detected " + << mergedCentres.size() + << " 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(); + + forAll(oldToMerged, oldEdgeI) + { + label triI = oldEdgeI / 3; + label edgeI = oldToMerged[oldEdgeI]; + + if (edgeFace0[edgeI] == -1) + { + edgeFace0[edgeI] = triI; + } + else if (edgeFace1[edgeI] == -1) + { + edgeFace1[edgeI] = triI; + } + else + { + //WarningIn("orientSurface(triSurface&)") + // << "Edge " << edgeI << " with centre " << mergedCentres[edgeI] + // << " used by more than two triangles: " << edgeFace0[edgeI] + // << ", " + // << edgeFace1[edgeI] << " and " << triI << endl; + Map<labelList>::iterator iter = edgeFacesRest.find(edgeI); + + if (iter != edgeFacesRest.end()) + { + labelList& eFaces = iter(); + label sz = eFaces.size(); + eFaces.setSize(sz+1); + eFaces[sz] = triI; + } + else + { + edgeFacesRest.insert(edgeI, labelList(1, triI)); + } + } + } +} + + +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() > 0) + { + 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) + { + FatalErrorIn + ( + "isoSurfaceCell::orientSurface(triSurface&, const label)" + ) << "problem" << abort(FatalError); + } + } +} + + +// Checks if triangle is connected through edgeI only. +bool Foam::isoSurfaceCell::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::isoSurfaceCell::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::isoSurfaceCell::subsetMesh +( + const triSurface& s, + const labelList& newToOldFaces, + labelList& oldToNewPoints, + labelList& newToOldPoints +) +{ + const boolList include + ( + createWithValues<boolList> + ( + s.size(), + false, + newToOldFaces, + true + ) + ); + + newToOldPoints.setSize(s.points().size()); + oldToNewPoints.setSize(s.points().size()); + oldToNewPoints = -1; + { + label pointI = 0; + + forAll(include, oldFacei) + { + if (include[oldFacei]) + { + // Renumber labels for triangle + const labelledTri& tri = s[oldFacei]; + + forAll(tri, fp) + { + label oldPointI = tri[fp]; + + if (oldToNewPoints[oldPointI] == -1) + { + oldToNewPoints[oldPointI] = pointI; + newToOldPoints[pointI++] = oldPointI; + } + } + } + } + newToOldPoints.setSize(pointI); + } + + // Extract points + pointField newPoints(newToOldPoints.size()); + forAll(newToOldPoints, i) + { + newPoints[i] = s.points()[newToOldPoints[i]]; + } + // Extract faces + List<labelledTri> newTriangles(newToOldFaces.size()); + + forAll(newToOldFaces, i) + { + // Get old vertex labels + const labelledTri& tri = s[newToOldFaces[i]]; + + newTriangles[i][0] = oldToNewPoints[tri[0]]; + newTriangles[i][1] = oldToNewPoints[tri[1]]; + newTriangles[i][2] = oldToNewPoints[tri[2]]; + newTriangles[i].region() = tri.region(); + } + + // Reuse storage. + return triSurface(newTriangles, s.patches(), newPoints, true); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::isoSurfaceCell::isoSurfaceCell +( + const polyMesh& mesh, + const scalarField& cVals, + const scalarField& pVals, + const scalar iso, + const bool regularise, + const scalar mergeTol +) +: + mesh_(mesh), + iso_(iso), + mergeDistance_(mergeTol*mag(mesh.bounds().max()-mesh.bounds().min())) +{ + // Determine if cell is tet + PackedList<1> isTet(mesh_.nCells()); + { + tetMatcher tet; + + forAll(isTet, cellI) + { + if (tet.isA(mesh_, cellI)) + { + isTet.set(cellI, 1); + } + } + } + + // Determine if point is on boundary. Points on boundaries are never + // snapped. Coupled boundaries are handled explicitly so not marked here. + PackedList<1> isBoundaryPoint(mesh_.nPoints()); + const polyBoundaryMesh& patches = mesh_.boundaryMesh(); + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + + if (!pp.coupled()) + { + label faceI = pp.start(); + forAll(pp, i) + { + const face& f = mesh_.faces()[faceI++]; + + forAll(f, fp) + { + isBoundaryPoint.set(f[fp], 1); + } + } + } + } + + + // Determine if any cut through cell + calcCutTypes(isTet, cVals, pVals); + + DynamicList<point> snappedPoints(nCutCells_); + + // Per cc -1 or a point inside snappedPoints. + labelList snappedCc; + if (regularise) + { + calcSnappedCc + ( + isTet, + cVals, + pVals, + snappedPoints, + snappedCc + ); + } + else + { + snappedCc.setSize(mesh_.nCells()); + snappedCc = -1; + } + + snappedPoints.shrink(); + + if (debug) + { + Pout<< "isoSurfaceCell : shifted " << snappedPoints.size() + << " cell centres to intersection." << endl; + } + + label nCellSnaps = snappedPoints.size(); + + // Per point -1 or a point inside snappedPoints. + labelList snappedPoint; + if (regularise) + { + calcSnappedPoint + ( + isBoundaryPoint, + isTet, + cVals, + pVals, + snappedPoints, + snappedPoint + ); + } + else + { + snappedPoint.setSize(mesh_.nPoints()); + snappedPoint = -1; + } + + if (debug) + { + Pout<< "isoSurfaceCell : shifted " << snappedPoints.size()-nCellSnaps + << " vertices to intersection." << endl; + } + + + + DynamicList<point> triPoints(nCutCells_); + DynamicList<label> triMeshCells(nCutCells_); + + generateTriPoints + ( + cVals, + pVals, + + mesh_.cellCentres(), + mesh_.points(), + + snappedPoints, + snappedCc, + snappedPoint, + + triPoints, + triMeshCells + ); + + 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 + ( + true, // check for duplicate tris + triPoints, + triPointMergeMap_, // unmerged to merged point + triMap + ) + ); + + if (debug) + { + Pout<< "isoSurfaceCell : generated " << triMap.size() + << " merged triangles." << endl; + } + + meshCells_.setSize(triMap.size()); + forAll(triMap, i) + { + meshCells_[i] = triMeshCells[triMap[i]]; + } + + if (debug) + { + Pout<< "isoSurfaceCell : checking " << size() + << " triangles for validity." << endl; + + forAll(*this, triI) + { + // Copied from surfaceCheck + validTri(*this, triI); + } + } + + + 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<< "isoSurfaceCell : 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); + + //combineCellTriangles(); +} + + +////XXXXXXX +//// Experimental retriangulation of triangles per cell. Problem is that +//// -it is very expensive -only gets rid of internal points, not of boundary +//// ones so limited benefit (e.g. 60 v.s. 88 triangles) +//void Foam::isoSurfaceCell::combineCellTriangles() +//{ +// if (size() > 0) +// { +// DynamicList<labelledTri> newTris(size()); +// DynamicList<label> newTriToCell(size()); +// +// label startTriI = 0; +// +// DynamicList<labelledTri> tris; +// +// for (label triI = 1; triI <= meshCells_.size(); triI++) +// { +// if +// ( +// triI == meshCells_.size() +// || meshCells_[triI] != meshCells_[startTriI] +// ) +// { +// label nTris = triI-startTriI; +// +// if (nTris == 1) +// { +// newTris.append(operator[](startTriI)); +// newTriToCell.append(meshCells_[startTriI]); +// } +// else +// { +// // Collect from startTriI to triI in a triSurface +// tris.clear(); +// for (label i = startTriI; i < triI; i++) +// { +// tris.append(operator[](i)); +// } +// triSurface cellTris(tris, patches(), points()); +// tris.clear(); +// +// // Get outside +// const labelListList& loops = cellTris.edgeLoops(); +// +// forAll(loops, i) +// { +// // Do proper triangulation of loop +// face loop(renumber(cellTris.meshPoints(), loops[i])); +// +// faceTriangulation faceTris +// ( +// points(), +// loop, +// true +// ); +// +// // Copy into newTris +// forAll(faceTris, faceTriI) +// { +// const triFace& tri = faceTris[faceTriI]; +// +// newTris.append +// ( +// labelledTri +// ( +// tri[0], +// tri[1], +// tri[2], +// operator[](startTriI).region() +// ) +// ); +// newTriToCell.append(meshCells_[startTriI]); +// } +// } +// } +// +// startTriI = triI; +// } +// } +// newTris.shrink(); +// newTriToCell.shrink(); +// +// // Compact +// pointField newPoints(points().size()); +// label newPointI = 0; +// labelList oldToNewPoint(points().size(), -1); +// +// forAll(newTris, i) +// { +// labelledTri& tri = newTris[i]; +// forAll(tri, j) +// { +// label pointI = tri[j]; +// +// if (oldToNewPoint[pointI] == -1) +// { +// oldToNewPoint[pointI] = newPointI; +// newPoints[newPointI++] = points()[pointI]; +// } +// tri[j] = oldToNewPoint[pointI]; +// } +// } +// newPoints.setSize(newPointI); +// +// triSurface::operator= +// ( +// triSurface +// ( +// newTris, +// patches(), +// newPoints, +// true +// ) +// ); +// meshCells_.transfer(newTriToCell); +// } +//} +////XXXXXXX + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H new file mode 100644 index 0000000000000000000000000000000000000000..40e2b6d79378ab15f5552616785c7528717a3d5e --- /dev/null +++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H @@ -0,0 +1,372 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::isoSurfaceCell + +Description + A surface formed by the iso value. + After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke + and + "Regularised Marching Tetrahedra: improved iso-surface extraction", + G.M. Treece, R.W. Prager and A.H. Gee. + + See isoSurface. This is a variant which does tetrahedrisation from + triangulation of face to cell centre instead of edge of face to two + neighbouring cell centres. This gives much lower quality triangles + but they are local to a cell. + +SourceFiles + isoSurfaceCell.C + +\*---------------------------------------------------------------------------*/ + +#ifndef isoSurfaceCell_H +#define isoSurfaceCell_H + +#include "triSurface.H" +#include "labelPair.H" +#include "pointIndexHit.H" +#include "PackedList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class polyMesh; + +/*---------------------------------------------------------------------------*\ + Class isoSurfaceCell Declaration +\*---------------------------------------------------------------------------*/ + +class isoSurfaceCell +: + public triSurface +{ + // Private data + + enum segmentCutType + { + NEARFIRST, // intersection close to e.first() + NEARSECOND, // ,, e.second() + NOTNEAR // no intersection + }; + + enum cellCutType + { + NOTCUT, // not cut + SPHERE, // all edges to cell centre cut + CUT // normal cut + }; + + + //- Reference to mesh + const polyMesh& mesh_; + + //- isoSurfaceCell value + const scalar iso_; + + //- When to merge points + const scalar mergeDistance_; + + //- Whether cell might be cut + List<cellCutType> cellCutType_; + + //- Estimated number of cut cells + label nCutCells_; + + //- For every triangle the original cell in mesh + labelList meshCells_; + + //- For every unmerged triangle point the point in the triSurface + labelList triPointMergeMap_; + + + // Private Member Functions + + //- Get location of iso value as fraction inbetween s0,s1 + scalar isoFraction + ( + const scalar s0, + const scalar s1 + ) const; + + //List<triFace> triangulate(const face& f) const; + + //- Determine whether cell is cut + cellCutType calcCutType + ( + const PackedList<1>&, + const scalarField& cellValues, + const scalarField& pointValues, + const label + ) const; + + void calcCutTypes + ( + const PackedList<1>&, + const scalarField& cellValues, + const scalarField& pointValues + ); + + 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 + ( + const PackedList<1>& isTet, + const scalarField& cVals, + const scalarField& pVals, + DynamicList<point>& triPoints, + labelList& snappedCc + ) const; + + //- Generate triangles for face connected to pointI + void genPointTris + ( + const scalarField& cellValues, + const scalarField& pointValues, + const label pointI, + const face& f, + const label cellI, + DynamicList<point, 64>& localTriPoints + ) const; + + //- Generate triangles for tet connected to pointI + void genPointTris + ( + const scalarField& pointValues, + const label pointI, + const label cellI, + DynamicList<point, 64>& localTriPoints + ) const; + + //- Determine per point whether all near cuts can be snapped to single + // point. + void calcSnappedPoint + ( + const PackedList<1>& isBoundaryPoint, + const PackedList<1>& isTet, + const scalarField& cVals, + const scalarField& pVals, + DynamicList<point>& triPoints, + labelList& snappedPoint + ) const; + + //- Generate single point by interpolation or snapping + template<class Type> + Type generatePoint + ( + const DynamicList<Type>& snappedPoints, + const scalar s0, + const Type& p0, + const label p0Index, + const scalar s1, + const Type& p1, + const label p1Index + ) const; + + template<class Type> + void generateTriPoints + ( + const DynamicList<Type>& snapped, + const scalar s0, + const Type& p0, + const label p0Index, + const scalar s1, + const Type& p1, + const label p1Index, + const scalar s2, + const Type& p2, + const label p2Index, + const scalar s3, + const Type& p3, + const label p3Index, + DynamicList<Type>& points + ) const; + + template<class Type> + void generateTriPoints + ( + const scalarField& cVals, + const scalarField& pVals, + + const Field<Type>& cCoords, + const Field<Type>& pCoords, + + const DynamicList<Type>& snappedPoints, + const labelList& snappedCc, + const labelList& snappedPoint, + + DynamicList<Type>& triPoints, + DynamicList<label>& triMeshCells + ) const; + + triSurface stitchTriPoints + ( + const bool checkDuplicates, + const List<point>& triPoints, + labelList& triPointReverseMap, // unmerged to merged point + labelList& triMap // merged to unmerged triangle + ) const; + + //- 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, + const labelList& newToOldFaces, + labelList& oldToNewPoints, + labelList& newToOldPoints + ); + + //- Combine all triangles inside a cell into a minimal triangulation + void combineCellTriangles(); + +public: + + //- Runtime type information + TypeName("isoSurfaceCell"); + + + // Constructors + + //- Construct from dictionary + isoSurfaceCell + ( + const polyMesh& mesh, + const scalarField& cellValues, + const scalarField& pointValues, + const scalar iso, + const bool regularise, + const scalar mergeTol = 1E-6 // fraction of bounding box + ); + + + // Member Functions + + //- For every face 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. 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; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "isoSurfaceCellTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..748e1a3311dbc19e8e26ba0dea9d1b38ddd07faf --- /dev/null +++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C @@ -0,0 +1,384 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "isoSurfaceCell.H" +#include "polyMesh.H" +#include "tetMatcher.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template<class Type> +Type Foam::isoSurfaceCell::generatePoint +( + const DynamicList<Type>& snappedPoints, + + const scalar s0, + const Type& p0, + const label p0Index, + + const scalar s1, + const Type& p1, + const label p1Index +) const +{ + scalar d = s1-s0; + + if (mag(d) > VSMALL) + { + scalar s = (iso_-s0)/d; + + if (s >= 0.5 && s <= 1 && p1Index != -1) + { + return snappedPoints[p1Index]; + } + else if (s >= 0.0 && s <= 0.5 && p0Index != -1) + { + return snappedPoints[p0Index]; + } + else + { + return s*p1 + (1.0-s)*p0; + } + } + else + { + scalar s = 0.4999; + + return s*p1 + (1.0-s)*p0; + } +} + + +template<class Type> +void Foam::isoSurfaceCell::generateTriPoints +( + const DynamicList<Type>& snapped, + + const scalar s0, + const Type& p0, + const label p0Index, + + const scalar s1, + const Type& p1, + const label p1Index, + + const scalar s2, + const Type& p2, + const label p2Index, + + const scalar s3, + const Type& p3, + const label p3Index, + + DynamicList<Type>& points +) const +{ + int triIndex = 0; + if (s0 < iso_) + { + triIndex |= 1; + } + if (s1 < iso_) + { + triIndex |= 2; + } + if (s2 < iso_) + { + triIndex |= 4; + } + if (s3 < iso_) + { + triIndex |= 8; + } + + /* Form the vertices of the triangles for each case */ + switch (triIndex) + { + case 0x00: + case 0x0F: + break; + + case 0x0E: + case 0x01: + points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index)); + points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)); + points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); + break; + + case 0x0D: + case 0x02: + points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index)); + points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)); + points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)); + break; + + case 0x0C: + case 0x03: + { + Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index); + Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index); + + points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); + points.append(tp1); + points.append(tp2); + points.append(tp2); + points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)); + points.append(tp1); + } + break; + + case 0x0B: + case 0x04: + { + points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index)); + points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)); + points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)); + } + break; + + case 0x0A: + case 0x05: + { + Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index); + Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index); + + points.append(tp0); + points.append(tp1); + points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); + points.append(tp0); + points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)); + points.append(tp1); + } + break; + + case 0x09: + case 0x06: + { + Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index); + Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index); + + points.append(tp0); + points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)); + points.append(tp1); + points.append(tp0); + points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)); + points.append(tp1); + } + break; + + case 0x07: + case 0x08: + points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index)); + points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)); + points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)); + break; + } +} + + +template<class Type> +void Foam::isoSurfaceCell::generateTriPoints +( + const scalarField& cVals, + const scalarField& pVals, + + const Field<Type>& cCoords, + const Field<Type>& pCoords, + + const DynamicList<Type>& snappedPoints, + const labelList& snappedCc, + const labelList& snappedPoint, + + DynamicList<Type>& triPoints, + DynamicList<label>& triMeshCells +) const +{ + tetMatcher tet; + + forAll(mesh_.cells(), cellI) + { + if (cellCutType_[cellI] != NOTCUT) + { + label oldNPoints = triPoints.size(); + + const cell& cFaces = mesh_.cells()[cellI]; + + if (tet.isA(mesh_, cellI)) + { + // For tets don't do cell-centre decomposition, just use the + // tet points and values + + const face& f0 = mesh_.faces()[cFaces[0]]; + + // Get the other point + const face& f1 = mesh_.faces()[cFaces[1]]; + label oppositeI = -1; + forAll(f1, fp) + { + oppositeI = f1[fp]; + + if (findIndex(f0, oppositeI) == -1) + { + break; + } + } + + generateTriPoints + ( + snappedPoints, + pVals[f0[0]], + pCoords[f0[0]], + snappedPoint[f0[0]], + + pVals[f0[1]], + pCoords[f0[1]], + snappedPoint[f0[1]], + + pVals[f0[2]], + pCoords[f0[2]], + snappedPoint[f0[2]], + + pVals[oppositeI], + pCoords[oppositeI], + snappedPoint[oppositeI], + + triPoints + ); + } + else + { + const cell& cFaces = mesh_.cells()[cellI]; + + forAll(cFaces, cFaceI) + { + label faceI = cFaces[cFaceI]; + const face& f = mesh_.faces()[faceI]; + + for (label fp = 1; fp < f.size() - 1; fp++) + { + triFace tri(f[0], f[fp], f[f.fcIndex(fp)]); + //List<triFace> tris(triangulate(f)); + //forAll(tris, i) + //{ + // const triFace& tri = tris[i]; + + + generateTriPoints + ( + snappedPoints, + + pVals[tri[0]], + pCoords[tri[0]], + snappedPoint[tri[0]], + + pVals[tri[1]], + pCoords[tri[1]], + snappedPoint[tri[1]], + + pVals[tri[2]], + pCoords[tri[2]], + snappedPoint[tri[2]], + + cVals[cellI], + cCoords[cellI], + snappedCc[cellI], + + triPoints + ); + } + } + } + + + // Every three triPoints is a cell + label nCells = (triPoints.size()-oldNPoints)/3; + for (label i = 0; i < nCells; i++) + { + triMeshCells.append(cellI); + } + } + } + + triPoints.shrink(); + triMeshCells.shrink(); +} + + +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; +} + + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C index e1500ca6e5e04459445da1bb42be8ad9ffa94740..b4d440c4fde3059145f3b69a7173a330b0ee6b3f 100644 --- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C +++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C @@ -26,69 +26,90 @@ License #include "isoSurface.H" #include "polyMesh.H" -#include "tetMatcher.H" +#include "syncTools.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template<class Type> -Type Foam::isoSurface::vertexInterp +Type Foam::isoSurface::generatePoint ( - const scalar iso, + const DynamicList<Type>& snappedPoints, + + const scalar s0, const Type& p0, + const label p0Index, + + const scalar s1, const Type& p1, - const scalar s0, - const scalar s1 -) + const label p1Index +) const { scalar d = s1-s0; if (mag(d) > VSMALL) { - return (iso-s0)/d*p1 + (s1-iso)/d*p0; + scalar s = (iso_-s0)/d; + + if (s >= 0.5 && s <= 1 && p1Index != -1) + { + return snappedPoints[p1Index]; + } + else if (s >= 0.0 && s <= 0.5 && p0Index != -1) + { + return snappedPoints[p0Index]; + } + else + { + return s*p1 + (1.0-s)*p0; + } } else { - return 0.5*(p0+p1); + scalar s = 0.4999; + + return s*p1 + (1.0-s)*p0; } } -// After "Polygonising A Scalar Field Using Tetrahedrons" -// by Paul Bourke -// Get value consistent with uncompacted triangle points. -// Given tet corner sample values s0..s3 interpolate the corresponding -// values p0..p3 to construct the surface corresponding to sample value iso. template<class Type> -void Foam::isoSurface::vertexInterp +void Foam::isoSurface::generateTriPoints ( - const scalar iso, - const scalar s0, - const scalar s1, - const scalar s2, - const scalar s3, + const DynamicList<Type>& snapped, + const scalar s0, const Type& p0, + const label p0Index, + + const scalar s1, const Type& p1, + const label p1Index, + + const scalar s2, const Type& p2, + const label p2Index, + + const scalar s3, const Type& p3, + const label p3Index, DynamicList<Type>& points -) +) const { int triIndex = 0; - if (s0 < iso) + if (s0 < iso_) { triIndex |= 1; } - if (s1 < iso) + if (s1 < iso_) { triIndex |= 2; } - if (s2 < iso) + if (s2 < iso_) { triIndex |= 4; } - if (s3 < iso) + if (s3 < iso_) { triIndex |= 8; } @@ -102,29 +123,29 @@ void Foam::isoSurface::vertexInterp case 0x0E: case 0x01: - points.append(vertexInterp(iso,p0,p1,s0,s1)); - points.append(vertexInterp(iso,p0,p2,s0,s2)); - points.append(vertexInterp(iso,p0,p3,s0,s3)); + points.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index)); + points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)); + points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); break; case 0x0D: case 0x02: - points.append(vertexInterp(iso,p1,p0,s1,s0)); - points.append(vertexInterp(iso,p1,p3,s1,s3)); - points.append(vertexInterp(iso,p1,p2,s1,s2)); + points.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index)); + points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)); + points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)); break; case 0x0C: case 0x03: { - const Type tp1 = vertexInterp(iso,p0,p2,s0,s2); - const Type tp2 = vertexInterp(iso,p1,p3,s1,s3); + Type tp1 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index); + Type tp2 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index); - points.append(vertexInterp(iso,p0,p3,s0,s3)); + points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); points.append(tp1); points.append(tp2); points.append(tp2); - points.append(vertexInterp(iso,p1,p2,s1,s2)); + points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)); points.append(tp1); } break; @@ -132,23 +153,23 @@ void Foam::isoSurface::vertexInterp case 0x0B: case 0x04: { - points.append(vertexInterp(iso,p2,p0,s2,s0)); - points.append(vertexInterp(iso,p2,p1,s2,s1)); - points.append(vertexInterp(iso,p2,p3,s2,s3)); + points.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index)); + points.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)); + points.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)); } break; case 0x0A: case 0x05: { - const Type tp0 = vertexInterp(iso,p0,p1,s0,s1); - const Type tp1 = vertexInterp(iso,p2,p3,s2,s3); + Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index); + Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index); points.append(tp0); points.append(tp1); - points.append(vertexInterp(iso,p0,p3,s0,s3)); + points.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); points.append(tp0); - points.append(vertexInterp(iso,p1,p2,s1,s2)); + points.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)); points.append(tp1); } break; @@ -156,133 +177,305 @@ void Foam::isoSurface::vertexInterp case 0x09: case 0x06: { - const Type tp0 = vertexInterp(iso,p0,p1,s0,s1); - const Type tp1 = vertexInterp(iso,p2,p3,s2,s3); + Type tp0 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index); + Type tp1 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index); points.append(tp0); - points.append(vertexInterp(iso,p1,p3,s1,s3)); + points.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)); points.append(tp1); points.append(tp0); - points.append(vertexInterp(iso,p0,p2,s0,s2)); + points.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)); points.append(tp1); } break; case 0x07: case 0x08: - points.append(vertexInterp(iso,p3,p0,s3,s0)); - points.append(vertexInterp(iso,p3,p2,s3,s2)); - points.append(vertexInterp(iso,p3,p1,s3,s1)); + points.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index)); + points.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)); + points.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)); break; } } -template <class Type> -Foam::tmp<Foam::Field<Type> > -Foam::isoSurface::sample(const Field<Type>& vField) const +template<class Type> +Foam::label Foam::isoSurface::generateTriPoints +( + const volScalarField& cVals, + const scalarField& pVals, + + const GeometricField<Type, fvPatchField, volMesh>& cCoords, + const Field<Type>& pCoords, + + const DynamicList<Type>& snappedPoints, + const labelList& snappedCc, + const labelList& snappedPoint, + const label faceI, + + const scalar neiVal, + const Type& neiPt, + const label neiSnap, + + DynamicList<Type>& triPoints, + DynamicList<label>& triMeshCells +) const { - return tmp<Field<Type> >(new Field<Type>(vField, meshCells())); + label own = mesh_.faceOwner()[faceI]; + + label oldNPoints = triPoints.size(); + + const face& f = mesh_.faces()[faceI]; + + forAll(f, fp) + { + label pointI = f[fp]; + label nextPointI = f[f.fcIndex(fp)]; + + generateTriPoints + ( + snappedPoints, + + pVals[pointI], + pCoords[pointI], + snappedPoint[pointI], + + pVals[nextPointI], + pCoords[nextPointI], + snappedPoint[nextPointI], + + cVals[own], + cCoords[own], + snappedCc[own], + + neiVal, + neiPt, + neiSnap, + + triPoints + ); + } + + // Every three triPoints is a triangle + label nTris = (triPoints.size()-oldNPoints)/3; + for (label i = 0; i < nTris; i++) + { + triMeshCells.append(own); + } + + return nTris; } -template <class Type> -Foam::tmp<Foam::Field<Type> > -Foam::isoSurface::interpolate +template<class Type> +void Foam::isoSurface::generateTriPoints ( - const Field<Type>& sampleCellValues, - const Field<Type>& samplePointValues + const volScalarField& cVals, + const scalarField& pVals, + + const GeometricField<Type, fvPatchField, volMesh>& cCoords, + const Field<Type>& pCoords, + + const DynamicList<Type>& snappedPoints, + const labelList& snappedCc, + const labelList& snappedPoint, + + DynamicList<Type>& triPoints, + DynamicList<label>& triMeshCells ) const { - tetMatcher tet; + const polyBoundaryMesh& patches = mesh_.boundaryMesh(); + const labelList& own = mesh_.faceOwner(); + const labelList& nei = mesh_.faceNeighbour(); - DynamicList<Type> triValues; - // Note: in same order as construction of triSurface - label oldCellI = -1; - forAll(meshCells_, triI) + // Determine neighbouring snap status + labelList neiSnappedCc(mesh_.nFaces()-mesh_.nInternalFaces(), -1); + forAll(patches, patchI) { - label cellI = meshCells_[triI]; + const polyPatch& pp = patches[patchI]; - if (cellI != oldCellI) + if (pp.coupled()) { - oldCellI = cellI; + label faceI = pp.start(); + forAll(pp, i) + { + neiSnappedCc[faceI-mesh_.nInternalFaces()] = + snappedCc[own[faceI]]; + faceI++; + } + } + } + syncTools::swapBoundaryFaceList(mesh_, neiSnappedCc, false); - const cell& cFaces = mesh_.cells()[cellI]; - if (tet.isA(mesh_, cellI)) - { - // For tets don't do cell-centre decomposition, just use the - // tet points and values - const face& f0 = mesh_.faces()[cFaces[0]]; + // Generate triangle points - // Get the other point - const face& f1 = mesh_.faces()[cFaces[1]]; - label oppositeI = -1; - forAll(f1, fp) - { - oppositeI = f1[fp]; + triPoints.clear(); + triMeshCells.clear(); - if (findIndex(f0, oppositeI) == -1) - { - break; - } - } + for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) + { + if (faceCutType_[faceI] != NOTCUT) + { + generateTriPoints + ( + cVals, + pVals, + + cCoords, + pCoords, + + snappedPoints, + snappedCc, + snappedPoint, + faceI, + + cVals[nei[faceI]], + cCoords[nei[faceI]], + snappedCc[nei[faceI]], + + triPoints, + triMeshCells + ); + } + } - vertexInterp - ( - iso_, - pointValues_[f0[0]], - pointValues_[f0[1]], - pointValues_[f0[2]], - pointValues_[oppositeI], - - samplePointValues[f0[0]], - samplePointValues[f0[1]], - samplePointValues[f0[2]], - samplePointValues[oppositeI], - - triValues - ); - } - else + + forAll(patches, patchI) + { + const polyPatch& pp = patches[patchI]; + + if (pp.coupled()) + { + if (refCast<const processorPolyPatch>(pp).owner()) { - forAll(cFaces, cFaceI) - { - label faceI = cFaces[cFaceI]; - const face& f = mesh_.faces()[faceI]; + label faceI = pp.start(); - for(label fp = 1; fp < f.size() - 1; fp++) + forAll(pp, i) + { + if (faceCutType_[faceI] != NOTCUT) { - vertexInterp + generateTriPoints ( - iso_, - pointValues_[f[0]], - pointValues_[f[fp]], - pointValues_[f[f.fcIndex(fp)]], - cellValues_[cellI], - - samplePointValues[f[0]], - samplePointValues[f[fp]], - samplePointValues[f[f.fcIndex(fp)]], - sampleCellValues[cellI], - - triValues + cVals, + pVals, + + cCoords, + pCoords, + + snappedPoints, + snappedCc, + snappedPoint, + faceI, + + cVals.boundaryField()[patchI][i], + cCoords.boundaryField()[patchI][i], + neiSnappedCc[faceI-mesh_.nInternalFaces()], + + triPoints, + triMeshCells ); } + faceI++; } } } + else + { + label faceI = pp.start(); + + forAll(pp, i) + { + if (faceCutType_[faceI] != NOTCUT) + { + generateTriPoints + ( + cVals, + pVals, + + cCoords, + pCoords, + + snappedPoints, + snappedCc, + snappedPoint, + faceI, + + cVals.boundaryField()[patchI][i], + cCoords.boundaryField()[patchI][i], + -1, // fc not snapped + + triPoints, + triMeshCells + ); + } + faceI++; + } + } } + triPoints.shrink(); + triMeshCells.shrink(); +} + + +//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 volScalarField& cVals, + const scalarField& pVals, + const GeometricField<Type, fvPatchField, volMesh>& 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(triValues, i) + forAll(triPoints, i) { - values[triPointMergeMap_[i]] = triValues[i]; + label mergedPointI = triPointMergeMap_[i]; + + if (mergedPointI >= 0) + { + values[mergedPointI] = triPoints[i]; + } } return tvalues; diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C index a5ad5dc0edab86dd6bc1f7f21bfa4dfd20a834bf..782615107ba2498a486f8202a4a09a83c3cd4eb2 100644 --- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C @@ -41,124 +41,192 @@ namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::sampledIsoSurface::createGeometry() const +void Foam::sampledIsoSurface::getIsoFields() const { const fvMesh& fvm = static_cast<const fvMesh&>(mesh()); - if (fvm.time().timeIndex() != storedTimeIndex_) - { - storedTimeIndex_ = fvm.time().timeIndex(); + // Get volField + // ~~~~~~~~~~~~ - // Clear any stored topo - facesPtr_.clear(); - - // Optionally read volScalarField - autoPtr<volScalarField> readFieldPtr_; - - // 1. see if field in database - // 2. see if field can be read - const volScalarField* cellFldPtr = NULL; - if (fvm.foundObject<volScalarField>(isoField_)) + if (fvm.foundObject<volScalarField>(isoField_)) + { + if (debug) { - if (debug) - { - Info<< "sampledIsoSurface::createGeometry() : lookup " - << isoField_ << endl; - } - - cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_); + Info<< "sampledIsoSurface::getIsoField() : lookup " + << isoField_ << endl; } - else - { - // Bit of a hack. Read field and store. + storedVolFieldPtr_.clear(); + volFieldPtr_ = &fvm.lookupObject<volScalarField>(isoField_); + } + else + { + // Bit of a hack. Read field and store. - if (debug) - { - Info<< "sampledIsoSurface::createGeometry() : reading " - << isoField_ << " from time " <<fvm.time().timeName() - << endl; - } + if (debug) + { + Info<< "sampledIsoSurface::getIsoField() : reading " + << isoField_ << " from time " <<fvm.time().timeName() + << endl; + } - readFieldPtr_.reset + storedVolFieldPtr_.reset + ( + new volScalarField ( - new volScalarField + IOobject ( - IOobject - ( - isoField_, - fvm.time().timeName(), - fvm, - IOobject::MUST_READ, - IOobject::NO_WRITE, - false - ), - fvm - ) - ); + isoField_, + fvm.time().timeName(), + fvm, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ), + fvm + ) + ); + volFieldPtr_ = storedVolFieldPtr_.operator->(); + } - cellFldPtr = readFieldPtr_.operator->(); + + + // Get pointField + // ~~~~~~~~~~~~~~ + + word pointFldName = "volPointInterpolate(" + isoField_ + ')'; + + if (fvm.foundObject<pointScalarField>(pointFldName)) + { + if (debug) + { + Info<< "sampledIsoSurface::getIsoField() : lookup " + << pointFldName << endl; + } + storedPointFieldPtr_.clear(); + pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName); + } + else + { + // Not in registry. Interpolate. + + if (debug) + { + Info<< "sampledIsoSurface::getIsoField() : interpolating " + << pointFldName << endl; } - const volScalarField& cellFld = *cellFldPtr; - tmp<pointScalarField> pointFld + storedPointFieldPtr_.reset ( - volPointInterpolation::New(fvm).interpolate(cellFld) + volPointInterpolation::New(fvm).interpolate(*volFieldPtr_).ptr() ); + pointFieldPtr_ = storedPointFieldPtr_.operator->(); + } - //- Direct from cell field and point field. Gives bad continuity. - //const isoSurface iso - //( - // fvm, - // cellFld.internalField(), - // pointFld().internalField(), - // isoVal_ - //); - - //- From point field and interpolated cell. - scalarField cellAvg(fvm.nCells(), scalar(0.0)); - labelField nPointCells(fvm.nCells(), 0); + if (debug) + { + Info<< "sampledIsoSurface::getIsoField() : obtained volField " + << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value() + << " max:" << max(*volFieldPtr_).value() << endl; + Info<< "sampledIsoSurface::getIsoField() : obtained pointField " + << pointFieldPtr_->name() + << " min:" << gMin(pointFieldPtr_->internalField()) + << " max:" << gMax(pointFieldPtr_->internalField()) << endl; + } +} + + +void Foam::sampledIsoSurface::createGeometry() const +{ + const fvMesh& fvm = static_cast<const fvMesh&>(mesh()); + + if (fvm.time().timeIndex() != storedTimeIndex_) + { + storedTimeIndex_ = fvm.time().timeIndex(); + + getIsoFields(); + + // Clear any stored topo + surfPtr_.clear(); + facesPtr_.clear(); + + if (average_) { - for (label pointI = 0; pointI < fvm.nPoints(); pointI++) + //- From point field and interpolated cell. + volScalarField cellAvg + ( + IOobject + ( + "cellAvg", + fvm.time().timeName(), + fvm.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + fvm, + dimensionedScalar("zero", dimless, scalar(0.0)) + ); + labelField nPointCells(fvm.nCells(), 0); { - const labelList& pCells = fvm.pointCells(pointI); - - forAll(pCells, i) + for (label pointI = 0; pointI < fvm.nPoints(); pointI++) { - label cellI = pCells[i]; + const labelList& pCells = fvm.pointCells(pointI); + + forAll(pCells, i) + { + label cellI = pCells[i]; - cellAvg[cellI] += pointFld().internalField()[pointI]; - nPointCells[cellI]++; + cellAvg[cellI] += (*pointFieldPtr_)[pointI]; + nPointCells[cellI]++; + } } } + forAll(cellAvg, cellI) + { + cellAvg[cellI] /= nPointCells[cellI]; + } + // Give value to calculatedFvPatchFields + cellAvg.correctBoundaryConditions(); + + surfPtr_.reset + ( + new isoSurface + ( + cellAvg, + *pointFieldPtr_, + isoVal_, + regularise_ + ) + ); } - forAll(cellAvg, cellI) + else { - cellAvg[cellI] /= nPointCells[cellI]; + //- Direct from cell field and point field. + surfPtr_.reset + ( + new isoSurface + ( + *volFieldPtr_, + *pointFieldPtr_, + isoVal_, + regularise_ + ) + ); } - const isoSurface iso - ( - fvm, - cellAvg, - pointFld().internalField(), - isoVal_ - ); - - - const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso); - meshCells_ = iso.meshCells(); - triPointMergeMap_ = iso.triPointMergeMap(); if (debug) { Pout<< "sampledIsoSurface::createGeometry() : constructed iso:" << nl + << " regularise : " << regularise_ << nl + << " average : " << average_ << nl << " isoField : " << isoField_ << nl << " isoValue : " << isoVal_ << nl - << " unmerged points: " << triPointMergeMap_.size() << nl << " points : " << points().size() << nl - << " tris : " << triSurface::size() << nl - << " cut cells : " << meshCells_.size() << endl; + << " tris : " << surface().size() << nl + << " cut cells : " << surface().meshCells().size() + << endl; } } } @@ -176,12 +244,27 @@ Foam::sampledIsoSurface::sampledIsoSurface sampledSurface(name, mesh, dict), isoField_(dict.lookup("isoField")), isoVal_(readScalar(dict.lookup("isoValue"))), + regularise_(dict.lookupOrDefault("regularise", true)), + average_(dict.lookupOrDefault("average", false)), zoneName_(word::null), + surfPtr_(NULL), facesPtr_(NULL), storedTimeIndex_(-1), - meshCells_(0), - triPointMergeMap_(0) + storedVolFieldPtr_(NULL), + volFieldPtr_(NULL), + storedPointFieldPtr_(NULL), + pointFieldPtr_(NULL) { + if (!sampledSurface::interpolate()) + { + FatalErrorIn + ( + "sampledIsoSurface::sampledIsoSurface" + "(const word&, const polyMesh&, const dictionary&)" + ) << "Non-interpolated iso surface not supported since triangles" + << " span across cells." << exit(FatalError); + } + // label zoneId = -1; // if (dict.readIfPresent("zone", zoneName_)) // { @@ -209,6 +292,7 @@ void Foam::sampledIsoSurface::correct(const bool meshChanged) // Only change of mesh changes plane - zone restriction gets lost if (meshChanged) { + surfPtr_.clear(); facesPtr_.clear(); } } @@ -317,9 +401,9 @@ void Foam::sampledIsoSurface::print(Ostream& os) const { os << "sampledIsoSurface: " << name() << " :" << " field:" << isoField_ - << " value:" << isoVal_ - << " faces:" << faces().size() - << " points:" << points().size(); + << " value:" << isoVal_; + //<< " faces:" << faces().size() // note: possibly no geom yet + //<< " points:" << points().size(); } diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H index a55cb6641a22c750dc4f8329953cc1fc4e891433..3ccd72f629ab71c11e4d205247d42f3cd9c6b912 100644 --- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H @@ -39,7 +39,7 @@ SourceFiles #define sampledIsoSurface_H #include "sampledSurface.H" -#include "triSurface.H" +#include "isoSurface.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -52,20 +52,27 @@ namespace Foam class sampledIsoSurface : - public sampledSurface, - public triSurface + public sampledSurface { // Private data //- Field to get isoSurface of - word isoField_; + const word isoField_; //- iso value - scalar isoVal_; + const scalar isoVal_; + + //- Whether to coarse + const Switch regularise_; + + //- Whether to recalculate cell values as average of point values + const Switch average_; //- zone name (if restricted to zones) word zoneName_; + mutable autoPtr<isoSurface> surfPtr_; + //- triangles converted to faceList mutable autoPtr<faceList> facesPtr_; @@ -75,15 +82,21 @@ class sampledIsoSurface //- Time at last call mutable label storedTimeIndex_; - //- For every triangle the original cell in mesh - mutable labelList meshCells_; + //- Cached volfield + mutable autoPtr<volScalarField> storedVolFieldPtr_; + mutable const volScalarField* volFieldPtr_; + + //- Cached pointfield + mutable autoPtr<pointScalarField> storedPointFieldPtr_; + mutable const pointScalarField* pointFieldPtr_; - //- For every unmerged triangle point the point in the triSurface - mutable labelList triPointMergeMap_; // Private Member Functions + //- Get fields needed to recreate iso surface. + void getIsoFields() const; + //- Create iso surface (if time has changed) void createGeometry() const; @@ -127,7 +140,7 @@ public: //- Points of surface virtual const pointField& points() const { - return triSurface::points(); + return surface().points(); } //- Faces of surface @@ -135,7 +148,7 @@ public: { if (!facesPtr_.valid()) { - const triSurface& s = *this; + const triSurface& s = surface(); facesPtr_.reset(new faceList(s.size())); @@ -147,11 +160,19 @@ public: return facesPtr_; } - //- Correct for mesh movement and/or field changes virtual void correct(const bool meshChanged); + const isoSurface& surface() const + { + return surfPtr_(); + } + + //- Lookup or read isoField. Sets volFieldPtr_ and pointFieldPtr_. + void getIsoField(); + + //- sample field on surface virtual tmp<scalarField> sample ( diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C new file mode 100644 index 0000000000000000000000000000000000000000..c258761ab7fabc03def302d94cb1f24ceb606313 --- /dev/null +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C @@ -0,0 +1,344 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "sampledIsoSurfaceCell.H" +#include "dictionary.H" +#include "volFields.H" +#include "volPointInterpolation.H" +#include "addToRunTimeSelectionTable.H" +#include "fvMesh.H" +#include "isoSurfaceCell.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(sampledIsoSurfaceCell, 0); + addNamedToRunTimeSelectionTable(sampledSurface, sampledIsoSurfaceCell, word, isoSurfaceCell); +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::sampledIsoSurfaceCell::createGeometry() const +{ + const fvMesh& fvm = static_cast<const fvMesh&>(mesh()); + + if (fvm.time().timeIndex() != storedTimeIndex_) + { + storedTimeIndex_ = fvm.time().timeIndex(); + + // Clear any stored topo + facesPtr_.clear(); + + // Optionally read volScalarField + autoPtr<volScalarField> readFieldPtr_; + + // 1. see if field in database + // 2. see if field can be read + const volScalarField* cellFldPtr = NULL; + if (fvm.foundObject<volScalarField>(isoField_)) + { + if (debug) + { + Info<< "sampledIsoSurfaceCell::createGeometry() : lookup " + << isoField_ << endl; + } + + cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_); + } + else + { + // Bit of a hack. Read field and store. + + if (debug) + { + Info<< "sampledIsoSurfaceCell::createGeometry() : reading " + << isoField_ << " from time " <<fvm.time().timeName() + << endl; + } + + readFieldPtr_.reset + ( + new volScalarField + ( + IOobject + ( + isoField_, + fvm.time().timeName(), + fvm, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ), + fvm + ) + ); + + cellFldPtr = readFieldPtr_.operator->(); + } + const volScalarField& cellFld = *cellFldPtr; + + tmp<pointScalarField> pointFld + ( + volPointInterpolation::New(fvm).interpolate(cellFld) + ); + + if (average_) + { + //- From point field and interpolated cell. + scalarField cellAvg(fvm.nCells(), scalar(0.0)); + labelField nPointCells(fvm.nCells(), 0); + { + for (label pointI = 0; pointI < fvm.nPoints(); pointI++) + { + const labelList& pCells = fvm.pointCells(pointI); + + forAll(pCells, i) + { + label cellI = pCells[i]; + + cellAvg[cellI] += pointFld().internalField()[pointI]; + nPointCells[cellI]++; + } + } + } + forAll(cellAvg, cellI) + { + cellAvg[cellI] /= nPointCells[cellI]; + } + + const isoSurfaceCell iso + ( + fvm, + cellAvg, + pointFld().internalField(), + isoVal_, + regularise_ + ); + + const_cast<sampledIsoSurfaceCell&> + ( + *this + ).triSurface::operator=(iso); + meshCells_ = iso.meshCells(); + } + else + { + //- Direct from cell field and point field. Gives bad continuity. + const isoSurfaceCell iso + ( + fvm, + cellFld.internalField(), + pointFld().internalField(), + isoVal_, + regularise_ + ); + + const_cast<sampledIsoSurfaceCell&> + ( + *this + ).triSurface::operator=(iso); + meshCells_ = iso.meshCells(); + } + + + if (debug) + { + Pout<< "sampledIsoSurfaceCell::createGeometry() : constructed iso:" + << nl + << " regularise : " << regularise_ << nl + << " average : " << average_ << nl + << " isoField : " << isoField_ << nl + << " isoValue : " << isoVal_ << nl + << " points : " << points().size() << nl + << " tris : " << triSurface::size() << nl + << " cut cells : " << meshCells_.size() << endl; + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell +( + const word& name, + const polyMesh& mesh, + const dictionary& dict +) +: + sampledSurface(name, mesh, dict), + isoField_(dict.lookup("isoField")), + isoVal_(readScalar(dict.lookup("isoValue"))), + regularise_(dict.lookupOrDefault("regularise", true)), + average_(dict.lookupOrDefault("average", true)), + zoneName_(word::null), + facesPtr_(NULL), + storedTimeIndex_(-1), + meshCells_(0) +{ +// label zoneId = -1; +// if (dict.readIfPresent("zone", zoneName_)) +// { +// zoneId = mesh.cellZones().findZoneID(zoneName_); +// if (debug && zoneId < 0) +// { +// Info<< "cellZone \"" << zoneName_ +// << "\" not found - using entire mesh" +// << endl; +// } +// } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sampledIsoSurfaceCell::~sampledIsoSurfaceCell() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::sampledIsoSurfaceCell::correct(const bool meshChanged) +{ + // Only change of mesh changes plane - zone restriction gets lost + if (meshChanged) + { + facesPtr_.clear(); + } +} + + +Foam::tmp<Foam::scalarField> +Foam::sampledIsoSurfaceCell::sample +( + const volScalarField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::vectorField> +Foam::sampledIsoSurfaceCell::sample +( + const volVectorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::sphericalTensorField> +Foam::sampledIsoSurfaceCell::sample +( + const volSphericalTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::symmTensorField> +Foam::sampledIsoSurfaceCell::sample +( + const volSymmTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::tensorField> +Foam::sampledIsoSurfaceCell::sample +( + const volTensorField& vField +) const +{ + return sampleField(vField); +} + + +Foam::tmp<Foam::scalarField> +Foam::sampledIsoSurfaceCell::interpolate +( + const interpolation<scalar>& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp<Foam::vectorField> +Foam::sampledIsoSurfaceCell::interpolate +( + const interpolation<vector>& interpolator +) const +{ + return interpolateField(interpolator); +} + +Foam::tmp<Foam::sphericalTensorField> +Foam::sampledIsoSurfaceCell::interpolate +( + const interpolation<sphericalTensor>& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp<Foam::symmTensorField> +Foam::sampledIsoSurfaceCell::interpolate +( + const interpolation<symmTensor>& interpolator +) const +{ + return interpolateField(interpolator); +} + + +Foam::tmp<Foam::tensorField> +Foam::sampledIsoSurfaceCell::interpolate +( + const interpolation<tensor>& interpolator +) const +{ + return interpolateField(interpolator); +} + + +void Foam::sampledIsoSurfaceCell::print(Ostream& os) const +{ + os << "sampledIsoSurfaceCell: " << name() << " :" + << " field:" << isoField_ + << " value:" << isoVal_; + //<< " faces:" << faces().size() // possibly no geom yet + //<< " points:" << points().size(); +} + + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H new file mode 100644 index 0000000000000000000000000000000000000000..7d6868ae98e8515302b6666292cab0e362893a96 --- /dev/null +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H @@ -0,0 +1,238 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sampledIsoSurfaceCell + +Description + A sampledSurface defined by a surface of iso value. Always triangulated. + To be used in sampleSurfaces / functionObjects. Recalculates iso surface + only if time changes. + +SourceFiles + sampledIsoSurfaceCell.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sampledIsoSurfaceCell_H +#define sampledIsoSurfaceCell_H + +#include "sampledSurface.H" +#include "triSurface.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class sampledIsoSurfaceCell Declaration +\*---------------------------------------------------------------------------*/ + +class sampledIsoSurfaceCell +: + public sampledSurface, + public triSurface +{ + // Private data + + //- Field to get isoSurface of + const word isoField_; + + //- iso value + const scalar isoVal_; + + //- Whether to coarse + const Switch regularise_; + + //- Whether to recalculate cell values as average of point values + const Switch average_; + + //- zone name (if restricted to zones) + word zoneName_; + + //- triangles converted to faceList + mutable autoPtr<faceList> facesPtr_; + + + // Recreated for every isoSurface + + //- Time at last call + mutable label storedTimeIndex_; + + //- For every triangle the original cell in mesh + mutable labelList meshCells_; + + + // Private Member Functions + + //- Create iso surface (if time has changed) + void createGeometry() const; + + //- sample field on faces + template <class Type> + tmp<Field<Type> > sampleField + ( + const GeometricField<Type, fvPatchField, volMesh>& vField + ) const; + + + template <class Type> + tmp<Field<Type> > + interpolateField(const interpolation<Type>&) const; + + +public: + + //- Runtime type information + TypeName("sampledIsoSurfaceCell"); + + + // Constructors + + //- Construct from dictionary + sampledIsoSurfaceCell + ( + const word& name, + const polyMesh& mesh, + const dictionary& dict + ); + + + // Destructor + + virtual ~sampledIsoSurfaceCell(); + + + // Member Functions + + //- Points of surface + virtual const pointField& points() const + { + return triSurface::points(); + } + + //- Faces of surface + virtual const faceList& faces() const + { + if (!facesPtr_.valid()) + { + const triSurface& s = *this; + + facesPtr_.reset(new faceList(s.size())); + + forAll(s, i) + { + facesPtr_()[i] = s[i].triFaceFace(); + } + } + return facesPtr_; + } + + + //- Correct for mesh movement and/or field changes + virtual void correct(const bool meshChanged); + + + //- sample field on surface + virtual tmp<scalarField> sample + ( + const volScalarField& + ) const; + + //- sample field on surface + virtual tmp<vectorField> sample + ( + const volVectorField& + ) const; + + //- sample field on surface + virtual tmp<sphericalTensorField> sample + ( + const volSphericalTensorField& + ) const; + + //- sample field on surface + virtual tmp<symmTensorField> sample + ( + const volSymmTensorField& + ) const; + + //- sample field on surface + virtual tmp<tensorField> sample + ( + const volTensorField& + ) const; + + + //- interpolate field on surface + virtual tmp<scalarField> interpolate + ( + const interpolation<scalar>& + ) const; + + //- interpolate field on surface + virtual tmp<vectorField> interpolate + ( + const interpolation<vector>& + ) const; + + //- interpolate field on surface + virtual tmp<sphericalTensorField> interpolate + ( + const interpolation<sphericalTensor>& + ) const; + + //- interpolate field on surface + virtual tmp<symmTensorField> interpolate + ( + const interpolation<symmTensor>& + ) const; + + //- interpolate field on surface + virtual tmp<tensorField> interpolate + ( + const interpolation<tensor>& + ) const; + + //- Write + virtual void print(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "sampledIsoSurfaceCellTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..144a0ddf41b8f4cc4e8ec8945c8b0bed26894ced --- /dev/null +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCellTemplates.C @@ -0,0 +1,89 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "sampledIsoSurfaceCell.H" +#include "isoSurface.H" +#include "volFieldsFwd.H" +#include "pointFields.H" +#include "volPointInterpolation.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +template <class Type> +Foam::tmp<Foam::Field<Type> > +Foam::sampledIsoSurfaceCell::sampleField +( + const GeometricField<Type, fvPatchField, volMesh>& vField +) const +{ + // Recreate geometry if time has changed + createGeometry(); + + return tmp<Field<Type> >(new Field<Type>(vField, meshCells_)); +} + + +template <class Type> +Foam::tmp<Foam::Field<Type> > +Foam::sampledIsoSurfaceCell::interpolateField +( + const interpolation<Type>& interpolator +) const +{ + // Recreate geometry if time has changed + createGeometry(); + + // One value per point + tmp<Field<Type> > tvalues(new Field<Type>(points().size())); + Field<Type>& values = tvalues(); + + boolList pointDone(points().size(), false); + + forAll(faces(), cutFaceI) + { + const face& f = faces()[cutFaceI]; + + forAll(f, faceVertI) + { + label pointI = f[faceVertI]; + + if (!pointDone[pointI]) + { + values[pointI] = interpolator.interpolate + ( + points()[pointI], + meshCells_[cutFaceI] + ); + pointDone[pointI] = true; + } + } + } + + return tvalues; +} + + +// ************************************************************************* // diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C index dc2abbd74b7aaa50d33edec4c0a0ef9141820280..4d4f75c739c2ab441d2ee3c124a69435e1fea991 100644 --- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C @@ -25,7 +25,6 @@ License \*---------------------------------------------------------------------------*/ #include "sampledIsoSurface.H" -#include "isoSurface.H" #include "volFieldsFwd.H" #include "pointFields.H" #include "volPointInterpolation.H" @@ -41,8 +40,7 @@ Foam::sampledIsoSurface::sampleField { // Recreate geometry if time has changed createGeometry(); - - return tmp<Field<Type> >(new Field<Type>(vField, meshCells_)); + return tmp<Field<Type> >(new Field<Type>(vField, surface().meshCells())); } @@ -53,36 +51,28 @@ Foam::sampledIsoSurface::interpolateField const interpolation<Type>& interpolator ) const { - // Recreate geometry if time has changed - createGeometry(); - - // One value per point - tmp<Field<Type> > tvalues(new Field<Type>(points().size())); - Field<Type>& values = tvalues(); + const fvMesh& fvm = static_cast<const fvMesh&>(mesh()); - boolList pointDone(points().size(), false); + // Get fields to sample. Assume volPointInterpolation! + const GeometricField<Type, fvPatchField, volMesh>& volFld = + interpolator.psi(); - forAll(faces(), cutFaceI) - { - const face& f = faces()[cutFaceI]; + tmp<GeometricField<Type, pointPatchField, pointMesh> > pointFld + ( + volPointInterpolation::New(fvm).interpolate(volFld) + ); - forAll(f, faceVertI) - { - label pointI = f[faceVertI]; - - if (!pointDone[pointI]) - { - values[pointI] = interpolator.interpolate - ( - points()[pointI], - meshCells_[cutFaceI] - ); - pointDone[pointI] = true; - } - } - } + // Recreate geometry if time has changed + createGeometry(); - return tvalues; + // Sample. + return surface().interpolate + ( + *volFieldPtr_, + *pointFieldPtr_, + volFld, + pointFld() + ); }