From 65e925f1fef7a39aaf1f296c96b7fb2018c3bfad Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 21 Aug 2009 17:36:27 +0100 Subject: [PATCH] extrudeMesh does extrusion from existing mesh --- .../mesh/generation/extrudeMesh/Make/options | 4 + .../mesh/generation/extrudeMesh/extrudeMesh.C | 548 ++++++++++++++--- .../extrudeModel/extrudeModel/extrudeModel.C | 24 + .../extrudeModel/extrudeModel/extrudeModel.H | 8 + .../extrudeModel/linearNormal/linearNormal.C | 3 +- .../extrudeModel/linearNormal/linearNormal.H | 2 +- .../extrudeModel/linearRadial/linearRadial.C | 3 +- .../extrudeModel/linearRadial/linearRadial.H | 2 +- .../extrudeModel/sigmaRadial/sigmaRadial.C | 8 +- .../extrudeModel/sigmaRadial/sigmaRadial.H | 2 +- .../extrudeMesh/extrudeModel/wedge/wedge.C | 4 +- .../extrudeMesh/extrudeModel/wedge/wedge.H | 2 +- .../generation/extrudeMesh/extrudeProperties | 39 +- .../perfectInterface/perfectInterface.C | 577 +++++++++--------- .../perfectInterface/perfectInterface.H | 14 +- .../polyTopoChange/addPatchCellLayer.C | 30 +- .../polyTopoChange/addPatchCellLayer.H | 24 +- .../polyTopoChange/polyTopoChange.C | 5 + 18 files changed, 870 insertions(+), 429 deletions(-) diff --git a/applications/utilities/mesh/generation/extrudeMesh/Make/options b/applications/utilities/mesh/generation/extrudeMesh/Make/options index ce0e27f401..1cdc68a2c7 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/Make/options +++ b/applications/utilities/mesh/generation/extrudeMesh/Make/options @@ -1,10 +1,14 @@ EXE_INC = \ -IextrudedMesh \ -IextrudeModel/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude EXE_LIBS = \ + -lfiniteVolume \ + -lsurfMesh \ -lmeshTools \ -ldynamicMesh \ -lextrudeModel diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C index d149f85923..4d5bca0e82 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeMesh.C @@ -36,13 +36,15 @@ Description #include "Time.H" #include "dimensionedTypes.H" #include "IFstream.H" -#include "faceMesh.H" #include "polyTopoChange.H" #include "polyTopoChanger.H" #include "edgeCollapser.H" #include "mathematicalConstants.H" #include "globalMeshData.H" #include "perfectInterface.H" +#include "addPatchCellLayer.H" +#include "fvMesh.H" +#include "MeshedSurfaces.H" #include "extrudedMesh.H" #include "extrudeModel.H" @@ -50,14 +52,148 @@ Description using namespace Foam; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +enum ExtrudeMode +{ + MESH, + PATCH, + SURFACE +}; + +template<> +const char* NamedEnum::names[] = +{ + "mesh", + "patch", + "surface" +}; +static const NamedEnum ExtrudeModeNames; + + +void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName) +{ + // Create dummy system/fv* + { + IOobject io + ( + "fvSchemes", + mesh.time().system(), + regionName, + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ); + + Info<< "Testing:" << io.objectPath() << endl; + + if (!io.headerOk()) + { + Info<< "Writing dummy " << regionName/io.name() << endl; + dictionary dummyDict; + dictionary divDict; + dummyDict.add("divSchemes", divDict); + dictionary gradDict; + dummyDict.add("gradSchemes", gradDict); + dictionary laplDict; + dummyDict.add("laplacianSchemes", laplDict); + + IOdictionary(io, dummyDict).regIOobject::write(); + } + } + { + IOobject io + ( + "fvSolution", + mesh.time().system(), + regionName, + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ); + + if (!io.headerOk()) + { + Info<< "Writing dummy " << regionName/io.name() << endl; + dictionary dummyDict; + IOdictionary(io, dummyDict).regIOobject::write(); + } + } +} + + +label findPatchID(const polyBoundaryMesh& patches, const word& name) +{ + label patchID = patches.findPatchID(name); + + if (patchID == -1) + { + FatalErrorIn("findPatchID(const polyBoundaryMesh&, const word&)") + << "Cannot find patch " << name + << " in the source mesh.\n" + << "Valid patch names are " << patches.names() + << exit(FatalError); + } + return patchID; +} + + +labelList patchFaces(const polyBoundaryMesh& patches, const word& name) +{ + label patchID = findPatchID(patches, name); + const polyPatch& pp = patches[patchID]; + + return identity(pp.size()) + pp.start(); +} + + +void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels) +{ + const labelList& reverseMap = map.reverseFaceMap(); + + labelList newFaceLabels(faceLabels.size()); + label newI = 0; + + forAll(faceLabels, i) + { + label oldFaceI = faceLabels[i]; + + if (reverseMap[oldFaceI] >= 0) + { + newFaceLabels[newI++] = reverseMap[oldFaceI]; + } + } + newFaceLabels.setSize(newI); + faceLabels.transfer(newFaceLabels); +} + + + // Main program: int main(int argc, char *argv[]) { + #include "addRegionOption.H" #include "setRootCase.H" #include "createTimeExtruded.H" - autoPtr meshPtr(NULL); + // Get optional regionName + word regionName; + word regionDir; + if (args.optionReadIfPresent("region", regionName)) + { + regionDir = regionName; + Info<< "Create mesh " << regionName << " for time = " + << runTimeExtruded.timeName() << nl << endl; + } + else + { + regionName = fvMesh::defaultRegion; + Info<< "Create mesh for time = " + << runTimeExtruded.timeName() << nl << endl; + } + IOdictionary dict ( @@ -65,26 +201,44 @@ int main(int argc, char *argv[]) ( "extrudeProperties", runTimeExtruded.constant(), + regionDir, runTimeExtruded, IOobject::MUST_READ ) ); + // Point generator autoPtr model(extrudeModel::New(dict)); - const word sourceType(dict.lookup("constructFrom")); + // Whether to flip normals + const Switch flipNormals(dict.lookup("flipNormals")); - autoPtr fMesh; + // What to extrude + const ExtrudeMode mode = ExtrudeModeNames.read + ( + dict.lookup("constructFrom") + ); - if (sourceType == "patch") + + // Generated mesh (one of either) + autoPtr meshFromMesh; + autoPtr meshFromSurface; + + // Faces on front and back for stitching (in case of mergeFaces) + word frontPatchName; + labelList frontPatchFaces; + word backPatchName; + labelList backPatchFaces; + + if (mode == PATCH || mode == MESH) { fileName sourceCasePath(dict.lookup("sourceCase")); sourceCasePath.expand(); fileName sourceRootDir = sourceCasePath.path(); fileName sourceCaseDir = sourceCasePath.name(); - word patchName(dict.lookup("sourcePatch")); + dict.lookup("sourcePatch") >> frontPatchName; - Info<< "Extruding patch " << patchName + Info<< "Extruding patch " << frontPatchName << " on mesh " << sourceCasePath << nl << endl; @@ -94,31 +248,183 @@ int main(int argc, char *argv[]) sourceRootDir, sourceCaseDir ); - #include "createPolyMesh.H" + #include "createMesh.H" - label patchID = mesh.boundaryMesh().findPatchID(patchName); + const polyBoundaryMesh& patches = mesh.boundaryMesh(); - if (patchID == -1) + + // Topo change container. Either copy an existing mesh or start + // with empty storage (number of patches only needed for checking) + autoPtr 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)); + + indirectPrimitivePatch extrudePatch + ( + IndirectList + ( + mesh.faces(), + patchFaces(patches, frontPatchName) + ), + mesh.points() + ); + + // Only used for addPatchCellLayer into new mesh + labelList exposedPatchIDs; + if (mode == PATCH) { - FatalErrorIn(args.executable()) - << "Cannot find patch " << patchName - << " in the source mesh.\n" - << "Valid patch names are " << mesh.boundaryMesh().names() - << exit(FatalError); + dict.lookup("exposedPatchName") >> backPatchName; + exposedPatchIDs.setSize + ( + extrudePatch.size(), + findPatchID(patches, backPatchName) + ); } - const polyPatch& pp = mesh.boundaryMesh()[patchID]; - fMesh.reset(new faceMesh(pp.localFaces(), pp.localPoints())); + pointField layer0Points(extrudePatch.nPoints()); + pointField displacement(extrudePatch.nPoints()); + forAll(displacement, pointI) { - fileName surfName(runTime.path()/patchName + ".sMesh"); - Info<< "Writing patch as surfaceMesh to " - << surfName << nl << endl; - OFstream os(surfName); - os << fMesh() << nl; + const vector& patchNormal = extrudePatch.pointNormals()[pointI]; + + // layer0 point + layer0Points[pointI] = model() + ( + extrudePatch.localPoints()[pointI], + patchNormal, + 0 + ); + // layerN point + point extrudePt = model() + ( + extrudePatch.localPoints()[pointI], + patchNormal, + model().nLayers() + ); + displacement[pointI] = extrudePt - layer0Points[pointI]; } + + if (flipNormals) + { + Info<< "Flipping faces." << nl << endl; + displacement = -displacement; + } + + // Check if wedge (has layer0 different from original patch points) + // If so move the mesh to starting position. + if (gMax(mag(layer0Points-extrudePatch.localPoints())) > SMALL) + { + Info<< "Moving mesh to layer0 points since differ from original" + << " points - this can happen for wedge extrusions." << nl + << endl; + + pointField newPoints(mesh.points()); + forAll(extrudePatch.meshPoints(), i) + { + newPoints[extrudePatch.meshPoints()[i]] = layer0Points[i]; + } + mesh.movePoints(newPoints); + } + + + // Layers per face + labelList nFaceLayers(extrudePatch.size(), model().nLayers()); + // Layers per point + labelList nPointLayers(extrudePatch.nPoints(), model().nLayers()); + // Displacement for first layer + vectorField firstLayerDisp = displacement*model().sumThickness(1); + // Expansion ratio not used. + scalarField ratio(extrudePatch.nPoints(), 1.0); + + layerExtrude.setRefinement + ( + ratio, // expansion ratio + extrudePatch, // patch faces to extrude + exposedPatchIDs, // if new mesh: patches for exposed faces + nFaceLayers, + nPointLayers, + firstLayerDisp, + meshMod() + ); + + // Reset points according to extrusion model + forAll(layerExtrude.addedPoints(), pointI) + { + const labelList& pPoints = layerExtrude.addedPoints()[pointI]; + forAll(pPoints, pPointI) + { + label meshPointI = pPoints[pPointI]; + + point modelPt + ( + model() + ( + extrudePatch.localPoints()[pointI], + extrudePatch.pointNormals()[pointI], + pPointI+1 // layer + ) + ); + + const_cast&> + ( + meshMod().points() + )[meshPointI] = modelPt; + } + } + + // Store faces on front and exposed patch (if mode=patch there are + // only added faces so cannot used map to old faces) + const labelListList& layerFaces = layerExtrude.layerFaces(); + backPatchFaces.setSize(layerFaces.size()); + frontPatchFaces.setSize(layerFaces.size()); + forAll(backPatchFaces, i) + { + backPatchFaces[i] = layerFaces[i][0]; + frontPatchFaces[i] = layerFaces[i][layerFaces[i].size()-1]; + } + + // Create dummy fvSchemes, fvSolution + createDummyFvMeshFiles(mesh, regionDir); + + // Create actual mesh from polyTopoChange container + autoPtr map = meshMod().makeMesh + ( + meshFromMesh, + IOobject + ( + regionName, + runTimeExtruded.constant(), + runTimeExtruded, + IOobject::NO_READ, + IOobject::AUTO_WRITE, + false + ), + mesh + ); + + // Calculate face labels for front and back. + frontPatchFaces = renumber + ( + map().reverseFaceMap(), + frontPatchFaces + ); + backPatchFaces = renumber + ( + map().reverseFaceMap(), + backPatchFaces + ); } - else if (sourceType == "surface") + else { // Read from surface fileName surfName(dict.lookup("surface")); @@ -126,52 +432,62 @@ int main(int argc, char *argv[]) Info<< "Extruding surfaceMesh read from file " << surfName << nl << endl; - IFstream is(surfName); + MeshedSurface fMesh(surfName); - fMesh.reset(new faceMesh(is)); - - Info<< "Read patch from file " << surfName << nl - << endl; - } - else - { - FatalErrorIn(args.executable()) - << "Illegal 'constructFrom' specification. Should either be " - << "patch or surface." << exit(FatalError); - } - - Switch flipNormals(dict.lookup("flipNormals")); - - if (flipNormals) - { - Info<< "Flipping faces." << nl << endl; - - faceList faces(fMesh().size()); - forAll(faces, i) + if (flipNormals) { - faces[i] = fMesh()[i].reverseFace(); + Info<< "Flipping faces." << nl << endl; + faceList& faces = const_cast(fMesh.faces()); + forAll(faces, i) + { + faces[i] = fMesh[i].reverseFace(); + } } - fMesh.reset(new faceMesh(faces, fMesh().localPoints())); + + Info<< "Extruding surface with :" << nl + << " points : " << fMesh.points().size() << nl + << " faces : " << fMesh.size() << nl + << " normals[0] : " << fMesh.faceNormals()[0] + << nl + << endl; + + meshFromSurface.reset + ( + new extrudedMesh + ( + IOobject + ( + extrudedMesh::defaultRegion, + runTimeExtruded.constant(), + runTimeExtruded + ), + fMesh, + model() + ) + ); + + + // Get the faces on front and back + frontPatchName = "originalPatch"; + frontPatchFaces = patchFaces + ( + meshFromSurface().boundaryMesh(), + frontPatchName + ); + backPatchName = "otherSide"; + backPatchFaces = patchFaces + ( + meshFromSurface().boundaryMesh(), + backPatchName + ); } - Info<< "Extruding patch with :" << nl - << " points : " << fMesh().points().size() << nl - << " faces : " << fMesh().size() << nl - << " normals[0] : " << fMesh().faceNormals()[0] - << nl - << endl; - - extrudedMesh mesh + polyMesh& mesh = ( - IOobject - ( - extrudedMesh::defaultRegion, - runTimeExtruded.constant(), - runTimeExtruded - ), - fMesh(), - model() + meshFromMesh.valid() + ? meshFromMesh() + : meshFromSurface() ); @@ -184,17 +500,6 @@ int main(int argc, char *argv[]) << "Merge distance : " << mergeDim << nl << endl; - const polyBoundaryMesh& patches = mesh.boundaryMesh(); - - const label origPatchID = patches.findPatchID("originalPatch"); - const label otherPatchID = patches.findPatchID("otherSide"); - - if (origPatchID == -1 || otherPatchID == -1) - { - FatalErrorIn(args.executable()) - << "Cannot find patch originalPatch or otherSide." << nl - << "Valid patches are " << patches.names() << exit(FatalError); - } // Collapse edges // ~~~~~~~~~~~~~~ @@ -237,6 +542,10 @@ int main(int argc, char *argv[]) // Update fields mesh.updateMesh(map); + // Update stored data + updateFaceLabels(map(), frontPatchFaces); + updateFaceLabels(map(), backPatchFaces); + // Move mesh (if inflation used) if (map().hasMotionPoints()) { @@ -252,22 +561,33 @@ int main(int argc, char *argv[]) Switch mergeFaces(dict.lookup("mergeFaces")); if (mergeFaces) { + if (mode == MESH) + { + FatalErrorIn(args.executable()) + << "Cannot stitch front and back of extrusion since" + << " in 'mesh' mode (extrusion appended to mesh)." + << exit(FatalError); + } + Info<< "Assuming full 360 degree axisymmetric case;" << " stitching faces on patches " - << patches[origPatchID].name() << " and " - << patches[otherPatchID].name() << " together ..." << nl << endl; + << frontPatchName << " and " + << backPatchName << " together ..." << nl << endl; + + if (frontPatchFaces.size() != backPatchFaces.size()) + { + FatalErrorIn(args.executable()) + << "Differing number of faces on front (" + << frontPatchFaces.size() << ") and back (" + << backPatchFaces.size() << ")" + << exit(FatalError); + } + + polyTopoChanger stitcher(mesh); stitcher.setSize(1); - // Make list of masterPatch faces - labelList isf(patches[origPatchID].size()); - - forAll (isf, i) - { - isf[i] = patches[origPatchID].start() + i; - } - const word cutZoneName("originalCutFaceZone"); List fz @@ -276,8 +596,8 @@ int main(int argc, char *argv[]) new faceZone ( cutZoneName, - isf, - boolList(isf.size(), false), + frontPatchFaces, + boolList(frontPatchFaces.size(), false), 0, mesh.faceZones() ) @@ -286,26 +606,58 @@ int main(int argc, char *argv[]) mesh.addZones(List(0), fz, List(0)); // Add the perfect interface mesh modifier - stitcher.set + perfectInterface perfectStitcher ( + "couple", 0, - new perfectInterface - ( - "couple", - 0, - stitcher, - cutZoneName, - patches[origPatchID].name(), - patches[otherPatchID].name() - ) + stitcher, + cutZoneName, + word::null, // dummy patch name + word::null // dummy patch name ); - // Execute all polyMeshModifiers - autoPtr morphMap = stitcher.changeMesh(true); + // Topo change container + polyTopoChange meshMod(mesh); - mesh.movePoints(morphMap->preMotionPoints()); + perfectStitcher.setRefinement + ( + indirectPrimitivePatch + ( + IndirectList + ( + mesh.faces(), + frontPatchFaces + ), + mesh.points() + ), + indirectPrimitivePatch + ( + IndirectList + ( + mesh.faces(), + backPatchFaces + ), + mesh.points() + ), + meshMod + ); + + // Construct new mesh from polyTopoChange. + autoPtr map = meshMod.changeMesh(mesh, false); + + // Update fields + mesh.updateMesh(map); + + // Move mesh (if inflation used) + if (map().hasMotionPoints()) + { + mesh.movePoints(map().preMotionPoints()); + } } + mesh.setInstance(runTimeExtruded.constant()); + Info<< "Writing mesh to " << mesh.objectPath() << nl << endl; + if (!mesh.write()) { FatalErrorIn(args.executable()) << "Failed writing mesh" diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C index 7305431b96..5550c56b15 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.C @@ -43,6 +43,7 @@ Foam::extrudeModel::extrudeModel ) : nLayers_(readLabel(dict.lookup("nLayers"))), + expansionRatio_(readScalar(dict.lookup("expansionRatio"))), dict_(dict), coeffDict_(dict.subDict(modelType + "Coeffs")) {} @@ -62,5 +63,28 @@ Foam::label Foam::extrudeModel::nLayers() const } +Foam::scalar Foam::extrudeModel::expansionRatio() const +{ + return expansionRatio_; +} + + +Foam::scalar Foam::extrudeModel::sumThickness(const label layer) const +{ + // 1+r+r^2+ .. +r^(n-1) = (1-r^n)/(1-r) + + if (mag(1.0-expansionRatio_) < SMALL) + { + return scalar(layer)/nLayers_; + } + else + { + return + (1.0-pow(expansionRatio_, layer)) + / (1.0-pow(expansionRatio_, nLayers_)); + } +} + + // ************************************************************************* // diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H index d8f3699fd5..ac96bbeea6 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H @@ -58,6 +58,8 @@ protected: const label nLayers_; + const scalar expansionRatio_; + const dictionary& dict_; const dictionary& coeffDict_; @@ -113,9 +115,15 @@ public: label nLayers() const; + scalar expansionRatio() const; + // Member Operators + //- Helper: calculate cumulative relative thickness for layer. + // (layer=0 -> 0; layer=nLayers -> 1) + scalar sumThickness(const label layer) const; + virtual point operator() ( const point& surfacePoint, diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C index c36098ab1e..cd42e91816 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C @@ -72,7 +72,8 @@ point linearNormal::operator() const label layer ) const { - scalar d = thickness_*layer/nLayers_; + //scalar d = thickness_*layer/nLayers_; + scalar d = thickness_*sumThickness(layer); return surfacePoint + d*surfaceNormal; } diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H index 949547311b..2dbcdccb16 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H @@ -64,7 +64,7 @@ public: // Constructors - //- Construct from components + //- Construct from dictionary linearNormal(const dictionary& dict); diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C index e8d406b110..7d4fcff9ba 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C @@ -67,8 +67,7 @@ point linearRadial::operator() scalar rs = mag(surfacePoint); vector rsHat = surfacePoint/rs; - scalar delta = (R_ - rs)/nLayers_; - scalar r = rs + layer*delta; + scalar r = rs + (R_ - rs)*sumThickness(layer); return r*rsHat; } diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H index 0fda12e1b1..66758c1b12 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H @@ -61,7 +61,7 @@ public: // Constructors - //- Construct from components + //- Construct from dictionary linearRadial(const dictionary& dict); diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C index fc6f040305..d35ad4e539 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C @@ -49,7 +49,13 @@ sigmaRadial::sigmaRadial(const dictionary& dict) RTbyg_(readScalar(coeffDict_.lookup("RTbyg"))), pRef_(readScalar(coeffDict_.lookup("pRef"))), pStrat_(readScalar(coeffDict_.lookup("pStrat"))) -{} +{ + if (mag(expansionRatio() - 1.0) > SMALL) + { + WarningIn("sigmaRadial::sigmaRadial(const dictionary&)") + << "Ignoring expansionRatio setting." << endl; + } +} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H index 525e89d11c..c74ad02457 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H @@ -63,7 +63,7 @@ public: // Constructors - //- Construct from components + //- Construct from dictionary sigmaRadial(const dictionary& dict); diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C index 3d2c883ea1..5d0a3621b4 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C @@ -88,8 +88,8 @@ point wedge::operator() } else { - //sliceAngle = angle_*(layer + 1)/nLayers_; - sliceAngle = angle_*layer/nLayers_; + //sliceAngle = angle_*layer/nLayers_; + sliceAngle = angle_*sumThickness(layer); } // Find projection onto axis (or rather decompose surfacePoint diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H index 39c4b4a2d3..18707c32d1 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H @@ -77,7 +77,7 @@ public: // Constructors - //- Construct from components + //- Construct from dictionary wedge(const dictionary& dict); diff --git a/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties b/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties index 45618034de..e1283ab8b5 100644 --- a/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties +++ b/applications/utilities/mesh/generation/extrudeMesh/extrudeProperties @@ -14,24 +14,28 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Where to get surface from: either from surface ('surface') or -// from (flipped) patch of existing case ('patch') -constructFrom patch; //surface; +// What to extrude: +// patch : from patch of another case ('sourceCase') +// mesh : as above but with original case included +// surface : from externally read surface -// If construct from (flipped) patch -sourceCase "$FOAM_RUN/icoFoam/cavity"; +//constructFrom mesh; +constructFrom patch; +//constructFrom surface; + +// If construct from patch/mesh: +sourceCase "../cavity"; sourcePatch movingWall; +// If construct from patch: patch to use for back (can be same as sourcePatch) +exposedPatchName movingWall; +// If construct from surface: +surface "movingWall.stl"; + + // Flip surface normals before usage. flipNormals false; -// If construct from surface -surface "movingWall.sMesh"; - - -// Do front and back need to be merged? Usually only makes sense for 360 -// degree wedges. -mergeFaces true; //- Linear extrusion in point-normal direction @@ -46,11 +50,13 @@ extrudeModel wedge; //- Extrudes into sphere with grading according to pressure (atmospherics) //extrudeModel sigmaRadial; -nLayers 20; +nLayers 10; + +expansionRatio 1.0; //0.9; wedgeCoeffs { - axisPt (0 0.1 0); + axisPt (0 0.1 -0.05); axis (-1 0 0); angle 360; // For nLayers=1 assume symmetry so angle/2 on each side } @@ -73,5 +79,10 @@ sigmaRadialCoeffs } +// Do front and back need to be merged? Usually only makes sense for 360 +// degree wedges. +mergeFaces true; + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/dynamicMesh/perfectInterface/perfectInterface.C b/src/dynamicMesh/perfectInterface/perfectInterface.C index 80c9ed6583..ee892b7e87 100644 --- a/src/dynamicMesh/perfectInterface/perfectInterface.C +++ b/src/dynamicMesh/perfectInterface/perfectInterface.C @@ -38,6 +38,7 @@ Description #include "polyModifyFace.H" #include "polyRemovePoint.H" #include "polyRemoveFace.H" +#include "indirectPrimitivePatch.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -61,7 +62,7 @@ const Foam::scalar Foam::perfectInterface::tol_ = 1E-3; Foam::pointField Foam::perfectInterface::calcFaceCentres ( - const primitivePatch& pp + const indirectPrimitivePatch& pp ) { const pointField& points = pp.points(); @@ -155,6 +156,295 @@ bool Foam::perfectInterface::changeTopology() const } +void Foam::perfectInterface::setRefinement +( + const indirectPrimitivePatch& pp0, + const indirectPrimitivePatch& pp1, + polyTopoChange& ref +) const +{ + const polyMesh& mesh = topoChanger().mesh(); + + const polyBoundaryMesh& patches = mesh.boundaryMesh(); + + // Some aliases + const edgeList& edges0 = pp0.edges(); + const pointField& pts0 = pp0.localPoints(); + const pointField& pts1 = pp1.localPoints(); + const labelList& meshPts0 = pp0.meshPoints(); + const labelList& meshPts1 = pp1.meshPoints(); + + + // Get local dimension as fraction of minimum edge length + + scalar minLen = GREAT; + + forAll(edges0, edgeI) + { + minLen = min(minLen, edges0[edgeI].mag(pts0)); + } + scalar typDim = tol_*minLen; + + if (debug) + { + Pout<< "typDim:" << typDim << " edges0:" << edges0.size() + << " pts0:" << pts0.size() << " pts1:" << pts1.size() + << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl; + } + + + // Determine pointMapping in mesh point labels. Uses geometric + // comparison to find correspondence between patch points. + + labelList renumberPoints(mesh.points().size()); + forAll(renumberPoints, i) + { + renumberPoints[i] = i; + } + { + labelList from1To0Points(pts1.size()); + + bool matchOk = matchPoints + ( + pts1, + pts0, + scalarField(pts1.size(), typDim), // tolerance + true, // verbose + from1To0Points + ); + + if (!matchOk) + { + FatalErrorIn + ( + "perfectInterface::setRefinement(polyTopoChange& ref) const" + ) << "Points on patch sides do not match to within tolerance " + << typDim << exit(FatalError); + } + + forAll(pts1, i) + { + renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]]; + } + } + + + + // Calculate correspondence between patch faces + + labelList from0To1Faces(pp1.size()); + + bool matchOk = matchPoints + ( + calcFaceCentres(pp0), + calcFaceCentres(pp1), + scalarField(pp0.size(), typDim), // tolerance + true, // verbose + from0To1Faces + ); + + if (!matchOk) + { + FatalErrorIn + ( + "perfectInterface::setRefinement(polyTopoChange& ref) const" + ) << "Face centres of patch sides do not match to within tolerance " + << typDim << exit(FatalError); + } + + + + // Now + // - renumber faces using pts1 (except patch1 faces) + // - remove patch1 faces. Remember cell label on owner side. + // - modify patch0 faces to be internal. + + // 1. Get faces to be renumbered + labelHashSet affectedFaces(2*pp1.size()); + forAll(meshPts1, i) + { + label meshPointI = meshPts1[i]; + + if (meshPointI != renumberPoints[meshPointI]) + { + const labelList& pFaces = mesh.pointFaces()[meshPointI]; + + forAll(pFaces, pFaceI) + { + affectedFaces.insert(pFaces[pFaceI]); + } + } + } + forAll(pp1, i) + { + affectedFaces.erase(pp1.addressing()[i]); + } + // Remove patch0 from renumbered faces. Should not be nessecary since + // patch0 and 1 should not share any point (if created by mergeMeshing) + // so affectedFaces should not contain any patch0 faces but you can + // never be sure what the user is doing. + forAll(pp0, i) + { + label faceI = pp0.addressing()[i]; + + if (affectedFaces.erase(faceI)) + { + WarningIn + ( + "perfectInterface::setRefinement(polyTopoChange&) const" + ) << "Found face " << faceI << " vertices " + << mesh.faces()[faceI] << " whose points are" + << " used both by master patch and slave patch" << endl; + } + } + + + // 2. Renumber (non patch0/1) faces. + for + ( + labelHashSet::const_iterator iter = affectedFaces.begin(); + iter != affectedFaces.end(); + ++iter + ) + { + label faceI = iter.key(); + + const face& f = mesh.faces()[faceI]; + + face newFace(f.size()); + + forAll(newFace, fp) + { + newFace[fp] = renumberPoints[f[fp]]; + } + + label nbr = -1; + + label patchI = -1; + + if (mesh.isInternalFace(faceI)) + { + nbr = mesh.faceNeighbour()[faceI]; + } + else + { + patchI = patches.whichPatch(faceI); + } + + label zoneID = mesh.faceZones().whichZone(faceI); + + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = mesh.faceZones()[zoneID]; + + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + ref.setAction + ( + polyModifyFace + ( + newFace, // modified face + faceI, // label of face being modified + mesh.faceOwner()[faceI], // owner + nbr, // neighbour + false, // face flip + patchI, // patch for face + false, // remove from zone + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + } + + + // 3. Remove patch1 points + forAll(meshPts1, i) + { + label meshPointI = meshPts1[i]; + + if (meshPointI != renumberPoints[meshPointI]) + { + ref.setAction(polyRemovePoint(meshPointI)); + } + } + + + // 4. Remove patch1 faces + forAll(pp1, i) + { + label faceI = pp1.addressing()[i]; + ref.setAction(polyRemoveFace(faceI)); + } + + + // 5. Modify patch0 faces for new points (not really nessecary; see + // comment above about patch1 and patch0 never sharing points) and + // becoming internal. + const boolList& mfFlip = + mesh.faceZones()[faceZoneID_.index()].flipMap(); + + forAll(pp0, i) + { + label faceI = pp0.addressing()[i]; + + const face& f = mesh.faces()[faceI]; + + face newFace(f.size()); + + forAll(newFace, fp) + { + newFace[fp] = renumberPoints[f[fp]]; + } + + label own = mesh.faceOwner()[faceI]; + + label pp1FaceI = pp1.addressing()[from0To1Faces[i]]; + + label nbr = mesh.faceOwner()[pp1FaceI]; + + if (own < nbr) + { + ref.setAction + ( + polyModifyFace + ( + newFace, // modified face + faceI, // label of face being modified + own, // owner + nbr, // neighbour + false, // face flip + -1, // patch for face + false, // remove from zone + faceZoneID_.index(), // zone for face + mfFlip[i] // face flip in zone + ) + ); + } + else + { + ref.setAction + ( + polyModifyFace + ( + newFace.reverseFace(), // modified face + faceI, // label of face being modified + nbr, // owner + own, // neighbour + true, // face flip + -1, // patch for face + false, // remove from zone + faceZoneID_.index(), // zone for face + !mfFlip[i] // face flip in zone + ) + ); + } + } +} + + void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const { if (debug) @@ -176,286 +466,25 @@ void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const const polyMesh& mesh = topoChanger().mesh(); const polyBoundaryMesh& patches = mesh.boundaryMesh(); - - const polyPatch& pp0 = patches[masterPatchID_.index()]; - const polyPatch& pp1 = patches[slavePatchID_.index()]; - - // Some aliases - const edgeList& edges0 = pp0.edges(); - const pointField& pts0 = pp0.localPoints(); - const pointField& pts1 = pp1.localPoints(); - const labelList& meshPts0 = pp0.meshPoints(); - const labelList& meshPts1 = pp1.meshPoints(); - - - // Get local dimension as fraction of minimum edge length - - scalar minLen = GREAT; - - forAll(edges0, edgeI) - { - minLen = min(minLen, edges0[edgeI].mag(pts0)); - } - scalar typDim = tol_*minLen; - - if (debug) - { - Pout<< "typDim:" << typDim << " edges0:" << edges0.size() - << " pts0:" << pts0.size() << " pts1:" << pts1.size() - << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl; - } + const polyPatch& patch0 = patches[masterPatchID_.index()]; + const polyPatch& patch1 = patches[slavePatchID_.index()]; - // Determine pointMapping in mesh point labels. Uses geometric - // comparison to find correspondence between patch points. - - labelList renumberPoints(mesh.points().size()); - forAll(renumberPoints, i) - { - renumberPoints[i] = i; - } - { - labelList from1To0Points(pts1.size()); - - bool matchOk = matchPoints - ( - pts1, - pts0, - scalarField(pts1.size(), typDim), // tolerance - true, // verbose - from1To0Points - ); - - if (!matchOk) - { - FatalErrorIn - ( - "perfectInterface::setRefinement(polyTopoChange& ref) const" - ) << "Points on patches " << pp0.name() << " and " - << pp1.name() << " do not match to within tolerance " - << typDim << exit(FatalError); - } - - forAll(pts1, i) - { - renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]]; - } - } - - - - // Calculate correspondence between patch faces - - labelList from0To1Faces(pp1.size()); - - bool matchOk = matchPoints + labelList pp0Labels(identity(patch0.size())+patch0.start()); + indirectPrimitivePatch pp0 ( - calcFaceCentres(pp0), - calcFaceCentres(pp1), - scalarField(pp0.size(), typDim), // tolerance - true, // verbose - from0To1Faces + IndirectList(mesh.faces(), pp0Labels), + mesh.points() ); - if (!matchOk) - { - FatalErrorIn - ( - "perfectInterface::setRefinement(polyTopoChange& ref) const" - ) << "Face centres of patches " << pp0.name() << " and " - << pp1.name() << " do not match to within tolerance " << typDim - << exit(FatalError); - } - - - - // Now - // - renumber faces using pts1 (except patch1 faces) - // - remove patch1 faces. Remember cell label on owner side. - // - modify patch0 faces to be internal. - - // 1. Get faces to be renumbered - labelHashSet affectedFaces(2*pp1.size()); - forAll(meshPts1, i) - { - label meshPointI = meshPts1[i]; - - if (meshPointI != renumberPoints[meshPointI]) - { - const labelList& pFaces = mesh.pointFaces()[meshPointI]; - - forAll(pFaces, pFaceI) - { - affectedFaces.insert(pFaces[pFaceI]); - } - } - } - forAll(pp1, i) - { - affectedFaces.erase(pp1.start() + i); - } - // Remove patch0 from renumbered faces. Should not be nessecary since - // patch0 and 1 should not share any point (if created by mergeMeshing) - // so affectedFaces should not contain any patch0 faces but you can - // never be sure what the user is doing. - forAll(pp0, i) - { - if (affectedFaces.erase(pp0.start() + i)) - { - WarningIn - ( - "perfectInterface::setRefinement(polyTopoChange&) const" - ) << "Found face " << pp0.start() + i << " vertices " - << mesh.faces()[pp0.start() + i] << " whose points are" - << " used both by master patch " << pp0.name() - << " and slave patch " << pp1.name() - << endl; - } - } - - - // 2. Renumber (non patch0/1) faces. - for + labelList pp1Labels(identity(patch1.size())+patch1.start()); + indirectPrimitivePatch pp1 ( - labelHashSet::const_iterator iter = affectedFaces.begin(); - iter != affectedFaces.end(); - ++iter - ) - { - label faceI = iter.key(); + IndirectList(mesh.faces(), pp1Labels), + mesh.points() + ); - const face& f = mesh.faces()[faceI]; - - face newFace(f.size()); - - forAll(newFace, fp) - { - newFace[fp] = renumberPoints[f[fp]]; - } - - label nbr = -1; - - label patchI = -1; - - if (mesh.isInternalFace(faceI)) - { - nbr = mesh.faceNeighbour()[faceI]; - } - else - { - patchI = patches.whichPatch(faceI); - } - - label zoneID = mesh.faceZones().whichZone(faceI); - - bool zoneFlip = false; - - if (zoneID >= 0) - { - const faceZone& fZone = mesh.faceZones()[zoneID]; - - zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; - } - - ref.setAction - ( - polyModifyFace - ( - newFace, // modified face - faceI, // label of face being modified - mesh.faceOwner()[faceI], // owner - nbr, // neighbour - false, // face flip - patchI, // patch for face - false, // remove from zone - zoneID, // zone for face - zoneFlip // face flip in zone - ) - ); - } - - - // 3. Remove patch1 points - forAll(meshPts1, i) - { - label meshPointI = meshPts1[i]; - - if (meshPointI != renumberPoints[meshPointI]) - { - ref.setAction(polyRemovePoint(meshPointI)); - } - } - - - // 4. Remove patch1 faces - forAll(pp1, i) - { - ref.setAction(polyRemoveFace(pp1.start() + i)); - } - - - // 5. Modify patch0 faces for new points (not really nessecary; see - // comment above about patch1 and patch0 never sharing points) and - // becoming internal. - const boolList& mfFlip = - mesh.faceZones()[faceZoneID_.index()].flipMap(); - - forAll(pp0, i) - { - label faceI = pp0.start() + i; - - const face& f = mesh.faces()[faceI]; - - face newFace(f.size()); - - forAll(newFace, fp) - { - newFace[fp] = renumberPoints[f[fp]]; - } - - label own = mesh.faceOwner()[faceI]; - - label pp1FaceI = pp1.start() + from0To1Faces[i]; - - label nbr = mesh.faceOwner()[pp1FaceI]; - - if (own < nbr) - { - ref.setAction - ( - polyModifyFace - ( - newFace, // modified face - faceI, // label of face being modified - own, // owner - nbr, // neighbour - false, // face flip - -1, // patch for face - false, // remove from zone - faceZoneID_.index(), // zone for face - mfFlip[i] // face flip in zone - ) - ); - } - else - { - ref.setAction - ( - polyModifyFace - ( - newFace.reverseFace(), // modified face - faceI, // label of face being modified - nbr, // owner - own, // neighbour - false, // face flip - -1, // patch for face - false, // remove from zone - faceZoneID_.index(), // zone for face - !mfFlip[i] // face flip in zone - ) - ); - } - } + setRefinement(pp0, pp1, ref); } } diff --git a/src/dynamicMesh/perfectInterface/perfectInterface.H b/src/dynamicMesh/perfectInterface/perfectInterface.H index 3916bb3e74..dec4f6a48d 100644 --- a/src/dynamicMesh/perfectInterface/perfectInterface.H +++ b/src/dynamicMesh/perfectInterface/perfectInterface.H @@ -40,6 +40,7 @@ SourceFiles #include "polyMeshModifier.H" #include "polyPatchID.H" #include "ZoneIDs.H" +#include "indirectPrimitivePatch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -75,7 +76,7 @@ class perfectInterface // Private Member Functions //- Calculate face centres on patch - static pointField calcFaceCentres(const primitivePatch&); + static pointField calcFaceCentres(const indirectPrimitivePatch&); //- Disallow default bitwise copy construct @@ -128,6 +129,17 @@ public: // into the topological change virtual void setRefinement(polyTopoChange&) const; + //- Insert the layer addition/removal instructions + // into the topological change. Uses only mesh, not any of the + // patch and zone indices. Bit of a workaround - used in extruding + // a mesh. + virtual void setRefinement + ( + const indirectPrimitivePatch& pp0, + const indirectPrimitivePatch& pp1, + polyTopoChange& + ) const; + //- Modify motion points to comply with the topological change virtual void modifyMotionPoints(pointField& motionPoints) const; diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C index c3b104134c..1a725f0ebf 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C @@ -558,6 +558,7 @@ void Foam::addPatchCellLayer::setRefinement ( const scalarField& expansionRatio, const indirectPrimitivePatch& pp, + const labelList& exposedPatchID, const labelList& nFaceLayers, const labelList& nPointLayers, const vectorField& firstLayerDisp, @@ -973,11 +974,6 @@ void Foam::addPatchCellLayer::setRefinement ownZoneI // zone for cell ) ); - - //Pout<< "For patchFace:" << patchFaceI - // << " meshFace:" << pp.addressing()[patchFaceI] - // << " layer:" << i << " added cell:" - // << addedCells[patchFaceI][i] << endl; } } } @@ -998,7 +994,6 @@ void Foam::addPatchCellLayer::setRefinement if (addedCells[patchFaceI].size()) { layerFaces_[patchFaceI].setSize(addedCells[patchFaceI].size() + 1); - layerFaces_[patchFaceI][0] = meshFaceI; label zoneI = mesh_.faceZones().whichZone(meshFaceI); @@ -1059,18 +1054,13 @@ void Foam::addPatchCellLayer::setRefinement nei, // neighbour -1, // master point -1, // master edge - meshFaceI, // master face for addition + (addToMesh_ ? meshFaceI : -1), // master face false, // flux flip patchI, // patch for face zoneI, // zone for face false // face zone flip ) ); - //Pout<< "Added inbetween face " << newFace - // << " own:" << addedCells[patchFaceI][i] - // << " nei:" << nei - // << " patch:" << patchI - // << endl; } } } @@ -1079,7 +1069,6 @@ void Foam::addPatchCellLayer::setRefinement // Modify old patch faces to be on the inside // - labelList copiedPatchFaces; if (addToMesh_) { forAll(pp, patchFaceI) @@ -1088,6 +1077,8 @@ void Foam::addPatchCellLayer::setRefinement { label meshFaceI = pp.addressing()[patchFaceI]; + layerFaces_[patchFaceI][0] = meshFaceI; + label zoneI = mesh_.faceZones().whichZone(meshFaceI); meshMod.setAction @@ -1105,17 +1096,13 @@ void Foam::addPatchCellLayer::setRefinement false // face flip in zone ) ); - //Pout<< "Modified old patch face " << meshFaceI - // << " own:" << mesh_.faceOwner()[meshFaceI] - // << " nei:" << addedCells[patchFaceI][0] - // << endl; } } } else { - // If creating new mesh: copy existing patch points - copiedPatchFaces.setSize(pp.size()); + // If creating new mesh: reverse original faces and put them + // in the exposed patch ID. forAll(pp, patchFaceI) { if (nFaceLayers[patchFaceI] > 0) @@ -1136,7 +1123,7 @@ void Foam::addPatchCellLayer::setRefinement f[fp] = copiedPatchPoints[f[fp]]; } - copiedPatchFaces[patchFaceI] = meshMod.setAction + layerFaces_[patchFaceI][0] = meshMod.setAction ( polyAddFace ( @@ -1147,7 +1134,7 @@ void Foam::addPatchCellLayer::setRefinement -1, // masterEdge -1, // masterFace true, // face flip - patchID[patchFaceI], // patch for face + exposedPatchID[patchFaceI], // patch for face zoneI, // zone for face zoneFlip // face flip in zone ) @@ -1158,7 +1145,6 @@ void Foam::addPatchCellLayer::setRefinement - // // Create 'side' faces, one per edge that is being extended. // diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H index f950ca2619..f4b5bb2e79 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H @@ -295,8 +295,10 @@ public: //- Play commands into polyTopoChange to create layers on top // of indirectPrimitivePatch (have to be outside faces). // Gets displacement per patch point. - // - nPointLayers : number of layers per (patch)point - // - nFaceLayers : number of layers per (patch) face + // - exposedPatchID : only used if creating a new mesh (addToMesh=false) + // gives per pp face the patch the exposed face should get. + // - nPointLayers : number of layers per (patch)point. + // - nFaceLayers : number of layers per (patch) face. // - firstDisplacement : displacement per point for first // layer of points (i.e. nearest to original mesh). If zero // do not add point. @@ -309,14 +311,15 @@ public: // get a cell should firstDisplacement be <> 0 // Note: cells get added from owner cells of patch faces // (instead of e.g. from patch faces) - void setRefinement - ( - const scalarField& expansionRatio, - const indirectPrimitivePatch& pp, - const labelList& nFaceLayers, - const labelList& nPointLayers, - const vectorField& firstLayerDisp, - polyTopoChange& meshMod + void setRefinement + ( + const scalarField& expansionRatio, + const indirectPrimitivePatch& pp, + const labelList& exposedPatchID, + const labelList& nFaceLayers, + const labelList& nPointLayers, + const vectorField& firstLayerDisp, + polyTopoChange& meshMod ); @@ -333,6 +336,7 @@ public: ( scalarField(pp.nPoints(), 1.0), // expansion ration pp, + labelList(0), labelList(pp.size(), nLayers), labelList(pp.nPoints(), nLayers), overallDisplacement / nLayers, diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C index 849481c2eb..6b1874dc4b 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C @@ -749,6 +749,11 @@ void Foam::polyTopoChange::getFaceOrder " const" ) << "Did not determine new position" << " for face " << faceI + << " owner " << faceOwner_[faceI] + << " neighbour " << faceNeighbour_[faceI] + << " region " << region_[faceI] << endl + << "This is usually caused by not specifying a patch for" + << " a boundary face." << abort(FatalError); } }