diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C index c369d7246ebb6e674018d8d98ed7134bd171211b..b6cb0055d06801a5725ca6dcf3af21943ae1e106 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2015-2022 OpenCFD Ltd. + Copyright (C) 2015-2022,2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -548,61 +548,61 @@ void Foam::addPatchCellLayer::setFaceProps } -void Foam::addPatchCellLayer::findZoneFace -( - const bool useInternalFaces, - const bool useBoundaryFaces, - - const polyMesh& mesh, - const indirectPrimitivePatch& pp, - const label ppEdgeI, - const labelUIndList& excludeFaces, - const labelList& meshFaces, - - label& inflateFaceI, - label& patchI, - label& zoneI, - bool& zoneFlip -) -{ - inflateFaceI = -1; - patchI = -1; - zoneI = -1; - zoneFlip = false; - - forAll(meshFaces, k) - { - label faceI = meshFaces[k]; - - if - ( - !excludeFaces.found(faceI) - && ( - (mesh.isInternalFace(faceI) && useInternalFaces) - || (!mesh.isInternalFace(faceI) && useBoundaryFaces) - ) - ) - { - setFaceProps - ( - mesh, - pp, - ppEdgeI, - faceI, - - patchI, - zoneI, - zoneFlip, - inflateFaceI - ); - - if (zoneI != -1 || patchI != -1) - { - break; - } - } - } -} +//void Foam::addPatchCellLayer::findZoneFace +//( +// const bool useInternalFaces, +// const bool useBoundaryFaces, +// +// const polyMesh& mesh, +// const indirectPrimitivePatch& pp, +// const label ppEdgeI, +// const labelUIndList& excludeFaces, +// const labelList& meshFaces, +// +// label& inflateFaceI, +// label& patchI, +// label& zoneI, +// bool& zoneFlip +//) +//{ +// inflateFaceI = -1; +// patchI = -1; +// zoneI = -1; +// zoneFlip = false; +// +// forAll(meshFaces, k) +// { +// label faceI = meshFaces[k]; +// +// if +// ( +// !excludeFaces.found(faceI) +// && ( +// (mesh.isInternalFace(faceI) && useInternalFaces) +// || (!mesh.isInternalFace(faceI) && useBoundaryFaces) +// ) +// ) +// { +// setFaceProps +// ( +// mesh, +// pp, +// ppEdgeI, +// faceI, +// +// patchI, +// zoneI, +// zoneFlip, +// inflateFaceI +// ); +// +// if (zoneI != -1 || patchI != -1) +// { +// break; +// } +// } +// } +//} // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -611,11 +611,13 @@ void Foam::addPatchCellLayer::findZoneFace Foam::addPatchCellLayer::addPatchCellLayer ( const polyMesh& mesh, - const bool addToMesh + const bool addToMesh, + const bool extrude ) : mesh_(mesh), addToMesh_(addToMesh), + extrude_(extrude), addedPoints_(0), layerFaces_(0) {} @@ -799,6 +801,11 @@ void Foam::addPatchCellLayer::globalEdgeInfo } + bitSet isPpFace(mesh.nFaces()); + isPpFace.set(pp.addressing()); + // Note: no need to sync isPpFace since does not include processor patches + + const faceZoneMesh& fzs = mesh.faceZones(); // Extract zone info into mesh face indexing for ease of addressing @@ -839,6 +846,11 @@ void Foam::addPatchCellLayer::globalEdgeInfo for (const label facei : isInternalOrCoupled) { + if (isPpFace[facei]) + { + continue; + } + const face& f = mesh.faces()[facei]; label prevPointi = f.last(); @@ -891,9 +903,6 @@ void Foam::addPatchCellLayer::globalEdgeInfo const polyBoundaryMesh& patches = mesh.boundaryMesh(); - bitSet isPpFace(mesh.nFaces()); - isPpFace.set(pp.addressing()); - // Note: no need to sync ppFace since does not include processor patches for (const polyPatch& pp : patches) { @@ -1861,15 +1870,18 @@ void Foam::addPatchCellLayer::setRefinement baseF.setSize(f.size(), labelMax); baseF[fp] = pp.meshPoints()[patchPointi]; - if (ppFlip[patchFacei]) + if (extrude_ == ppFlip[patchFacei]) { - // Neighbour stays. Affected points on the owner side. + // either: + // - extrude out of neighbour so new points connect + // to owner + // - or intrude into owner const label celli = mesh_.faceOwner()[meshFacei]; isAffectedCell.set(celli); } else if (mesh_.isInternalFace(meshFacei)) { - // Owner unaffected. Unaffected points on neighbour side + // Owner unaffected. Affected points on neighbour side const label celli = mesh_.faceNeighbour()[meshFacei]; isAffectedCell.set(celli); } @@ -2120,13 +2132,13 @@ void Foam::addPatchCellLayer::setRefinement } else if (!ppFlip[patchFacei]) { - // Normal: extrude from owner face + // Normal: extrude from owner side extrudeCelli = mesh_.faceOwner()[meshFacei]; extrudeZonei = mesh_.cellZones().whichZone(extrudeCelli); } else if (mesh_.isInternalFace(meshFacei)) { - // Extrude from neighbour face (if internal). Might be + // Extrude from neighbour side (if internal). Might be // that it is a coupled face and the other side is // extruded extrudeCelli = mesh_.faceNeighbour()[meshFacei]; @@ -2164,8 +2176,9 @@ void Foam::addPatchCellLayer::setRefinement - // Create faces on top of the original patch faces. - // These faces are created from original patch faces outwards so in order + // Create faces from the original patch faces. + // These faces are created from original patch faces outwards (or inwards) + // so in order // of increasing cell number. So orientation should be same as original // patch face for them to have owner<neighbour. @@ -2214,74 +2227,97 @@ void Foam::addPatchCellLayer::setRefinement // << endl; // Get new neighbour - label own = addedCells[patchFacei][i]; - label nei; - label patchi; - label zoneI = -1; - bool flip = false; + label own; + label nei = -1; + label patchi = -1; + label zonei = -1; + bool zoneFlip = false; bool fluxFlip = false; if (i == addedCells[patchFacei].size()-1) { - // Top layer so is either patch face or connects to - // the other cell - patchi = patchID[patchFacei]; - if (patchi == -1) - { - // Internal face - nei = - ( - !ppFlip[patchFacei] - ? mesh_.faceNeighbour()[meshFacei] - : mesh_.faceOwner()[meshFacei] - ); + // Last layer so + // - extrude : original patch or other cell + // - intrude : original cell - if (ppFlip[patchFacei]) + if (extrude_) + { + patchi = patchID[patchFacei]; + zonei = mesh_.faceZones().whichZone(meshFacei); + if (zonei != -1) { - newFace = newFace.reverseFace(); + const faceZone& fz = mesh_.faceZones()[zonei]; + zoneFlip = fz.flipMap()[fz.whichFace(meshFacei)]; + } + if (patchi == -1) + { + // Internal face between added cell and original + // neighbour / owner + own = + ( + !ppFlip[patchFacei] + ? mesh_.faceNeighbour()[meshFacei] + : mesh_.faceOwner()[meshFacei] + ); + nei = addedCells[patchFacei][i]; + if (ppFlip[patchFacei]) + { + newFace = newFace.reverseFace(); + if (zonei != -1) + { + zoneFlip = !zoneFlip; + } + fluxFlip = true; + } + } + else + { + // Boundary face + own = addedCells[patchFacei][i]; } - //Pout<< "** adding top (internal) face:" - // << " at:" << mesh_.faceCentres()[meshFacei] - // << " own:" << own << " nei:" << nei - // << " patchi:" << patchi - // << " newFace:" << newFace - // << endl; } else { - nei = -1; - } - zoneI = mesh_.faceZones().whichZone(meshFacei); - if (zoneI != -1) - { - const faceZone& fz = mesh_.faceZones()[zoneI]; - flip = fz.flipMap()[fz.whichFace(meshFacei)]; + // Intrude : face connects to original cell + //if (patchID[patchFacei] == -1) + { + // Internal face between added cell and original + // owner / neighbour. TBD: what if ppFlip set on + // boundary face (so with no neighbour) + own = + ( + !ppFlip[patchFacei] + ? mesh_.faceOwner()[meshFacei] + : mesh_.faceNeighbour()[meshFacei] + ); + nei = addedCells[patchFacei][i]; + if (ppFlip[patchFacei]) + { + // Since neighbour now owns the face flip it + newFace = newFace.reverseFace(); + if (zonei != -1) + { + zoneFlip = !zoneFlip; + } + fluxFlip = true; + } + } } } else { // Internal face between layer i and i+1 + own = addedCells[patchFacei][i]; nei = addedCells[patchFacei][i+1]; - patchi = -1; - } - if (nei != -1 && nei < own) - { - // Wrongly oriented internal face - newFace = newFace.reverseFace(); - std::swap(own, nei); - flip = !flip; - fluxFlip = true; - - //Pout<< "Flipped newFace:" - // << newFace.unitNormal(meshMod.points()) - // << " own:" << own - // << " nei:" << nei - // << endl; + if (ppFlip[patchFacei]) + { + // Since neighbour now owns the face flip it + newFace = newFace.reverseFace(); + fluxFlip = true; + } } - - layerFaces_[patchFacei][i+1] = meshMod.setAction ( polyAddFace @@ -2294,8 +2330,8 @@ void Foam::addPatchCellLayer::setRefinement (addToMesh_ ? meshFacei : -1), // master face fluxFlip, // flux flip patchi, // patch for face - zoneI, // zone for face - flip // face zone flip + zonei, // zone for face + zoneFlip // face zone flip ) ); @@ -2307,13 +2343,14 @@ void Foam::addPatchCellLayer::setRefinement // << " n:" << newFace.unitNormal(meshMod.points()) // << " own:" << own << " nei:" << nei // << " patchi:" << patchi + // << " zonei:" << zonei // << endl; } } } // - // Modify owner faces to have addedCells as neighbour + // Modify original face to have addedCells as neighbour // if (addToMesh_) @@ -2322,38 +2359,102 @@ void Foam::addPatchCellLayer::setRefinement { if (addedCells[patchFacei].size()) { - label meshFacei = pp.addressing()[patchFacei]; + const label meshFacei = pp.addressing()[patchFacei]; layerFaces_[patchFacei][0] = meshFacei; const face& f = pp[patchFacei]; - const label own = - ( - !ppFlip[patchFacei] - ? mesh_.faceOwner()[meshFacei] - : mesh_.faceNeighbour()[meshFacei] - ); - const label nei = addedCells[patchFacei][0]; + if (extrude_) + { + const label own = + ( + !ppFlip[patchFacei] + ? mesh_.faceOwner()[meshFacei] + : mesh_.faceNeighbour()[meshFacei] + ); + const label nei = addedCells[patchFacei][0]; - meshMod.setAction - ( - polyModifyFace + meshMod.setAction ( - (ppFlip[patchFacei] ? f.reverseFace() : f),// verts - meshFacei, // label of face - own, // owner - nei, // neighbour - ppFlip[patchFacei], // face flip - -1, // patch for face - true, //false, // remove from zone - -1, //zoneI, // zone for face - false // face flip in zone - ) - ); + polyModifyFace + ( + (ppFlip[patchFacei] ? f.reverseFace() : f),// verts + meshFacei, // label of face + own, // owner + nei, // neighbour + ppFlip[patchFacei], // face flip + -1, // patch for face + true, //false, // remove from zone + -1, //zoneI, // zone for face + false // face flip in zone + ) + ); + } + else + { + label patchi; + label zonei; + bool zoneOrient; + setFaceProps + ( + mesh_, + meshFacei, + + patchi, + zonei, + zoneOrient + ); - //Pout<< "Modified bottom face " << meshFacei + label own; + label nei = -1; + bool flip = false; + if (!ppFlip[patchFacei]) + { + if (patchi == -1) + { + own = mesh_.faceNeighbour()[meshFacei]; + nei = addedCells[patchFacei].first(); + flip = true; + } + else + { + own = addedCells[patchFacei].first(); + } + } + else + { + if (patchi == -1) + { + own = mesh_.faceOwner()[meshFacei]; + nei = addedCells[patchFacei].first(); + } + else + { + own = addedCells[patchFacei].first(); + } + } + + meshMod.setAction + ( + polyModifyFace + ( + (flip ? f.reverseFace() : f), // verts + meshFacei, // label of face + own, // owner + nei, // neighbour + flip, // face flip + patchi, // patch for face + false, // remove from zone + zonei, // zone for face + flip != zoneOrient // face flip in zone + ) + ); + } + + //Pout<< "Modified original face " << meshFacei // << " at:" << mesh_.faceCentres()[meshFacei] - // << " new own:" << own << " new nei:" << nei + // << " new own:" << meshMod.faceOwner()[meshFacei] + // << " new nei:" << meshMod.faceNeighbour()[meshFacei] // << " verts:" << meshMod.faces()[meshFacei] // << " n:" // << meshMod.faces()[meshFacei].unitNormal(meshMod.points()) @@ -2472,9 +2573,10 @@ void Foam::addPatchCellLayer::setRefinement // between the same two faces and extrude string into a single face. forAll(pp, patchFacei) { + // Get edges of face in vertex order const labelList& fEdges = faceEdges[patchFacei]; - forAll(fEdges, fp) + forAll(fEdges, i) { // Get string of edges that needs to be extruded as a single face. // Returned as indices in fEdges. @@ -2695,7 +2797,7 @@ void Foam::addPatchCellLayer::setRefinement // Walked edges as if owner face was extruded. Reverse // for neighbour face extrusion. - if (ppFlip[patchFacei]) + if (extrude_ == ppFlip[patchFacei]) { newFace = newFace.reverseFace(); } @@ -2797,6 +2899,10 @@ void Foam::addPatchCellLayer::setRefinement { face newFace; + // For intrusion correction later on: store displacement of modified + // points so they can be adapted on connected faces. + pointField displacement(mesh_.nPoints(), Zero); + forAll(baseFaces, facei) { const face& f = mesh_.faces()[facei]; @@ -2815,18 +2921,23 @@ void Foam::addPatchCellLayer::setRefinement const label meshPointi = f[fp]; if (baseF[fp] != labelMax) { - // Duplicated point - const label patchPointi = pp.meshPointMap()[meshPointi]; - const label addedPointi = addedPoints_[patchPointi].last(); - - //Pout<< " For point:" << meshPointi - // << " at:" << mesh_.points()[meshPointi] - // << " at:" << pp.localPoints()[patchPointi] - // << " using addedpoint:" << addedPointi - // << " at:" << meshMod.points()[addedPointi] - // << endl; + // Duplicated point. Might not be present on processor? + //const label patchPointi = pp.meshPointMap()[meshPointi]; + const auto pointFnd = pp.meshPointMap().find(meshPointi); + if (pointFnd) + { + const label patchPointi = pointFnd(); + const label addedPointi = + addedPoints_[patchPointi].last(); + + // Adapt vertices of face + newFace[fp] = addedPointi; - newFace[fp] = addedPointi; + // Store displacement for syncing later on + displacement[meshPointi] = + meshMod.points()[addedPointi] + -mesh_.points()[meshPointi]; + } } } @@ -2879,6 +2990,90 @@ void Foam::addPatchCellLayer::setRefinement zoneFlip // face flip in zone ); } + + + //Pout<< "Adapted point on existing face:" << facei + // << " at:" << mesh_.faceCentres()[facei] + // << " zone:" << zoneID << nl + // << " old:" << f + // << " n:" << newFace.unitNormal(meshMod.points()) + // << " coords:" << UIndirectList<point>(mesh_.points(), f) + // << nl + // << " new:" << newFace + // << " n:" << newFace.unitNormal(meshMod.points()) + // << " coords:" + // << UIndirectList<point>(meshMod.points(), newFace) + // << endl; + } + + if (!extrude_) + { + // Bit tricky: + // - if intruding we're modifying the vertices on existing + // coupled faces (in extrude mode we're adding additional points/ + // faces where we're already doing the right thing) + // - if the other side face is itself not on an extrudePatch + // it will not be adapted so you'll get a coupled point mismatch + // - so send over the new point position (or use some map on + // coupled faces to describe which vertex on the face is changed?) + // - and adapt any point that was marked but not an extrusion point + // - to test: take 2x2 blockMesh and remove one of the cells. Put + // all 3 cells on different processors. Now one of the processors + // will be coupled to extrude/intruded points but itself have + // no patch. + + syncTools::syncPointList + ( + mesh_, + displacement, + maxMagSqrEqOp<vector>(), + vector::zero + ); + + forAll(baseFaces, facei) + { + const face& f = mesh_.faces()[facei]; + const face& baseF = baseFaces[facei]; + + if (isBlockedFace(facei) || baseF.empty()) + { + // Either part of patch or no duplicated points on face + continue; + } + + // Start off from original face + forAll(f, fp) + { + if (baseF[fp] != labelMax) + { + // Duplicated point. Might not be present on processor? + const label meshPointi = f[fp]; + + if + ( + !pp.meshPointMap().found(meshPointi) + && displacement[meshPointi] != vector::zero + ) + { + const point newPt + ( + mesh_.points()[meshPointi] + + displacement[meshPointi] + ); + meshMod.modifyPoint + ( + meshPointi, + newPt, + mesh_.pointZones().whichZone(meshPointi), + true + ); + + // Unmark as being done + displacement[meshPointi] = Zero; + } + } + } + } } } } diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H index fe59a8f6080b4e613f750a0b203e090a98ced277..f42b63b7b7c72faf3643535d8236805602d9cb09 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2015 OpenFOAM Foundation - Copyright (C) 2020 OpenCFD Ltd. + Copyright (C) 2020,2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,8 +28,8 @@ Class Foam::addPatchCellLayer Description - Adds layers of cells to outside of polyPatch. Can optionally create - stand-alone extruded mesh (addToMesh=false). + Adds layers of cells to outside (or inside) of polyMesh. Can optionally + create stand-alone extruded mesh (addToMesh=false). Call setRefinement with offset vector for every patch point and number of layers per patch face and number of layers per patch point. @@ -51,24 +51,47 @@ Description \verbatim Was: - a b <- patch of boundary face - +------+------+ - | | | <- original cells - +------+------+ - - Becomes: - - a b <- patch of boundary face - +------+------+ - + +------+ - +------+------+ - +------+------+ - | | | <- original cells - +------+------+ + a b <- patch of boundary face + +------+------+ + | | | <- original cells + +------+------+ + + Extrusion: + + + added added + face cell + ---- ---- + a b <- patch of boundary face + +------+------+ 3 + | | | 2 + + +------+ 2 + | | | 1 + +------+------+ 1 + | | | 0 + +------+------+ 0 + | | | <- original cells + +------+------+ + + + + Intrusion: + + face cell + ---- ---- + a b <- patch of boundary face + +------+------+ 0 + | | | 0 + +------+------+ 1 + | | | 1 + + +------+ 2 + | | | 2 + +------+------+ 3 + | | | <- original cells + +------+------+ \endverbatim - - added faces get same patchID as face they are extruded from - 'side' faces (i.e. on the edge of pp) get the patchID/zoneID of the other patch/zone they are connected to (hopefully only 1) @@ -78,10 +101,15 @@ Description \verbatim a b b <- patch of boundary face +------+------+------+ + | | | | | | | | <- cells + | | | | +------+------+------+ + (shown for extrusion mode only): + + ^ ^ <- wanted extrusion vector (none at far right) a | b | b <- patch of boundary face +------+------+------+ @@ -123,7 +151,7 @@ class primitiveMesh; class globalIndex; /*---------------------------------------------------------------------------*\ - Class addPatchCellLayer Declaration + Class addPatchCellLayer Declaration \*---------------------------------------------------------------------------*/ class addPatchCellLayer @@ -136,15 +164,23 @@ class addPatchCellLayer //- Add layers to existing mesh or create new mesh const bool addToMesh_; + //- Add layers to outside of mesh or to inside + const bool extrude_; + + //- For all patchpoints: list of added points (size 0 or nLayers) // First point in list is one nearest to original point in patch, - // last one is the new point on the surface. + // last one is + // - extrude : the new point on the surface + // - intrude : the point connecting to the original cell. labelListList addedPoints_; //- For all patchfaces: list of layer faces. // - empty if no face extruded - // - first face is original boundary face - // - last one is new boundary face. + // - first element is original patch face + // - last element is + // - extrude : the new boundary face + // - intrude : the new internal face to the original cell. labelListList layerFaces_; @@ -226,25 +262,6 @@ class addPatchCellLayer label& inflateFaceI ); - //- Find internal or boundary face to get extrude properties - // from. zoneFlip consistent with ppEdge ordering - static void findZoneFace - ( - const bool useInternalFaces, - const bool useBoundaryFaces, - - const polyMesh& mesh, - const indirectPrimitivePatch& pp, - const label ppEdgeI, - const labelUIndList& excludeFaces, - const labelList& meshFaces, - - label& inflateFaceI, - label& patchI, - label& zoneI, - bool& zoneFlip - ); - //- Mark internal and boundary edges of patch. In mesh edges //- since processor might not have pp but does have edge. static void markPatchEdges @@ -294,7 +311,12 @@ public: // Constructors //- Construct from mesh. - explicit addPatchCellLayer(const polyMesh&, const bool addToMesh=true); + explicit addPatchCellLayer + ( + const polyMesh&, + const bool addToMesh=true, + const bool extrude=true + ); // Member Functions @@ -314,8 +336,13 @@ public: } //- Helper: get added cells per patch face. - // addedCells[patchFace] is list of cells added. Last element is - // the top cells (i.e. the boundary cell) + // addedCells[patchFace] is list of cells added. + // extrude : + // first element : next to original cell + // last element : is the top cell (i.e. the boundary cell) + // intrude : + // first element : top cell + // last element : next to original cell static labelListList addedCells ( const polyMesh&,