From a9f25c19dc618e4d8b87bec7e1d22430652095c0 Mon Sep 17 00:00:00 2001 From: mattijs <mattijs> Date: Mon, 14 Feb 2011 14:52:51 +0000 Subject: [PATCH] ENH: extrudeMesh,snappyHexMesh: extrude across non-manifold boundaries --- ReleaseNotes-dev | 8 + .../extrude/extrudeMesh/extrudeMesh.C | 180 +++++++- .../polyTopoChange/addPatchCellLayer.C | 423 +++++++++++++----- .../polyTopoChange/addPatchCellLayer.H | 73 ++- .../autoHexMeshDriver/autoLayerDriver.C | 198 +++++--- .../autoHexMeshDriver/autoLayerDriver.H | 12 +- .../meshRefinement/meshRefinement.C | 108 +++-- .../meshRefinement/meshRefinement.H | 13 +- .../refinementSurfaces/refinementSurfaces.C | 78 +++- .../refinementSurfaces/refinementSurfaces.H | 15 +- 10 files changed, 847 insertions(+), 261 deletions(-) diff --git a/ReleaseNotes-dev b/ReleaseNotes-dev index 227041536f1..230a1168263 100644 --- a/ReleaseNotes-dev +++ b/ReleaseNotes-dev @@ -104,6 +104,11 @@ taking film into account + Parallel aware *** *New* ptscotch decomposition method +*** *Updated* scotch decomposition method to run in parallel by doing + decomposition on the master. Unfortunately scotch and ptscotch cannot + be linked in to the same executable. +*** *Updated* simple decomposition method to run in parallel by doing + decomposition on the master. *** *Updated* decomposePar maps polyPatches instead of recreating them so polyPatches holding data can map the data. *** *Updated* particle tracking algorithm @@ -183,6 +188,9 @@ e.g. pitzDailyDirectMapped tutorial. + =setSet=: allows time range (e.g. 0:100) in combination with -batch argument to execute the commands for multiple times. + + =extrudeMesh=: + - option to add extrusion to existing mesh. + - works in parallel * Post-processing + =paraFoam=, =foamToVTK=: full support for polyhedral cell type in recent Paraview versions. diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C index 8a7296eca00..5804d962c81 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C +++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C @@ -43,6 +43,7 @@ Description #include "addPatchCellLayer.H" #include "fvMesh.H" #include "MeshedSurfaces.H" +#include "globalIndex.H" #include "extrudedMesh.H" #include "extrudeModel.H" @@ -261,6 +262,12 @@ int main(int argc, char *argv[]) sourceCasePath.expand(); fileName sourceRootDir = sourceCasePath.path(); fileName sourceCaseDir = sourceCasePath.name(); + if (Pstream::parRun()) + { + sourceCaseDir = + sourceCaseDir + /"processor" + Foam::name(Pstream::myProcNo()); + } wordList sourcePatches; dict.lookup("sourcePatches") >> sourcePatches; @@ -279,54 +286,173 @@ int main(int argc, char *argv[]) sourceRootDir, sourceCaseDir ); + #include "createMesh.H" const polyBoundaryMesh& patches = mesh.boundaryMesh(); - // Topo change container. Either copy an existing mesh or start - // with empty storage (number of patches only needed for checking) - autoPtr<polyTopoChange> meshMod - ( - ( - mode == MESH - ? new polyTopoChange(mesh) - : new polyTopoChange(patches.size()) - ) - ); - // Extrusion engine. Either adding to existing mesh or // creating separate mesh. addPatchCellLayer layerExtrude(mesh, (mode == MESH)); + const labelList meshFaces(patchFaces(patches, sourcePatches)); indirectPrimitivePatch extrudePatch ( IndirectList<face> ( mesh.faces(), - patchFaces(patches, sourcePatches) + meshFaces ), mesh.points() ); + // Determine extrudePatch normal + pointField extrudePatchPointNormals + ( + PatchTools::pointNormals //calcNormals + ( + mesh, + extrudePatch, + meshFaces + ) + ); + + + // Precalculate mesh edges for pp.edges. + const labelList meshEdges + ( + extrudePatch.meshEdges + ( + mesh.edges(), + mesh.pointEdges() + ) + ); + + // Global face indices engine + const globalIndex globalFaces(mesh.nFaces()); + + // Determine extrudePatch.edgeFaces in global numbering (so across + // coupled patches) + labelListList edgeGlobalFaces + ( + addPatchCellLayer::globalEdgeFaces + ( + mesh, + globalFaces, + extrudePatch + ) + ); + + + // Determine what patches boundary edges need to get extruded into. + // This might actually cause edge-connected processors to become + // face-connected so might need to introduce new processor boundaries. + // Calculates: + // - per pp.edge the patch to extrude into + // - any additional processor boundaries (patchToNbrProc = map from + // new patchID to neighbour processor) + // - number of new patches (nPatches) + + labelList sidePatchID; + label nPatches; + Map<label> nbrProcToPatch; + Map<label> patchToNbrProc; + addPatchCellLayer::calcSidePatch + ( + mesh, + globalFaces, + edgeGlobalFaces, + extrudePatch, + + sidePatchID, + nPatches, + nbrProcToPatch, + patchToNbrProc + ); + + + // Add any patches. + + label nAdded = nPatches - mesh.boundaryMesh().size(); + reduce(nAdded, sumOp<label>()); + + Info<< "Adding overall " << nAdded << " processor patches." << endl; + + if (nAdded > 0) + { + DynamicList<polyPatch*> newPatches(nPatches); + forAll(mesh.boundaryMesh(), patchI) + { + newPatches.append + ( + mesh.boundaryMesh()[patchI].clone + ( + mesh.boundaryMesh() + ).ptr() + ); + } + for + ( + label patchI = mesh.boundaryMesh().size(); + patchI < nPatches; + patchI++ + ) + { + label nbrProcI = patchToNbrProc[patchI]; + word name = + "procBoundary" + + Foam::name(Pstream::myProcNo()) + + "to" + + Foam::name(nbrProcI); + + Pout<< "Adding patch " << patchI + << " name:" << name + << " between " << Pstream::myProcNo() + << " and " << nbrProcI + << endl; + + + newPatches.append + ( + new processorPolyPatch + ( + name, + 0, // size + mesh.nFaces(), // start + patchI, // index + mesh.boundaryMesh(),// polyBoundaryMesh + Pstream::myProcNo(),// myProcNo + nbrProcI // neighbProcNo + ) + ); + } + + // Add patches. Do no parallel updates. + mesh.removeFvBoundary(); + mesh.addFvPatches(newPatches, true); + } + + + // Only used for addPatchCellLayer into new mesh - labelList exposedPatchIDs; + labelList exposedPatchID; if (mode == PATCH) { dict.lookup("exposedPatchName") >> backPatchName; - exposedPatchIDs.setSize + exposedPatchID.setSize ( extrudePatch.size(), findPatchID(patches, backPatchName) ); } - + // Determine points and extrusion pointField layer0Points(extrudePatch.nPoints()); pointField displacement(extrudePatch.nPoints()); forAll(displacement, pointI) { - const vector& patchNormal = extrudePatch.pointNormals()[pointI]; + const vector& patchNormal = extrudePatchPointNormals[pointI]; // layer0 point layer0Points[pointI] = model() @@ -373,15 +499,31 @@ int main(int argc, char *argv[]) // Layers per point labelList nPointLayers(extrudePatch.nPoints(), model().nLayers()); // Displacement for first layer - vectorField firstLayerDisp(displacement*model().sumThickness(1)); + vectorField firstLayerDisp = displacement*model().sumThickness(1); + // Expansion ratio not used. scalarField ratio(extrudePatch.nPoints(), 1.0); + // Topo change container. Either copy an existing mesh or start + // with empty storage (number of patches only needed for checking) + autoPtr<polyTopoChange> meshMod + ( + ( + mode == MESH + ? new polyTopoChange(mesh) + : new polyTopoChange(patches.size()) + ) + ); + layerExtrude.setRefinement ( + globalFaces, + edgeGlobalFaces, + ratio, // expansion ratio extrudePatch, // patch faces to extrude - exposedPatchIDs, // if new mesh: patches for exposed faces + sidePatchID, // if boundary edge: patch to use + exposedPatchID, // if new mesh: patches for exposed faces nFaceLayers, nPointLayers, firstLayerDisp, @@ -401,7 +543,7 @@ int main(int argc, char *argv[]) model() ( extrudePatch.localPoints()[pointI], - extrudePatch.pointNormals()[pointI], + extrudePatchPointNormals[pointI], pPointI+1 // layer ) ); diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C index 0a906aa7bcf..f8f7ab41b00 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C @@ -43,75 +43,6 @@ defineTypeNameAndDebug(Foam::addPatchCellLayer, 0); // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -// Calculate global faces per pp edge. -Foam::labelListList Foam::addPatchCellLayer::calcGlobalEdgeFaces -( - const polyMesh& mesh, - const globalIndex& globalFaces, - const indirectPrimitivePatch& pp, - const labelList& meshEdges -) -{ - //// Determine coupled edges just so we don't have to have storage - //// for all non-coupled edges. - // - //PackedBoolList isCoupledEdge(mesh.nEdges()); - // - //const polyBoundaryMesh& patches = mesh.boundaryMesh(); - // - //forAll(patches, patchI) - //{ - // const polyPatch& pp = patches[patchI]; - // - // if (pp.coupled()) - // { - // const labelList& meshEdges = pp.meshEdges(); - // - // forAll(meshEdges, i) - // { - // isCoupledEdge.set(meshEdges[i], 1); - // } - // } - //} - - // From mesh edge to global face labels. Sized only for pp edges. - labelListList globalEdgeFaces(mesh.nEdges()); - - const labelListList& edgeFaces = pp.edgeFaces(); - - forAll(edgeFaces, edgeI) - { - label meshEdgeI = meshEdges[edgeI]; - - //if (isCoupledEdge.get(meshEdgeI) == 1) - { - const labelList& eFaces = edgeFaces[edgeI]; - - // Store face and processor as unique tag. - labelList& globalEFaces = globalEdgeFaces[meshEdgeI]; - globalEFaces.setSize(eFaces.size()); - forAll(eFaces, i) - { - globalEFaces[i] = - globalFaces.toGlobal(pp.addressing()[eFaces[i]]); - } - } - } - - // Synchronise across coupled edges. - syncTools::syncEdgeList - ( - mesh, - globalEdgeFaces, - uniqueEqOp(), - labelList() // null value - ); - - // Extract pp part - return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges)); -} - - Foam::label Foam::addPatchCellLayer::nbrFace ( const labelListList& edgeFaces, @@ -316,12 +247,12 @@ Foam::labelPair Foam::addPatchCellLayer::getEdgeString Foam::label Foam::addPatchCellLayer::addSideFace ( const indirectPrimitivePatch& pp, - const labelList& patchID, // prestored patch per pp face const labelListList& addedCells, // per pp face the new extruded cell const face& newFace, + const label newPatchID, + const label ownFaceI, // pp face that provides owner const label nbrFaceI, - const label patchEdgeI, // edge to add to const label meshEdgeI, // corresponding mesh edge const label layerI, // layer const label numEdgeFaces, // number of layers for edge @@ -329,8 +260,9 @@ Foam::label Foam::addPatchCellLayer::addSideFace polyTopoChange& meshMod ) const { - // Edge to 'inflate' from + // Face or edge to 'inflate' from label inflateEdgeI = -1; + label inflateFaceI = -1; // Check mesh faces using edge if (addToMesh_) @@ -346,8 +278,6 @@ Foam::label Foam::addPatchCellLayer::addSideFace } } - // Get my mesh face and its zone. - label meshFaceI = pp.addressing()[ownFaceI]; // Zone info comes from any side patch face. Otherwise -1 since we // don't know what to put it in - inherit from the extruded faces? label zoneI = -1; //mesh_.faceZones().whichZone(meshFaceI); @@ -358,14 +288,15 @@ Foam::label Foam::addPatchCellLayer::addSideFace // Is patch edge external edge of indirectPrimitivePatch? if (nbrFaceI == -1) { - // External edge so external face. Patch id is obtained from - // any other patch connected to edge. + // External edge so external face. const polyBoundaryMesh& patches = mesh_.boundaryMesh(); // Loop over all faces connected to edge to inflate and - // see if any boundary face (but not meshFaceI) - label otherPatchID = patchID[ownFaceI]; + // see if we can find a face that is otherPatchID + + // Get my mesh face and its zone. + label meshFaceI = pp.addressing()[ownFaceI]; forAll(meshFaces, k) { @@ -373,11 +304,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace if ( - faceI != meshFaceI - && !mesh_.isInternalFace(faceI) + (faceI != meshFaceI) + && (patches.whichPatch(faceI) == newPatchID) ) { - otherPatchID = patches.whichPatch(faceI); + // Found the patch face. Use it to inflate from + inflateEdgeI = -1; + inflateFaceI = faceI; + zoneI = mesh_.faceZones().whichZone(faceI); if (zoneI != -1) { @@ -414,7 +348,7 @@ Foam::label Foam::addPatchCellLayer::addSideFace //Pout<< "Added boundary face:" << newFace // << " own:" << addedCells[ownFaceI][layerOwn] - // << " patch:" << otherPatchID + // << " patch:" << newPatchID // << endl; addedFaceI = meshMod.setAction @@ -426,9 +360,9 @@ Foam::label Foam::addPatchCellLayer::addSideFace -1, // neighbour -1, // master point inflateEdgeI, // master edge - -1, // master face + inflateFaceI, // master face false, // flux flip - otherPatchID, // patch for face + newPatchID, // patch for face zoneI, // zone for face flip // face zone flip ) @@ -510,6 +444,51 @@ Foam::label Foam::addPatchCellLayer::addSideFace } +Foam::label Foam::addPatchCellLayer::findProcPatch +( + const polyMesh& mesh, + const label nbrProcID +) +{ + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + forAll(mesh.globalData().processorPatches(), i) + { + label patchI = mesh.globalData().processorPatches()[i]; + + if + ( + refCast<const processorPolyPatch>(patches[patchI]).neighbProcNo() + == nbrProcID + ) + { + return patchI; + } + } + return -1; +} + + +void Foam::addPatchCellLayer::setFaceProps +( + const polyMesh& mesh, + const label faceI, + + label& patchI, + label& zoneI, + bool& zoneFlip +) +{ + patchI = mesh.boundaryMesh().whichPatch(faceI); + zoneI = mesh.faceZones().whichZone(faceI); + if (zoneI != -1) + { + label index = mesh.faceZones()[zoneI].whichFace(faceI); + zoneFlip = mesh.faceZones()[zoneI].flipMap()[index]; + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct from mesh @@ -561,10 +540,251 @@ Foam::labelListList Foam::addPatchCellLayer::addedCells() const } +// Calculate global faces per pp edge. +Foam::labelListList Foam::addPatchCellLayer::globalEdgeFaces +( + const polyMesh& mesh, + const globalIndex& globalFaces, + const indirectPrimitivePatch& pp +) +{ + // Precalculate mesh edges for pp.edges. + const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges())); + + // From mesh edge to global face labels. Non-empty sublists only for + // pp edges. + labelListList globalEdgeFaces(mesh.nEdges()); + + const labelListList& edgeFaces = pp.edgeFaces(); + + forAll(edgeFaces, edgeI) + { + label meshEdgeI = meshEdges[edgeI]; + + const labelList& eFaces = edgeFaces[edgeI]; + + // Store face and processor as unique tag. + labelList& globalEFaces = globalEdgeFaces[meshEdgeI]; + globalEFaces.setSize(eFaces.size()); + forAll(eFaces, i) + { + globalEFaces[i] = globalFaces.toGlobal(pp.addressing()[eFaces[i]]); + } + } + + // Synchronise across coupled edges. + syncTools::syncEdgeList + ( + mesh, + globalEdgeFaces, + uniqueEqOp(), + labelList() // null value + ); + + // Extract pp part + return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges)); +} + + +void Foam::addPatchCellLayer::calcSidePatch +( + const polyMesh& mesh, + const globalIndex& globalFaces, + const labelListList& globalEdgeFaces, + const indirectPrimitivePatch& pp, + + labelList& sidePatchID, + label& nPatches, + Map<label>& nbrProcToPatch, + Map<label>& patchToNbrProc +) +{ + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + // Precalculate mesh edges for pp.edges. + const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges())); + + sidePatchID.setSize(pp.nEdges()); + sidePatchID = -1; + + // These also get determined but not (yet) exported: + // - whether face is created from other face or edge + // - what zone&orientation face should have + + labelList inflateEdgeI(pp.nEdges(), -1); + labelList inflateFaceI(pp.nEdges(), -1); + labelList sideZoneID(pp.nEdges(), -1); + boolList sideFlip(pp.nEdges(), false); + + nPatches = patches.size(); + + forAll(globalEdgeFaces, edgeI) + { + const labelList& eGlobalFaces = globalEdgeFaces[edgeI]; + if + ( + eGlobalFaces.size() == 2 + && pp.edgeFaces()[edgeI].size() == 1 + ) + { + // Locally but not globally a boundary edge. Hence a coupled + // edge. Find the patch to use if on different + // processors. + + label f0 = eGlobalFaces[0]; + label f1 = eGlobalFaces[1]; + + label otherProcI = -1; + if (globalFaces.isLocal(f0) && !globalFaces.isLocal(f1)) + { + otherProcI = globalFaces.whichProcID(f1); + } + else if (!globalFaces.isLocal(f0) && globalFaces.isLocal(f1)) + { + otherProcI = globalFaces.whichProcID(f0); + } + + + if (otherProcI != -1) + { + sidePatchID[edgeI] = findProcPatch(mesh, otherProcI); + if (sidePatchID[edgeI] == -1) + { + // Cannot find a patch to processor. See if already + // marked for addition + if (nbrProcToPatch.found(otherProcI)) + { + sidePatchID[edgeI] = nbrProcToPatch[otherProcI]; + } + else + { + sidePatchID[edgeI] = nPatches; + nbrProcToPatch.insert(otherProcI, nPatches); + patchToNbrProc.insert(nPatches, otherProcI); + nPatches++; + } + } + } + } + } + + + + // Determine face properties for all other boundary edges + // ------------------------------------------------------ + + const labelListList& edgeFaces = pp.edgeFaces(); + forAll(edgeFaces, edgeI) + { + if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1) + { + // Proper, uncoupled patch edge. + + label myFaceI = pp.addressing()[edgeFaces[edgeI][0]]; + + // Pick up any boundary face on this edge and use its properties + label meshEdgeI = meshEdges[edgeI]; + const labelList& meshFaces = mesh.edgeFaces()[meshEdgeI]; + + forAll(meshFaces, k) + { + label faceI = meshFaces[k]; + + if (faceI != myFaceI && !mesh.isInternalFace(faceI)) + { + setFaceProps + ( + mesh, + faceI, + + sidePatchID[edgeI], + sideZoneID[edgeI], + sideFlip[edgeI] + ); + inflateFaceI[edgeI] = faceI; + inflateEdgeI[edgeI] = -1; + + break; + } + } + } + } + + + + // Now hopefully every boundary edge has a side patch. Check + forAll(edgeFaces, edgeI) + { + if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1) + { + const edge& e = pp.edges()[edgeI]; + FatalErrorIn("addPatchCellLayer::calcSidePatch(..)") + << "Have no sidePatchID for edge " << edgeI << " points " + << pp.points()[pp.meshPoints()[e[0]]] + << pp.points()[pp.meshPoints()[e[1]]] + << abort(FatalError); + } + } + + + + // Now we have sidepatch see if we have patchface or edge to inflate + // from. + forAll(edgeFaces, edgeI) + { + if (edgeFaces[edgeI].size() == 1 && inflateFaceI[edgeI] == -1) + { + // 1. Do we have a boundary face to inflate from + + label myFaceI = pp.addressing()[edgeFaces[edgeI][0]]; + + // Pick up any boundary face on this edge and use its properties + label meshEdgeI = meshEdges[edgeI]; + const labelList& meshFaces = mesh.edgeFaces()[meshEdgeI]; + + forAll(meshFaces, k) + { + label faceI = meshFaces[k]; + + if (faceI != myFaceI) + { + if (mesh.isInternalFace(faceI)) + { + inflateEdgeI[edgeI] = meshEdgeI; + } + else + { + if (patches.whichPatch(faceI) == sidePatchID[edgeI]) + { + setFaceProps + ( + mesh, + faceI, + + sidePatchID[edgeI], + sideZoneID[edgeI], + sideFlip[edgeI] + ); + inflateFaceI[edgeI] = faceI; + inflateEdgeI[edgeI] = -1; + + break; + } + } + } + } + } + } +} + + void Foam::addPatchCellLayer::setRefinement ( + const globalIndex& globalFaces, + const labelListList& globalEdgeFaces, const scalarField& expansionRatio, const indirectPrimitivePatch& pp, + const labelList& sidePatchID, const labelList& exposedPatchID, const labelList& nFaceLayers, const labelList& nPointLayers, @@ -575,7 +795,7 @@ void Foam::addPatchCellLayer::setRefinement if (debug) { Pout<< "addPatchCellLayer::setRefinement : Adding up to " - << max(nPointLayers) + << gMax(nPointLayers) << " layers of cells to indirectPrimitivePatch with " << pp.nPoints() << " points" << endl; } @@ -788,8 +1008,6 @@ void Foam::addPatchCellLayer::setRefinement label meshEdgeI = meshEdges[edgeI]; - // Mesh faces using edge - // Mesh faces using edge const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef); @@ -1213,22 +1431,6 @@ void Foam::addPatchCellLayer::setRefinement } - // Global indices engine - const globalIndex globalFaces(mesh_.nFaces()); - - // Get for all pp edgeFaces a unique faceID - labelListList globalEdgeFaces - ( - calcGlobalEdgeFaces - ( - mesh_, - globalFaces, - pp, - meshEdges - ) - ); - - // Mark off which edges have been extruded boolList doneEdge(pp.nEdges(), false); @@ -1474,16 +1676,17 @@ void Foam::addPatchCellLayer::setRefinement addSideFace ( pp, - patchID, addedCells, - newFace, + + newFace, // vertices of new face + sidePatchID[startEdgeI],// -1 or patch for face + patchFaceI, nbrFaceI, - startEdgeI, // edge to inflate from - meshEdgeI, // corresponding mesh edge - i, - numEdgeSideFaces, - meshFaces, + meshEdgeI, // (mesh) edge to inflate + i, // layer + numEdgeSideFaces, // num layers + meshFaces, // edgeFaces meshMod ); } diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H index 4835ac68e39..2458de9fa45 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H @@ -181,16 +181,6 @@ class addPatchCellLayer // Private Member Functions - //- Per patch edge the pp faces (in global indices) using it. Uses - // uniqueEqOp() to remove duplicates. - labelListList calcGlobalEdgeFaces - ( - const polyMesh& mesh, - const globalIndex& globalFaces, - const indirectPrimitivePatch& pp, - const labelList& meshEdges - ); - //- Get the face on the other side of the edge. static label nbrFace ( @@ -226,12 +216,13 @@ class addPatchCellLayer label addSideFace ( const indirectPrimitivePatch&, - const labelList& patchID, const labelListList& addedCells, + const face& newFace, + const label newPatchID, + const label ownFaceI, const label nbrFaceI, - const label patchEdgeI, const label meshEdgeI, const label layerI, const label numEdgeFaces, @@ -239,6 +230,18 @@ class addPatchCellLayer polyTopoChange& ) const; + //- Find patch to neighbouring processor + static label findProcPatch(const polyMesh&, const label nbrProcID); + + //- Extract properties from mesh face + static void setFaceProps + ( + const polyMesh&, + const label, + label&, + label&, + bool& + ); //- Disallow default bitwise copy construct addPatchCellLayer(const addPatchCellLayer&); @@ -256,7 +259,7 @@ public: // Constructors //- Construct from mesh. - addPatchCellLayer(const polyMesh& mesh, const bool addToMesh = true); + addPatchCellLayer(const polyMesh&, const bool addToMesh = true); // Member Functions @@ -291,6 +294,33 @@ public: // Edit + //- Per patch edge the pp faces (in global indices) using it. Uses + // uniqueEqOp() to remove duplicates. + static labelListList globalEdgeFaces + ( + const polyMesh&, + const globalIndex& globalFaces, + const indirectPrimitivePatch& pp + ); + + //- Boundary edges get extruded into boundary faces. Determine patch + // for these faces. This might be a to-be-created processor patch + // (patchI >= mesh.boundaryMesh().size()) in which case the + // nbrProcToPatch, patchToNbrProc give the correspondence. nPatches + // is the new number of patches. + static void calcSidePatch + ( + const polyMesh&, + const globalIndex& globalFaces, + const labelListList& globalEdgeFaces, + const indirectPrimitivePatch& pp, + + labelList& sidePatchID, + label& nPatches, + Map<label>& nbrProcToPatch, + Map<label>& patchToNbrProc + ); + //- Play commands into polyTopoChange to create layers on top // of indirectPrimitivePatch (have to be outside faces). // Gets displacement per patch point. @@ -313,8 +343,11 @@ public: // (instead of e.g. from patch faces) void setRefinement ( + const globalIndex& globalFaces, + const labelListList& globalEdgeFaces, const scalarField& expansionRatio, const indirectPrimitivePatch& pp, + const labelList& sidePatchID, const labelList& exposedPatchID, const labelList& nFaceLayers, const labelList& nPointLayers, @@ -326,20 +359,26 @@ public: //- Add with constant expansion ratio and same nLayers everywhere void setRefinement ( + const globalIndex& globalFaces, + const labelListList& globalEdgeFaces, const label nLayers, const indirectPrimitivePatch& pp, + const labelList& sidePatchID, const vectorField& overallDisplacement, polyTopoChange& meshMod ) { setRefinement ( - scalarField(pp.nPoints(), 1.0), // expansion ration + globalFaces, + globalEdgeFaces, + scalarField(pp.nPoints(), 1.0), // expansion ration pp, + sidePatchID, labelList(0), - labelList(pp.size(), nLayers), - labelList(pp.nPoints(), nLayers), - overallDisplacement / nLayers, + labelList(pp.size(), nLayers), // nFaceLayers + labelList(pp.nPoints(), nLayers), // nPointLayers + overallDisplacement / nLayers, // firstLayerDisp meshMod ); } diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index 46d4ae13431..164da09c5ba 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -45,6 +45,7 @@ Description #include "layerParameters.H" #include "combineFaces.H" #include "IOmanip.H" +#include "globalIndex.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -888,68 +889,6 @@ void Foam::autoLayerDriver::handleNonManifolds Info<< "Outside of local patch is multiply connected across edges or" << " points at " << nNonManif << " points." << endl; - - const labelList& meshPoints = pp.meshPoints(); - - - // 2. Parallel check - // For now disable extrusion at any shared edges. - const labelHashSet sharedEdgeSet(mesh.globalData().sharedEdgeLabels()); - - forAll(pp.edges(), edgeI) - { - if (sharedEdgeSet.found(meshEdges[edgeI])) - { - const edge& e = mesh.edges()[meshEdges[edgeI]]; - - Pout<< "Disabling extrusion at edge " - << mesh.points()[e[0]] << mesh.points()[e[1]] - << " since it is non-manifold across coupled patches." - << endl; - - nonManifoldPoints.insert(e[0]); - nonManifoldPoints.insert(e[1]); - } - } - - // 3b. extrusion can produce multiple faces between 2 cells - // across processor boundary - // This occurs when a coupled face shares more than 1 edge with a - // non-coupled boundary face. - // This is now correctly handled by addPatchCellLayer in that it - // extrudes a single face from the stringed up edges. - - - nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>()); - - if (nNonManif > 0) - { - Info<< "Outside of patches is multiply connected across edges or" - << " points at " << nNonManif << " points." << nl - << "Writing " << nNonManif - << " points where this happens to pointSet " - << nonManifoldPoints.name() << nl - << "and setting layer thickness to zero on these points." - << endl; - - nonManifoldPoints.instance() = mesh.time().timeName(); - nonManifoldPoints.write(); - - forAll(meshPoints, patchPointI) - { - if (nonManifoldPoints.found(meshPoints[patchPointI])) - { - unmarkExtrusion - ( - patchPointI, - patchDisp, - patchNLayers, - extrudeStatus - ); - } - } - } - Info<< "Set displacement to zero for all " << nNonManif << " non-manifold points" << endl; } @@ -1348,7 +1287,6 @@ void Foam::autoLayerDriver::setNumLayers } -// Grow no-extrusion layer void Foam::autoLayerDriver::growNoExtrusion ( const indirectPrimitivePatch& pp, @@ -1439,6 +1377,103 @@ void Foam::autoLayerDriver::growNoExtrusion } +void Foam::autoLayerDriver::determineSidePatches +( + const globalIndex& globalFaces, + const labelListList& edgeGlobalFaces, + const indirectPrimitivePatch& pp, + + labelList& sidePatchID +) +{ + // Sometimes edges-to-be-extruded are on more than 2 processors. + // Work out which 2 hold the faces to be extruded and thus which procpatch + // the side-face should be in. As an additional complication this might + // mean that 2 procesors that were only edge-connected now suddenly need + // to become face-connected i.e. have a processor patch between them. + + fvMesh& mesh = meshRefiner_.mesh(); + + // Determine sidePatchID. Any additional processor boundary gets added to + // patchToNbrProc,nbrProcToPatch and nPatches gets set to the new number + // of patches. + label nPatches; + Map<label> nbrProcToPatch; + Map<label> patchToNbrProc; + addPatchCellLayer::calcSidePatch + ( + mesh, + globalFaces, + edgeGlobalFaces, + pp, + + sidePatchID, + nPatches, + nbrProcToPatch, + patchToNbrProc + ); + + label nOldPatches = mesh.boundaryMesh().size(); + label nAdded = returnReduce(nPatches-nOldPatches, sumOp<label>()); + Info<< nl << "Adding in total " << nAdded/2 << " inter-processor patches to" + << " handle extrusion of non-manifold processor boundaries." + << endl; + + if (nAdded > 0) + { + // We might not add patches in same order as in patchToNbrProc + // so prepare to renumber sidePatchID + Map<label> wantedToAddedPatch; + + for (label patchI = nOldPatches; patchI < nPatches; patchI++) + { + label nbrProcI = patchToNbrProc[patchI]; + word name = + "procBoundary" + + Foam::name(Pstream::myProcNo()) + + "to" + + Foam::name(nbrProcI); + + dictionary patchDict; + patchDict.add("type", processorPolyPatch::typeName); + patchDict.add("myProcNo", Pstream::myProcNo()); + patchDict.add("neighbProcNo", nbrProcI); + patchDict.add("nFaces", 0); + patchDict.add("startFace", mesh.nFaces()); + + Pout<< "Adding patch " << patchI + << " name:" << name + << " between " << Pstream::myProcNo() + << " and " << nbrProcI + << endl; + + label procPatchI = meshRefiner_.appendPatch + ( + mesh, + mesh.boundaryMesh().size(), // new patch index + name, + patchDict + ); + wantedToAddedPatch.insert(patchI, procPatchI); + } + + // Renumber sidePatchID + forAll(sidePatchID, i) + { + label patchI = sidePatchID[i]; + Map<label>::const_iterator fnd = wantedToAddedPatch.find(patchI); + if (fnd != wantedToAddedPatch.end()) + { + sidePatchID[i] = fnd(); + } + } + + mesh.clearOut(); + const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()).updateMesh(); + } +} + + void Foam::autoLayerDriver::calculateLayerThickness ( const indirectPrimitivePatch& pp, @@ -2553,6 +2588,37 @@ void Foam::autoLayerDriver::addLayers ) ); + + // Global face indices engine + const globalIndex globalFaces(mesh.nFaces()); + + // Determine extrudePatch.edgeFaces in global numbering (so across + // coupled patches). This is used only to string up edges between coupled + // faces (all edges between same (global)face indices get extruded). + labelListList edgeGlobalFaces + ( + addPatchCellLayer::globalEdgeFaces + ( + mesh, + globalFaces, + pp + ) + ); + + // Determine patches for extruded boundary edges. Calculates if any + // additional processor patches need to be constructed. + + labelList sidePatchID; + determineSidePatches + ( + globalFaces, + edgeGlobalFaces, + pp, + + sidePatchID + ); + + // Construct iterative mesh mover. Info<< "Constructing mesh displacer ..." << endl; @@ -3038,8 +3104,12 @@ void Foam::autoLayerDriver::addLayers // Not add layer if patchDisp is zero. addLayer.setRefinement ( + globalFaces, + edgeGlobalFaces, + invExpansionRatio, pp(), + sidePatchID, // boundary patch for extruded boundary edges labelList(0), // exposed patchIDs, not used for adding layers nPatchFaceLayers, // layers per face nPatchPointLayers, // layers per point diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H index c96e304983e..fd5a1f3aa88 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -243,6 +243,16 @@ class autoLayerDriver List<extrudeMode>& extrudeStatus ) const; + //- See what patches boundaryedges should be extruded into + void determineSidePatches + ( + const globalIndex& globalFaces, + const labelListList& edgeGlobalFaces, + const indirectPrimitivePatch& pp, + + labelList& sidePatchID + ); + //- Calculate pointwise wanted and minimum thickness. // thickness: wanted thickness // minthickness: when to give up and not extrude diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C index 8a749cd60b6..de0ba875879 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C @@ -1549,60 +1549,28 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh) } -// Adds patch if not yet there. Returns patchID. -Foam::label Foam::meshRefinement::addPatch +Foam::label Foam::meshRefinement::appendPatch ( fvMesh& mesh, + const label insertPatchI, const word& patchName, - const dictionary& patchInfo + const dictionary& patchDict ) { - polyBoundaryMesh& polyPatches = - const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); - - const label patchI = polyPatches.findPatchID(patchName); - if (patchI != -1) - { - // Already there - return patchI; - } - - - label insertPatchI = polyPatches.size(); - label startFaceI = mesh.nFaces(); - - forAll(polyPatches, patchI) - { - const polyPatch& pp = polyPatches[patchI]; - - if (isA<processorPolyPatch>(pp)) - { - insertPatchI = patchI; - startFaceI = pp.start(); - break; - } - } - - - // Below is all quite a hack. Feel free to change once there is a better - // mechanism to insert and reorder patches. - // Clear local fields and e.g. polyMesh parallelInfo. mesh.clearOut(); - label sz = polyPatches.size(); - + polyBoundaryMesh& polyPatches = + const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary()); - dictionary patchDict(patchInfo); - patchDict.set("nFaces", 0); - patchDict.set("startFace", startFaceI); + label patchI = polyPatches.size(); // Add polyPatch at the end - polyPatches.setSize(sz+1); + polyPatches.setSize(patchI+1); polyPatches.set ( - sz, + patchI, polyPatch::New ( patchName, @@ -1611,13 +1579,13 @@ Foam::label Foam::meshRefinement::addPatch polyPatches ) ); - fvPatches.setSize(sz+1); + fvPatches.setSize(patchI+1); fvPatches.set ( - sz, + patchI, fvPatch::New ( - polyPatches[sz], // point to newly added polyPatch + polyPatches[patchI], // point to newly added polyPatch mesh.boundary() ) ); @@ -1675,21 +1643,69 @@ Foam::label Foam::meshRefinement::addPatch mesh, calculatedFvPatchField<tensor>::typeName ); + return patchI; +} + + +// Adds patch if not yet there. Returns patchID. +Foam::label Foam::meshRefinement::addPatch +( + fvMesh& mesh, + const word& patchName, + const dictionary& patchInfo +) +{ + polyBoundaryMesh& polyPatches = + const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()); + fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary()); + + const label patchI = polyPatches.findPatchID(patchName); + if (patchI != -1) + { + // Already there + return patchI; + } + + + label insertPatchI = polyPatches.size(); + label startFaceI = mesh.nFaces(); + + forAll(polyPatches, patchI) + { + const polyPatch& pp = polyPatches[patchI]; + + if (isA<processorPolyPatch>(pp)) + { + insertPatchI = patchI; + startFaceI = pp.start(); + break; + } + } + + dictionary patchDict(patchInfo); + patchDict.set("nFaces", 0); + patchDict.set("startFace", startFaceI); + + // Below is all quite a hack. Feel free to change once there is a better + // mechanism to insert and reorder patches. + + label addedPatchI = appendPatch(mesh, insertPatchI, patchName, patchDict); + // Create reordering list // patches before insert position stay as is - labelList oldToNew(sz+1); + labelList oldToNew(addedPatchI+1); for (label i = 0; i < insertPatchI; i++) { oldToNew[i] = i; } // patches after insert position move one up - for (label i = insertPatchI; i < sz; i++) + for (label i = insertPatchI; i < addedPatchI; i++) { oldToNew[i] = i+1; } // appended patch gets moved to insert position - oldToNew[sz] = insertPatchI; + oldToNew[addedPatchI] = insertPatchI; // Shuffle into place polyPatches.reorder(oldToNew); diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H index 7337161bb75..5a1d8ff8211 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -747,8 +747,17 @@ public: // Other topo changes + //- Helper:append patch to end of mesh. + static label appendPatch + ( + fvMesh&, + const label insertPatchI, + const word&, + const dictionary& + ); + //- Helper:add patch to mesh. Update all registered fields. - // Use addMeshedPatch to add patches originating from surfaces. + // Used by addMeshedPatch to add patches originating from surfaces. static label addPatch(fvMesh&, const word& name, const dictionary&); //- Add patch originating from meshing. Update meshedPatches_. diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C index 18e5bf01b46..b9e2d4d3e07 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C +++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -1007,6 +1007,82 @@ void Foam::refinementSurfaces::findNearestRegion } +void Foam::refinementSurfaces::findNearestRegion +( + const labelList& surfacesToTest, + const pointField& samples, + const scalarField& nearestDistSqr, + labelList& hitSurface, + List<pointIndexHit>& hitInfo, + labelList& hitRegion, + vectorField& hitNormal +) const +{ + labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest)); + + // Do the tests. Note that findNearest returns index in geometries. + searchableSurfacesQueries::findNearest + ( + allGeometry_, + geometries, + samples, + nearestDistSqr, + hitSurface, + hitInfo + ); + + // Rework the hitSurface to be surface (i.e. index into surfaces_) + forAll(hitSurface, pointI) + { + if (hitSurface[pointI] != -1) + { + hitSurface[pointI] = surfacesToTest[hitSurface[pointI]]; + } + } + + // Collect the region + hitRegion.setSize(hitSurface.size()); + hitRegion = -1; + hitNormal.setSize(hitSurface.size()); + hitNormal = vector::zero; + + forAll(surfacesToTest, i) + { + label surfI = surfacesToTest[i]; + + // Collect hits for surfI + const labelList localIndices(findIndices(hitSurface, surfI)); + + List<pointIndexHit> localHits + ( + UIndirectList<pointIndexHit> + ( + hitInfo, + localIndices + ) + ); + + // Region + labelList localRegion; + allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion); + + forAll(localIndices, i) + { + hitRegion[localIndices[i]] = localRegion[i]; + } + + // Normal + vectorField localNormal; + allGeometry_[surfaces_[surfI]].getNormal(localHits, localNormal); + + forAll(localIndices, i) + { + hitNormal[localIndices[i]] = localNormal[i]; + } + } +} + + //// Find intersection with max of edge. Return -1 or the surface //// with the highest maxLevel above currentLevel //Foam::label Foam::refinementSurfaces::findHighestIntersection diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H index dc243f316fb..a7728b90824 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H +++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -327,6 +327,19 @@ public: labelList& hitRegion ) const; + //- Find nearest point on surfaces. Return surface, region and + // normal on surface (so not global surface) + void findNearestRegion + ( + const labelList& surfacesToTest, + const pointField& samples, + const scalarField& nearestDistSqr, + labelList& hitSurface, + List<pointIndexHit>& hitInfo, + labelList& hitRegion, + vectorField& hitNormal + ) const; + //- Detect if a point is 'inside' (closed) surfaces. // Returns -1 if not, returns first surface it is. void findInside -- GitLab