diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C index 859a018ef7324fd5a3a11960d481ffb704376ed0..299e045c9254e6f7d128b59807c56467b6773c6f 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -30,6 +30,8 @@ License #include "mapPolyMesh.H" #include "OFstream.H" #include "EdgeMap.H" +#include "syncTools.H" +#include "triPointRef.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -46,6 +48,7 @@ Foam::tetDecomposer::decompositionTypeNames { { decompositionType::FACE_CENTRE_TRIS, "faceCentre" }, { decompositionType::FACE_DIAG_TRIS, "faceDiagonal" }, + { decompositionType::PYRAMID, "pyramid" }, }; @@ -183,46 +186,93 @@ Foam::tetDecomposer::tetDecomposer(const polyMesh& mesh) void Foam::tetDecomposer::setRefinement ( const decompositionType decomposeType, + const PackedBoolList& decomposeCell, polyTopoChange& meshMod ) { - cellToPoint_.setSize(mesh_.nCells()); + cellToPoint_.setSize(mesh_.nCells(), -1); forAll(mesh_.cellCentres(), celli) { - // Any point on the cell - label masterPointi = mesh_.faces()[mesh_.cells()[celli][0]][0]; + if (decomposeCell[celli]) + { + // Any point on the cell + label masterPointi = mesh_.faces()[mesh_.cells()[celli][0]][0]; + + cellToPoint_[celli] = meshMod.addPoint + ( + mesh_.cellCentres()[celli], + masterPointi, + -1, + true + ); + } + } + - cellToPoint_[celli] = meshMod.addPoint + // Determine for every face whether it borders a cell that is decomposed + PackedBoolList decomposeFace(mesh_.nFaces()); + { + for (label facei = 0; facei < mesh_.nInternalFaces(); facei++) + { + label own = mesh_.faceOwner()[facei]; + label nei = mesh_.faceNeighbour()[facei]; + if (decomposeCell[own] || decomposeCell[nei]) + { + decomposeFace[facei] = true; + } + } + + boolList neiDecomposeCell(mesh_.nFaces()-mesh_.nInternalFaces()); + forAll(neiDecomposeCell, bFacei) + { + label facei = mesh_.nInternalFaces()+bFacei; + label own = mesh_.faceOwner()[facei]; + neiDecomposeCell[bFacei] = decomposeCell[own]; + } + syncTools::swapBoundaryFaceList(mesh_, neiDecomposeCell); + + for ( - mesh_.cellCentres()[celli], - masterPointi, - -1, - true - ); + label facei = mesh_.nInternalFaces(); + facei < mesh_.nFaces(); + facei++ + ) + { + label own = mesh_.faceOwner()[facei]; + label bFacei = facei-mesh_.nInternalFaces(); + if (decomposeCell[own] || neiDecomposeCell[bFacei]) + { + decomposeFace[facei] = true; + } + } } // Add face centre points if (decomposeType == FACE_CENTRE_TRIS) { - faceToPoint_.setSize(mesh_.nFaces()); + faceToPoint_.setSize(mesh_.nFaces(), -1); forAll(mesh_.faceCentres(), facei) { - // Any point on the face - const label masterPointi = mesh_.faces()[facei][0]; - - faceToPoint_[facei] = meshMod.addPoint - ( - mesh_.faceCentres()[facei], - masterPointi, - -1, - true - ); + if (decomposeFace[facei]) + { + // Any point on the face + const label masterPointi = mesh_.faces()[facei][0]; + + faceToPoint_[facei] = meshMod.addPoint + ( + mesh_.faceCentres()[facei], + masterPointi, + -1, + true + ); + } } } - // Per face, per point (faceCentre) or triangle (faceDiag) the added cell + // Per face, per point (faceCentre) or triangle (faceDiag) the (existing + // or added) cell on either side faceOwnerCells_.setSize(mesh_.nFaces()); faceNeighbourCells_.setSize(mesh_.nFaces()); @@ -235,7 +285,7 @@ void Foam::tetDecomposer::setRefinement faceNeighbourCells_[facei].setSize(f.size(), -1); } } - else + else if (decomposeType == FACE_DIAG_TRIS) { // Force construction of diagonal decomposition (void)mesh_.tetBasePtIs(); @@ -247,13 +297,24 @@ void Foam::tetDecomposer::setRefinement faceNeighbourCells_[facei].setSize(f.size()-2, -1); } } + else + { + forAll(faceOwnerCells_, facei) + { + faceOwnerCells_[facei].setSize(1, -1); + faceNeighbourCells_[facei].setSize(1, -1); + } + } + // Add internal cells. Note: done in same order as pyramid triangle + // creation later to maintain same ordering. forAll(mesh_.cells(), celli) { const cell& cFaces = mesh_.cells()[celli]; - EdgeMap<label> edgeToFace(8*cFaces.size()); + // Whether cell has already been modified (all other cells get added) + bool modifiedCell = false; forAll(cFaces, cFacei) { @@ -268,40 +329,73 @@ void Foam::tetDecomposer::setRefinement : faceNeighbourCells_[facei] ); - if (decomposeType == FACE_CENTRE_TRIS) + if (decomposeCell[celli]) { - forAll(f, fp) + if (decomposeType == FACE_CENTRE_TRIS) { - if (cFacei == 0 && fp == 0) + forAll(f, fp) { - // Reuse cell itself - added[fp] = celli; + if (!modifiedCell) + { + // Reuse cell itself + added[fp] = celli; + modifiedCell = true; + } + else + { + added[fp] = meshMod.addCell + ( + -1, // masterPoint + -1, // masterEdge + -1, // masterFace + celli, // masterCell + mesh_.cellZones().whichZone(celli) + ); + } } - else + } + else if (decomposeType == FACE_DIAG_TRIS) + { + for (label triI = 0; triI < f.size()-2; triI++) { - added[fp] = meshMod.addCell - ( - -1, // masterPoint - -1, // masterEdge - -1, // masterFace - celli, // masterCell - mesh_.cellZones().whichZone(celli) - ); + if (!modifiedCell) + { + // Reuse cell itself + added[triI] = celli; + modifiedCell = true; + } + else + { + added[triI] = meshMod.addCell + ( + -1, // masterPoint + -1, // masterEdge + -1, // masterFace + celli, // masterCell + mesh_.cellZones().whichZone(celli) + ); + //Pout<< "For cell:" << celli + // << " at:" << mesh_.cellCentres()[celli] + // << " face:" << facei + // << " at:" << mesh_.faceCentres()[facei] + // << " tri:" << triI + // << " added cell:" << added[triI] << endl; + } } } - } - else - { - for (label triI = 0; triI < f.size()-2; triI++) + else // if (decomposeType == PYRAMID) { - if (cFacei == 0 && triI == 0) + // Pyramidal decomposition. + // Assign same cell to all slots + if (!modifiedCell) { // Reuse cell itself - added[triI] = celli; + added = celli; + modifiedCell = true; } else { - added[triI] = meshMod.addCell + added = meshMod.addCell ( -1, // masterPoint -1, // masterEdge @@ -312,6 +406,11 @@ void Foam::tetDecomposer::setRefinement } } } + else + { + // All vertices/triangles address to original cell + added = celli; + } } } @@ -341,8 +440,10 @@ void Foam::tetDecomposer::setRefinement zoneFlip = fz.flipMap()[fz.whichFace(facei)]; } + //Pout<< "Face:" << facei << " at:" << mesh_.faceCentres()[facei] + // << endl; - if (decomposeType == FACE_CENTRE_TRIS) + if (decomposeType == FACE_CENTRE_TRIS && decomposeFace[facei]) { forAll(f, fp) { @@ -353,6 +454,13 @@ void Foam::tetDecomposer::setRefinement triangle[1] = f[f.fcIndex(fp)]; triangle[2] = faceToPoint_[facei]; + //Pout<< " triangle:" << triangle + // << " points:" + // << UIndirectList<point>(meshMod.points(), triangle) + // << " between:" << addedOwn[fp] + // << " and:" << addedNei[fp] << endl; + + if (fp == 0) { modifyFace @@ -387,6 +495,7 @@ void Foam::tetDecomposer::setRefinement // 2. Within owner cell - to cell centre + if (decomposeCell[own]) { label newOwn = addedOwn[f.rcIndex(fp)]; label newNei = addedOwn[fp]; @@ -410,7 +519,11 @@ void Foam::tetDecomposer::setRefinement ); } // 2b. Within neighbour cell - to cell centre - if (facei < mesh_.nInternalFaces()) + if + ( + facei < mesh_.nInternalFaces() + && decomposeCell[mesh_.faceNeighbour()[facei]] + ) { label newOwn = addedNei[f.rcIndex(fp)]; label newNei = addedNei[fp]; @@ -435,7 +548,7 @@ void Foam::tetDecomposer::setRefinement } } } - else + else if (decomposeType == FACE_DIAG_TRIS && decomposeFace[facei]) { label fp0 = max(mesh_.tetBasePtIs()[facei], 0); label fp = f.fcIndex(fp0); @@ -450,7 +563,7 @@ void Foam::tetDecomposer::setRefinement label nextFp = f.fcIndex(fp); - // Triangle triI consisting of f[fp0], f[fp], f[nextFp] + // Triangle triI consisiting of f[fp0], f[fp], f[nextFp] // 1. Front triangle (decomposition of face itself) @@ -496,29 +609,36 @@ void Foam::tetDecomposer::setRefinement // 2. Within owner cell - diagonal to cell centre if (triI < f.size()-3) { - label newOwn = addedOwn[triI]; - label newNei = addedOwn[nextTri]; + if (decomposeCell[own]) + { + label newOwn = addedOwn[triI]; + label newNei = addedOwn[nextTri]; - triangle[0] = f[fp0]; - triangle[1] = f[nextFp]; - triangle[2] = cellToPoint_[own]; + triangle[0] = f[fp0]; + triangle[1] = f[nextFp]; + triangle[2] = cellToPoint_[own]; - addFace - ( - meshMod, - triangle, - newOwn, - newNei, - f[fp], //point - -1, //edge - -1, //face - -1, //patchi - zoneI, - zoneFlip - ); + addFace + ( + meshMod, + triangle, + newOwn, + newNei, + f[fp], //point + -1, //edge + -1, //face + -1, //patchi + zoneI, + zoneFlip + ); + } // 2b. Within neighbour cell - to cell centre - if (facei < mesh_.nInternalFaces()) + if + ( + facei < mesh_.nInternalFaces() + && decomposeCell[mesh_.faceNeighbour()[facei]] + ) { label newOwn = addedNei[triI]; label newNei = addedNei[nextTri]; @@ -548,6 +668,21 @@ void Foam::tetDecomposer::setRefinement fp = nextFp; } } + else + { + // No decomposition. Use zero'th slot. + modifyFace + ( + meshMod, + f, + facei, + addedOwn[0], + addedNei[0], + patchi, + zoneI, + zoneFlip + ); + } } @@ -565,129 +700,153 @@ void Foam::tetDecomposer::setRefinement { label facei = cFaces[cFacei]; - label zoneI = mesh_.faceZones().whichZone(facei); - bool zoneFlip = false; - if (zoneI != -1) + if (decomposeCell[celli]) { - const faceZone& fz = mesh_.faceZones()[zoneI]; - zoneFlip = fz.flipMap()[fz.whichFace(facei)]; - } - - const face& f = mesh_.faces()[facei]; - //const labelList& fEdges = mesh_.faceEdges()[facei]; - forAll(f, fp) - { - label p0 = f[fp]; - label p1 = f[f.fcIndex(fp)]; - const edge e(p0, p1); - - EdgeMap<label>::const_iterator edgeFnd = edgeToFace.find(e); - if (edgeFnd == edgeToFace.end()) + label zoneI = mesh_.faceZones().whichZone(facei); + bool zoneFlip = false; + if (zoneI != -1) { - edgeToFace.insert(e, facei); + const faceZone& fz = mesh_.faceZones()[zoneI]; + zoneFlip = fz.flipMap()[fz.whichFace(facei)]; } - else - { - // Found the other face on the edge. - label otherFacei = edgeFnd(); - const face& otherF = mesh_.faces()[otherFacei]; - - // Found the other face on the edge. Note that since - // we are looping in the same order the tets added for - // otherFacei will be before those of facei - - label otherFp = otherF.find(p0); - if (otherF.nextLabel(otherFp) == p1) - { - // ok. otherFp is first vertex of edge. - } - else if (otherF.prevLabel(otherFp) == p1) - { - otherFp = otherF.rcIndex(otherFp); - } - else - { - FatalErrorInFunction - << "problem." << abort(FatalError); - } + const face& f = mesh_.faces()[facei]; + //const labelList& fEdges = mesh_.faceEdges()[facei]; + forAll(f, fp) + { + label p0 = f[fp]; + label p1 = f[f.fcIndex(fp)]; + const edge e(p0, p1); - // Triangle from edge to cell centre - if (mesh_.faceOwner()[facei] == celli) + EdgeMap<label>::const_iterator edgeFnd = edgeToFace.find(e); + if (edgeFnd == edgeToFace.end()) { - triangle[0] = p0; - triangle[1] = p1; - triangle[2] = cellToPoint_[celli]; + edgeToFace.insert(e, facei); } else { - triangle[0] = p1; - triangle[1] = p0; - triangle[2] = cellToPoint_[celli]; - } + // Found the other face on the edge. + label otherFacei = edgeFnd(); + const face& otherF = mesh_.faces()[otherFacei]; - // Determine tets on either side - label thisTet, otherTet; + // Found the other face on the edge. Note that since + // we are looping in the same order the tets added for + // otherFacei will be before those of facei - if (decomposeType == FACE_CENTRE_TRIS) - { - if (mesh_.faceOwner()[facei] == celli) - { - thisTet = faceOwnerCells_[facei][fp]; - } - else + label otherFp = otherF.find(p0); + if (otherF.nextLabel(otherFp) == p1) { - thisTet = faceNeighbourCells_[facei][fp]; + // ok. otherFp is first vertex of edge. } - - if (mesh_.faceOwner()[otherFacei] == celli) + else if (otherF.prevLabel(otherFp) == p1) { - otherTet = faceOwnerCells_[otherFacei][otherFp]; + otherFp = otherF.rcIndex(otherFp); } else { - otherTet = - faceNeighbourCells_[otherFacei][otherFp]; + FatalErrorInFunction + << "problem." << abort(FatalError); } - } - else - { - label thisTriI = triIndex(facei, fp); + + + // Triangle from edge to cell centre if (mesh_.faceOwner()[facei] == celli) { - thisTet = faceOwnerCells_[facei][thisTriI]; + triangle[0] = p0; + triangle[1] = p1; + triangle[2] = cellToPoint_[celli]; } else { - thisTet = faceNeighbourCells_[facei][thisTriI]; + triangle[0] = p1; + triangle[1] = p0; + triangle[2] = cellToPoint_[celli]; } - label otherTriI = triIndex(otherFacei, otherFp); - if (mesh_.faceOwner()[otherFacei] == celli) + // Determine tets on either side + label thisTet, otherTet; + + if (decomposeType == FACE_CENTRE_TRIS) + { + if (mesh_.faceOwner()[facei] == celli) + { + thisTet = faceOwnerCells_[facei][fp]; + } + else + { + thisTet = faceNeighbourCells_[facei][fp]; + } + + if (mesh_.faceOwner()[otherFacei] == celli) + { + otherTet = + faceOwnerCells_[otherFacei][otherFp]; + } + else + { + otherTet = + faceNeighbourCells_[otherFacei][otherFp]; + } + } + else if (decomposeType == FACE_DIAG_TRIS) { - otherTet = faceOwnerCells_[otherFacei][otherTriI]; + label thisTriI = triIndex(facei, fp); + if (mesh_.faceOwner()[facei] == celli) + { + thisTet = faceOwnerCells_[facei][thisTriI]; + } + else + { + thisTet = faceNeighbourCells_[facei][thisTriI]; + } + + label otherTriI = triIndex(otherFacei, otherFp); + if (mesh_.faceOwner()[otherFacei] == celli) + { + otherTet = + faceOwnerCells_[otherFacei][otherTriI]; + } + else + { + otherTet = + faceNeighbourCells_[otherFacei][otherTriI]; + } } else { - otherTet = - faceNeighbourCells_[otherFacei][otherTriI]; + if (mesh_.faceOwner()[facei] == celli) + { + thisTet = faceOwnerCells_[facei][0]; + } + else + { + thisTet = faceNeighbourCells_[facei][0]; + } + if (mesh_.faceOwner()[otherFacei] == celli) + { + otherTet = faceOwnerCells_[otherFacei][0]; + } + else + { + otherTet = + faceNeighbourCells_[otherFacei][0]; + } } - } - - addFace - ( - meshMod, - triangle, - otherTet, - thisTet, - -1, //masterPoint - -1, //fEdges[fp], //masterEdge - facei, //masterFace - -1, //patchi - zoneI, - zoneFlip - ); + addFace + ( + meshMod, + triangle, + otherTet, + thisTet, + -1, //masterPoint + -1, //fEdges[fp], //masterEdge + facei, //masterFace + -1, //patchi + zoneI, + zoneFlip + ); + } } } } diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H index 7bc338f331b1cbb3548a0f5c9b1382b585dc8415..d981f6e62ac2f52decfae31c8871eb4b49fcc88f 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2012-2016 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -25,7 +25,10 @@ Class Foam::tetDecomposer Description - Decomposes polyMesh into tets. + Decomposes polyMesh into tets or pyramids. + + Cells neighbouring decomposed cells are not decomposed themselves + so will be polyhedral. SourceFiles tetDecomposer.C @@ -36,7 +39,7 @@ SourceFiles #define tetDecomposer_H #include "DynamicList.H" -#include "bitSet.H" +#include "PackedBoolList.H" #include "boolList.H" #include "typeInfo.H" #include "Enum.H" @@ -65,8 +68,8 @@ public: { FACE_CENTRE_TRIS, //- Faces decomposed into triangles // using face-centre - - FACE_DIAG_TRIS //- Faces decomposed into triangles diagonally + FACE_DIAG_TRIS, //- Faces decomposed into triangles diagonally + PYRAMID //- Faces not decomposed (generates pyramids) }; static const Enum<decompositionType> decompositionTypeNames; @@ -155,7 +158,7 @@ public: return cellToPoint_; } - //- From face to tet point + //- From face to tet point (only for faceCentre) const labelList& faceToPoint() const { return faceToPoint_; @@ -163,14 +166,16 @@ public: //- Per face, per point (faceCentre) or triangle (faceDiag) - // the added tet on the owner side + // the added tet on the owner side. For non-face (pyramid) + // size 1. const labelListList& faceOwnerCells() const { return faceOwnerCells_; } //- Per face, per point (faceCentre) or triangle (faceDiag) - // the added tet on the neighbour side + // the added tet on the neighbour side. For non-face (pyramid) + // size 1. const labelListList& faceNeighbourCells() const { return faceNeighbourCells_; @@ -184,6 +189,7 @@ public: void setRefinement ( const decompositionType decomposeType, + const PackedBoolList& decomposeCell, polyTopoChange& meshMod ); diff --git a/src/meshTools/regionSplit/localPointRegion.C b/src/meshTools/regionSplit/localPointRegion.C index c0b05d90810d561779f49e03e574bb6500b5d527..c2ba3e5314d5e39f62091f83324c7bbce6a53638 100644 --- a/src/meshTools/regionSplit/localPointRegion.C +++ b/src/meshTools/regionSplit/localPointRegion.C @@ -574,7 +574,8 @@ Foam::labelList Foam::localPointRegion::findDuplicateFaces { FatalErrorInFunction << "Face:" << bFacei + mesh.nInternalFaces() - << " has local points:" << f + << " has local points:" << f << " at:" + << UIndirectList<point>(allPatch.localPoints(), f) << " which are in same order as face:" << otherFacei + mesh.nInternalFaces() << " with local points:" << otherF @@ -597,7 +598,8 @@ Foam::labelList Foam::localPointRegion::findDuplicateFaces << "This means that three or more faces share" << " the same points and this is illegal." << nl << "Face:" << meshFace0 - << " with local points:" << f + << " with local points:" << f << " at:" + << UIndirectList<point>(allPatch.localPoints(), f) << " which are in same order as face:" << meshFace1 << " with local points:" << otherF diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C index c79da8b2d25e5991469adf66babecae719a8638f..1a067714c2bd2a6b0217fca4210d8ad318bc9320 100644 --- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C +++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2015-2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -867,7 +867,7 @@ void Foam::surfaceIntersection::doCutEdges maskFaces.append(pHit.index()); - if (maskRegions[surf1[pHit.index()].region()]) + if (maskRegions.test(surf1[pHit.index()].region())) { continue; } diff --git a/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.C b/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.C index e71ed39b5dcf283c7a6ae8411e6b3e3097fe7e77..aaddd063a3023c17fba0043640c1862943d8629c 100644 --- a/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.C +++ b/src/parallel/decompose/decompositionMethods/geomDecomp/geomDecomp.C @@ -67,11 +67,12 @@ void Foam::geomDecomp::checkDecompositionDirections { if (n_[dir] > 1 && meshDirs[dir] == -1) { - FatalErrorInFunction << "Trying to decompose a 1/2D mesh" + WarningInFunction + << "Trying to decompose a 1/2D mesh" << " into " << n_[dir] << " parts in direction " << Vector<label>::componentNames[dir] - << exit(FatalError); + << endl; } } }