diff --git a/applications/utilities/mesh/generation/extrudeMesh/Make/options b/applications/utilities/mesh/generation/extrudeMesh/Make/options
index ce0e27f401b48c65f68ffed57b4edb710b8401a7..1cdc68a2c7d9ca4dde6664f2a567caf945ee6660 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 d149f85923cb03ebd38eadbbc3f0eafb8d4ccaf6..4d5bca0e8220e2c8691ea0e914e870c6581165ea 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<ExtrudeMode, 3>::names[] =
+{
+ "mesh",
+ "patch",
+ "surface"
+};
+static const NamedEnum<ExtrudeMode, 3> 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<extrudedMesh> 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<extrudeModel> model(extrudeModel::New(dict));
- const word sourceType(dict.lookup("constructFrom"));
+ // Whether to flip normals
+ const Switch flipNormals(dict.lookup("flipNormals"));
+
+ // What to extrude
+ const ExtrudeMode mode = ExtrudeModeNames.read
+ (
+ dict.lookup("constructFrom")
+ );
+
- autoPtr<faceMesh> fMesh;
+ // Generated mesh (one of either)
+ autoPtr<fvMesh> meshFromMesh;
+ autoPtr<polyMesh> meshFromSurface;
- if (sourceType == "patch")
+ // 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"
+
+ const polyBoundaryMesh& patches = mesh.boundaryMesh();
- label patchID = mesh.boundaryMesh().findPatchID(patchName);
- if (patchID == -1)
+ // 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));
+
+ indirectPrimitivePatch extrudePatch
+ (
+ IndirectList<face>
+ (
+ 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<DynamicList<point>&>
+ (
+ 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<mapPolyMesh> 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<face> fMesh(surfName);
- fMesh.reset(new faceMesh(is));
+ if (flipNormals)
+ {
+ Info<< "Flipping faces." << nl << endl;
+ faceList& faces = const_cast<faceList&>(fMesh.faces());
+ forAll(faces, i)
+ {
+ faces[i] = fMesh[i].reverseFace();
+ }
+ }
- Info<< "Read patch from file " << surfName << nl
- << endl;
- }
- else
- {
- FatalErrorIn(args.executable())
- << "Illegal 'constructFrom' specification. Should either be "
- << "patch or surface." << exit(FatalError);
- }
+ Info<< "Extruding surface with :" << nl
+ << " points : " << fMesh.points().size() << nl
+ << " faces : " << fMesh.size() << nl
+ << " normals[0] : " << fMesh.faceNormals()[0]
+ << nl
+ << endl;
- Switch flipNormals(dict.lookup("flipNormals"));
+ meshFromSurface.reset
+ (
+ new extrudedMesh
+ (
+ IOobject
+ (
+ extrudedMesh::defaultRegion,
+ runTimeExtruded.constant(),
+ runTimeExtruded
+ ),
+ fMesh,
+ model()
+ )
+ );
- if (flipNormals)
- {
- Info<< "Flipping faces." << nl << endl;
- faceList faces(fMesh().size());
- forAll(faces, i)
- {
- faces[i] = fMesh()[i].reverseFace();
- }
- fMesh.reset(new faceMesh(faces, fMesh().localPoints()));
+ // 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;
-
- polyTopoChanger stitcher(mesh);
- stitcher.setSize(1);
-
- // Make list of masterPatch faces
- labelList isf(patches[origPatchID].size());
+ << frontPatchName << " and "
+ << backPatchName << " together ..." << nl << endl;
- forAll (isf, i)
+ if (frontPatchFaces.size() != backPatchFaces.size())
{
- isf[i] = patches[origPatchID].start() + i;
+ FatalErrorIn(args.executable())
+ << "Differing number of faces on front ("
+ << frontPatchFaces.size() << ") and back ("
+ << backPatchFaces.size() << ")"
+ << exit(FatalError);
}
+
+
+ polyTopoChanger stitcher(mesh);
+ stitcher.setSize(1);
+
const word cutZoneName("originalCutFaceZone");
List<faceZone*> 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<pointZone*>(0), fz, List<cellZone*>(0));
// Add the perfect interface mesh modifier
- stitcher.set
+ perfectInterface perfectStitcher
(
+ "couple",
0,
- new perfectInterface
+ stitcher,
+ cutZoneName,
+ word::null, // dummy patch name
+ word::null // dummy patch name
+ );
+
+ // Topo change container
+ polyTopoChange meshMod(mesh);
+
+ perfectStitcher.setRefinement
+ (
+ indirectPrimitivePatch
(
- "couple",
- 0,
- stitcher,
- cutZoneName,
- patches[origPatchID].name(),
- patches[otherPatchID].name()
- )
+ IndirectList<face>
+ (
+ mesh.faces(),
+ frontPatchFaces
+ ),
+ mesh.points()
+ ),
+ indirectPrimitivePatch
+ (
+ IndirectList<face>
+ (
+ mesh.faces(),
+ backPatchFaces
+ ),
+ mesh.points()
+ ),
+ meshMod
);
- // Execute all polyMeshModifiers
- autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
+ // Construct new mesh from polyTopoChange.
+ autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+ // Update fields
+ mesh.updateMesh(map);
- mesh.movePoints(morphMap->preMotionPoints());
+ // 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 7305431b96971e591123ccaf5f7696528486db1f..5550c56b15249e78e55725d39914a3c6b1ea4fb8 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 d8f3699fd5a2d0d3fd8334444b0a63d4ef981dd0..ac96bbeea6324e18d81a3ed42a7ebb161c03045a 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 c36098ab1e6a3cf10b748a84e7062524ad9e3fd9..cd42e918169e0431d3c35ac84deed960bccc7711 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 949547311b39d8c8cb191f69be04d0cf7c34e948..2dbcdccb162592ee75691124905b1d48bdb4c2a1 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 e8d406b110ad12350b00f733de5f0ad30ac40771..7d4fcff9ba02e8a348be0baa89787155dbcdb7fb 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 0fda12e1b14e8c49c18d7d24cfee56ff9ed5fd8c..66758c1b12897e64ebc01d4d98b519b4ed1c3262 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 fc6f04030504f04e57a14c7d5d447a0a3e754510..d35ad4e539447972ef603a8a9469a6276b140990 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 525e89d11c94aa96a31ec11b65916b4c9ca06627..c74ad024575998d06364129abf3a1cf24e70c0c2 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 3d2c883ea1901edb740f8afdd7de3c9b634858a2..5d0a3621b43c8752a1da56dfd950a41900f51e2f 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 39c4b4a2d3073ade350e8c5db80f6faea29b5d5b..18707c32d11a560fde7ec582741318ec89bce799 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 45618034deb13b6a1e5c55bfc02e7bedd9928203..e1283ab8b59d0ba2220a815ee375f242d1e452f4 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";
+// Flip surface normals before usage.
+flipNormals false;
-// 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/applications/utilities/mesh/generation/extrudeMesh/faceMesh.H b/applications/utilities/mesh/generation/extrudeMesh/faceMesh.H
deleted file mode 100644
index c8e1075f52e5f31ceba4f8930ad6be4dbab9158e..0000000000000000000000000000000000000000
--- a/applications/utilities/mesh/generation/extrudeMesh/faceMesh.H
+++ /dev/null
@@ -1,114 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
- \\/ M anipulation |
--------------------------------------------------------------------------------
-License
- This file is part of OpenFOAM.
-
- OpenFOAM is free software; you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the
- Free Software Foundation; either version 2 of the License, or (at your
- option) any later version.
-
- OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- for more details.
-
- You should have received a copy of the GNU General Public License
- along with OpenFOAM; if not, write to the Free Software Foundation,
- Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
- Foam::faceMesh
-
-Description
- Storage for surface mesh i.e. points and faces.
-
-SourceFiles
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef faceMesh_H
-#define faceMesh_H
-
-#include "PrimitivePatch.H"
-#include "faceList.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
- Class faceMesh Declaration
-\*---------------------------------------------------------------------------*/
-
-class faceMesh
-:
- public PrimitivePatch<face, List, pointField>
-{
- // Private member functions
-
- PrimitivePatch<face, List, pointField> read(Istream& is)
- {
- pointField points(is);
- faceList faces(is);
- return PrimitivePatch<face, Foam::List, pointField>(faces, points);
- }
-
-
-public:
-
- // Constructors
-
- //- Construct from components
- faceMesh(const faceList& faces, const pointField& points)
- :
- PrimitivePatch<face, Foam::List, pointField>(faces, points)
- {}
-
- //- Construct from Istream
- faceMesh(Istream& is)
- :
- PrimitivePatch<face, Foam::List, pointField>(read(is))
- {}
-
-
- // Member Functions
-
- void flip()
- {
- forAll(*this, i)
- {
- face& f = operator[](i);
- f = f.reverseFace();
- }
- clearOut();
- }
-
- // IOstream Operators
-
- friend Ostream& operator<<(Ostream& os, const faceMesh& fm)
- {
- return os
- << fm.points()
- << token::NL
- << static_cast<PrimitivePatch<face, Foam::List, pointField> >
- (fm);
- }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatch.C b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
index 7d74e7208c2e552d4bf0cc3af9d05b19ec8c64ef..40ba0efadebc02dec26e5de383e6fef0634a6c93 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatch.C
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatch.C
@@ -208,8 +208,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
// Dump halves
{
OFstream str(prefix+cycPatch.name()+"_half0.obj");
- Pout<< "Dumping cycPatch.name() half0 faces to " << str.name()
- << endl;
+ Pout<< "Dumping " << cycPatch.name()
+ << " half0 faces to " << str.name() << endl;
meshTools::writeOBJ
(
str,
@@ -226,8 +226,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
}
{
OFstream str(prefix+cycPatch.name()+"_half1.obj");
- Pout<< "Dumping cycPatch.name() half1 faces to " << str.name()
- << endl;
+ Pout<< "Dumping " << cycPatch.name()
+ << " half1 faces to " << str.name() << endl;
meshTools::writeOBJ
(
str,
@@ -563,7 +563,7 @@ int main(int argc, char *argv[])
dumpCyclicMatch("initial_", mesh);
// Read patch construct info from dictionary
- PtrList<dictionary> patchSources(dict.lookup("patches"));
+ PtrList<dictionary> patchSources(dict.lookup("patchInfo"));
diff --git a/applications/utilities/mesh/manipulation/createPatch/createPatchDict b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
index 9704392a10ea95e23f5afbe02a9c45bbbc736406..b8676f134d3798461fd5d0a7b8bed1af21151d4b 100644
--- a/applications/utilities/mesh/manipulation/createPatch/createPatchDict
+++ b/applications/utilities/mesh/manipulation/createPatch/createPatchDict
@@ -46,7 +46,7 @@ matchTolerance 1E-3;
pointSync true;
// Patches to create.
-patches
+patchInfo
(
{
// Name of new patch
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index fea25e35eb8a6e64bfd01564cda14d42299d1403..716bdc6ddfacb967f4e97c2a2e7481327bb6af86 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2987,6 +2987,7 @@ void Foam::autoLayerDriver::addLayers
(
invExpansionRatio,
pp,
+ labelList(0), // exposed patchIDs, not used for adding layers
nPatchFaceLayers, // layers per face
nPatchPointLayers, // layers per point
firstDisp, // thickness of layer nearest internal mesh
diff --git a/src/dynamicMesh/perfectInterface/perfectInterface.C b/src/dynamicMesh/perfectInterface/perfectInterface.C
index 80c9ed65836a807e8d87988cb6738d49320ad3dc..ee892b7e873dee29a6ee6fc89e7a9aad1c50cdfa 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,307 +156,335 @@ bool Foam::perfectInterface::changeTopology() const
}
-void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
+void Foam::perfectInterface::setRefinement
+(
+ const indirectPrimitivePatch& pp0,
+ const indirectPrimitivePatch& pp1,
+ polyTopoChange& ref
+) const
{
- if (debug)
+ 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)
{
- Pout<< "bool perfectInterface::setRefinement(polyTopoChange&) const : "
- << "for object " << name() << " : "
- << "masterPatchID_:" << masterPatchID_
- << " slavePatchID_:" << slavePatchID_
- << " faceZoneID_:" << faceZoneID_ << endl;
+ minLen = min(minLen, edges0[edgeI].mag(pts0));
}
+ scalar typDim = tol_*minLen;
- if
- (
- masterPatchID_.active()
- && slavePatchID_.active()
- && faceZoneID_.active()
- )
+ if (debug)
{
- const polyMesh& mesh = topoChanger().mesh();
-
- const polyBoundaryMesh& patches = mesh.boundaryMesh();
+ Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
+ << " pts0:" << pts0.size() << " pts1:" << pts1.size()
+ << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
+ }
- 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();
-
+ // Determine pointMapping in mesh point labels. Uses geometric
+ // comparison to find correspondence between patch points.
- // Get local dimension as fraction of minimum edge length
+ labelList renumberPoints(mesh.points().size());
+ forAll(renumberPoints, i)
+ {
+ renumberPoints[i] = i;
+ }
+ {
+ labelList from1To0Points(pts1.size());
- scalar minLen = GREAT;
+ bool matchOk = matchPoints
+ (
+ pts1,
+ pts0,
+ scalarField(pts1.size(), typDim), // tolerance
+ true, // verbose
+ from1To0Points
+ );
- forAll(edges0, edgeI)
+ if (!matchOk)
{
- minLen = min(minLen, edges0[edgeI].mag(pts0));
+ FatalErrorIn
+ (
+ "perfectInterface::setRefinement(polyTopoChange& ref) const"
+ ) << "Points on patch sides do not match to within tolerance "
+ << typDim << exit(FatalError);
}
- scalar typDim = tol_*minLen;
- if (debug)
+ forAll(pts1, i)
{
- Pout<< "typDim:" << typDim << " edges0:" << edges0.size()
- << " pts0:" << pts0.size() << " pts1:" << pts1.size()
- << " pp0:" << pp0.size() << " pp1:" << pp1.size() << endl;
+ renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
}
+ }
- // 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());
+ // Calculate correspondence between patch faces
- bool matchOk = matchPoints
- (
- pts1,
- pts0,
- scalarField(pts1.size(), typDim), // tolerance
- true, // verbose
- from1To0Points
- );
+ labelList from0To1Faces(pp1.size());
- 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);
- }
+ bool matchOk = matchPoints
+ (
+ calcFaceCentres(pp0),
+ calcFaceCentres(pp1),
+ scalarField(pp0.size(), typDim), // tolerance
+ true, // verbose
+ from0To1Faces
+ );
- forAll(pts1, i)
- {
- renumberPoints[meshPts1[i]] = meshPts0[from1To0Points[i]];
- }
- }
+ if (!matchOk)
+ {
+ FatalErrorIn
+ (
+ "perfectInterface::setRefinement(polyTopoChange& ref) const"
+ ) << "Face centres of patch sides do not match to within tolerance "
+ << typDim << exit(FatalError);
+ }
- // Calculate correspondence between patch faces
+ // Now
+ // - renumber faces using pts1 (except patch1 faces)
+ // - remove patch1 faces. Remember cell label on owner side.
+ // - modify patch0 faces to be internal.
- labelList from0To1Faces(pp1.size());
+ // 1. Get faces to be renumbered
+ labelHashSet affectedFaces(2*pp1.size());
+ forAll(meshPts1, i)
+ {
+ label meshPointI = meshPts1[i];
- bool matchOk = matchPoints
- (
- calcFaceCentres(pp0),
- calcFaceCentres(pp1),
- scalarField(pp0.size(), typDim), // tolerance
- true, // verbose
- from0To1Faces
- );
+ if (meshPointI != renumberPoints[meshPointI])
+ {
+ const labelList& pFaces = mesh.pointFaces()[meshPointI];
- if (!matchOk)
+ 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))
{
- FatalErrorIn
+ WarningIn
(
- "perfectInterface::setRefinement(polyTopoChange& ref) const"
- ) << "Face centres of patches " << pp0.name() << " and "
- << pp1.name() << " do not match to within tolerance " << typDim
- << exit(FatalError);
+ "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();
- // Now
- // - renumber faces using pts1 (except patch1 faces)
- // - remove patch1 faces. Remember cell label on owner side.
- // - modify patch0 faces to be internal.
+ const face& f = mesh.faces()[faceI];
- // 1. Get faces to be renumbered
- labelHashSet affectedFaces(2*pp1.size());
- forAll(meshPts1, i)
+ face newFace(f.size());
+
+ forAll(newFace, fp)
{
- label meshPointI = meshPts1[i];
+ newFace[fp] = renumberPoints[f[fp]];
+ }
- if (meshPointI != renumberPoints[meshPointI])
- {
- const labelList& pFaces = mesh.pointFaces()[meshPointI];
+ label nbr = -1;
- forAll(pFaces, pFaceI)
- {
- affectedFaces.insert(pFaces[pFaceI]);
- }
- }
- }
- forAll(pp1, i)
+ label patchI = -1;
+
+ if (mesh.isInternalFace(faceI))
{
- affectedFaces.erase(pp1.start() + i);
+ nbr = mesh.faceNeighbour()[faceI];
}
- // 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)
+ else
{
- 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;
- }
+ patchI = patches.whichPatch(faceI);
}
+ label zoneID = mesh.faceZones().whichZone(faceI);
+
+ bool zoneFlip = false;
- // 2. Renumber (non patch0/1) faces.
- for
+ if (zoneID >= 0)
+ {
+ const faceZone& fZone = mesh.faceZones()[zoneID];
+
+ zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+ }
+
+ ref.setAction
(
- labelHashSet::const_iterator iter = affectedFaces.begin();
- iter != affectedFaces.end();
- ++iter
- )
+ 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])
{
- label faceI = iter.key();
+ ref.setAction(polyRemovePoint(meshPointI));
+ }
+ }
- const face& f = mesh.faces()[faceI];
- face newFace(f.size());
+ // 4. Remove patch1 faces
+ forAll(pp1, i)
+ {
+ label faceI = pp1.addressing()[i];
+ ref.setAction(polyRemoveFace(faceI));
+ }
- forAll(newFace, fp)
- {
- newFace[fp] = renumberPoints[f[fp]];
- }
- label nbr = -1;
+ // 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];
- label patchI = -1;
+ const face& f = mesh.faces()[faceI];
- if (mesh.isInternalFace(faceI))
- {
- nbr = mesh.faceNeighbour()[faceI];
- }
- else
- {
- patchI = patches.whichPatch(faceI);
- }
+ face newFace(f.size());
- label zoneID = mesh.faceZones().whichZone(faceI);
+ forAll(newFace, fp)
+ {
+ newFace[fp] = renumberPoints[f[fp]];
+ }
- bool zoneFlip = false;
+ label own = mesh.faceOwner()[faceI];
- if (zoneID >= 0)
- {
- const faceZone& fZone = mesh.faceZones()[zoneID];
+ label pp1FaceI = pp1.addressing()[from0To1Faces[i]];
- zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
- }
+ label nbr = mesh.faceOwner()[pp1FaceI];
+ if (own < nbr)
+ {
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
+ 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
)
);
}
-
-
- // 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)
+ else
{
- ref.setAction(polyRemoveFace(pp1.start() + i));
+ 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
+ )
+ );
}
+ }
+}
- // 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];
+void Foam::perfectInterface::setRefinement(polyTopoChange& ref) const
+{
+ if (debug)
+ {
+ Pout<< "bool perfectInterface::setRefinement(polyTopoChange&) const : "
+ << "for object " << name() << " : "
+ << "masterPatchID_:" << masterPatchID_
+ << " slavePatchID_:" << slavePatchID_
+ << " faceZoneID_:" << faceZoneID_ << endl;
+ }
- face newFace(f.size());
+ if
+ (
+ masterPatchID_.active()
+ && slavePatchID_.active()
+ && faceZoneID_.active()
+ )
+ {
+ const polyMesh& mesh = topoChanger().mesh();
- forAll(newFace, fp)
- {
- newFace[fp] = renumberPoints[f[fp]];
- }
+ const polyBoundaryMesh& patches = mesh.boundaryMesh();
+ const polyPatch& patch0 = patches[masterPatchID_.index()];
+ const polyPatch& patch1 = patches[slavePatchID_.index()];
- label own = mesh.faceOwner()[faceI];
- label pp1FaceI = pp1.start() + from0To1Faces[i];
+ labelList pp0Labels(identity(patch0.size())+patch0.start());
+ indirectPrimitivePatch pp0
+ (
+ IndirectList<face>(mesh.faces(), pp0Labels),
+ mesh.points()
+ );
- label nbr = mesh.faceOwner()[pp1FaceI];
+ labelList pp1Labels(identity(patch1.size())+patch1.start());
+ indirectPrimitivePatch pp1
+ (
+ IndirectList<face>(mesh.faces(), pp1Labels),
+ mesh.points()
+ );
- 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 3916bb3e74b6bfd647ddbb5a46a5b116a2d2910b..dec4f6a48d0448c894aff74d84410381acf0d78c 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 2f842e6d66127b62fb226d79899f1b2fd7876b3a..1a725f0ebfe6979252cc8b720f05feede1287e93 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
@@ -335,17 +335,19 @@ Foam::label Foam::addPatchCellLayer::addSideFace
label inflateEdgeI = -1;
// Check mesh faces using edge
- forAll(meshFaces, i)
+ if (addToMesh_)
{
- if (mesh_.isInternalFace(meshFaces[i]))
+ forAll(meshFaces, i)
{
- // meshEdge uses internal faces so ok to inflate from it
- inflateEdgeI = meshEdgeI;
- break;
+ if (mesh_.isInternalFace(meshFaces[i]))
+ {
+ // meshEdge uses internal faces so ok to inflate from it
+ inflateEdgeI = meshEdgeI;
+ break;
+ }
}
}
-
// Get my mesh face and its zone.
label meshFaceI = pp.addressing()[ownFaceI];
label zoneI = mesh_.faceZones().whichZone(meshFaceI);
@@ -504,9 +506,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from mesh
-Foam::addPatchCellLayer::addPatchCellLayer(const polyMesh& mesh)
+Foam::addPatchCellLayer::addPatchCellLayer
+(
+ const polyMesh& mesh,
+ const bool addToMesh
+)
:
mesh_(mesh),
+ addToMesh_(addToMesh),
addedPoints_(0),
layerFaces_(0)
{}
@@ -551,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,
@@ -827,6 +835,19 @@ void Foam::addPatchCellLayer::setRefinement
}
+ const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+ // Precalculated patchID for each patch face
+ labelList patchID(pp.size());
+
+ forAll(pp, patchFaceI)
+ {
+ label meshFaceI = pp.addressing()[patchFaceI];
+
+ patchID[patchFaceI] = patches.whichPatch(meshFaceI);
+ }
+
+
// From master point (in patch point label) to added points (in mesh point
// label)
addedPoints_.setSize(pp.nPoints());
@@ -857,6 +878,33 @@ void Foam::addPatchCellLayer::setRefinement
// Create new points
//
+ // If creating new mesh: copy existing patch points
+ labelList copiedPatchPoints;
+ if (!addToMesh_)
+ {
+ copiedPatchPoints.setSize(firstLayerDisp.size());
+ forAll(firstLayerDisp, patchPointI)
+ {
+ if (addedPoints_[patchPointI].size())
+ {
+ label meshPointI = meshPoints[patchPointI];
+ label zoneI = mesh_.pointZones().whichZone(meshPointI);
+ copiedPatchPoints[patchPointI] = meshMod.setAction
+ (
+ polyAddPoint
+ (
+ mesh_.points()[meshPointI], // point
+ -1, // master point
+ zoneI, // zone for point
+ true // supports a cell
+ )
+ );
+ }
+ }
+ }
+
+
+ // Create points for additional layers
forAll(firstLayerDisp, patchPointI)
{
if (addedPoints_[patchPointI].size())
@@ -878,7 +926,7 @@ void Foam::addPatchCellLayer::setRefinement
polyAddPoint
(
pt, // point
- meshPointI, // master point
+ (addToMesh_ ? meshPointI : -1), // master point
zoneI, // zone for point
true // supports a cell
)
@@ -922,34 +970,15 @@ void Foam::addPatchCellLayer::setRefinement
-1, // master point
-1, // master edge
-1, // master face
- mesh_.faceOwner()[meshFaceI], // master cell id
+ (addToMesh_ ? mesh_.faceOwner()[meshFaceI] : -1),//master
ownZoneI // zone for cell
)
);
-
- //Pout<< "For patchFace:" << patchFaceI
- // << " meshFace:" << pp.addressing()[patchFaceI]
- // << " layer:" << i << " added cell:"
- // << addedCells[patchFaceI][i]
- // << endl;
}
}
}
- const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
- // Precalculated patchID for each patch face
- labelList patchID(pp.size());
-
- forAll(pp, patchFaceI)
- {
- label meshFaceI = pp.addressing()[patchFaceI];
-
- patchID[patchFaceI] = patches.whichPatch(meshFaceI);
- }
-
-
// Create faces on top of the original patch faces.
// These faces are created from original patch faces outwards so in order
@@ -965,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);
@@ -981,7 +1009,12 @@ void Foam::addPatchCellLayer::setRefinement
if (addedPoints_[f[fp]].empty())
{
// Keep original point
- newFace[fp] = meshPoints[f[fp]];
+ newFace[fp] =
+ (
+ addToMesh_
+ ? meshPoints[f[fp]]
+ : copiedPatchPoints[f[fp]]
+ );
}
else
{
@@ -1021,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;
}
}
}
@@ -1040,37 +1068,83 @@ void Foam::addPatchCellLayer::setRefinement
//
// Modify old patch faces to be on the inside
//
- forAll(pp, patchFaceI)
+
+ if (addToMesh_)
{
- if (addedCells[patchFaceI].size())
+ forAll(pp, patchFaceI)
{
- label meshFaceI = pp.addressing()[patchFaceI];
+ if (addedCells[patchFaceI].size())
+ {
+ label meshFaceI = pp.addressing()[patchFaceI];
- label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+ layerFaces_[patchFaceI][0] = meshFaceI;
- meshMod.setAction
- (
- polyModifyFace
+ label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+
+ meshMod.setAction
(
- pp[patchFaceI], // modified face
- meshFaceI, // label of face
- mesh_.faceOwner()[meshFaceI], // owner
- addedCells[patchFaceI][0], // neighbour
- false, // face flip
- -1, // patch for face
- false, // remove from zone
- zoneI, // zone for face
- false // face flip in zone
- )
- );
- //Pout<< "Modified old patch face " << meshFaceI
- // << " own:" << mesh_.faceOwner()[meshFaceI]
- // << " nei:" << addedCells[patchFaceI][0]
- // << endl;
+ polyModifyFace
+ (
+ pp[patchFaceI], // modified face
+ meshFaceI, // label of face
+ mesh_.faceOwner()[meshFaceI], // owner
+ addedCells[patchFaceI][0], // neighbour
+ false, // face flip
+ -1, // patch for face
+ false, // remove from zone
+ zoneI, // zone for face
+ false // face flip in zone
+ )
+ );
+ }
+ }
+ }
+ else
+ {
+ // If creating new mesh: reverse original faces and put them
+ // in the exposed patch ID.
+ forAll(pp, patchFaceI)
+ {
+ if (nFaceLayers[patchFaceI] > 0)
+ {
+ label meshFaceI = pp.addressing()[patchFaceI];
+ label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+ bool zoneFlip = false;
+ if (zoneI != -1)
+ {
+ const faceZone& fz = mesh_.faceZones()[zoneI];
+ zoneFlip = !fz.flipMap()[fz.whichFace(meshFaceI)];
+ }
+
+ // Reverse and renumber old patch face.
+ face f(pp.localFaces()[patchFaceI].reverseFace());
+ forAll(f, fp)
+ {
+ f[fp] = copiedPatchPoints[f[fp]];
+ }
+
+ layerFaces_[patchFaceI][0] = meshMod.setAction
+ (
+ polyAddFace
+ (
+ f, // modified face
+ addedCells[patchFaceI][0], // owner
+ -1, // neighbour
+ -1, // masterPoint
+ -1, // masterEdge
+ -1, // masterFace
+ true, // face flip
+ exposedPatchID[patchFaceI], // patch for face
+ zoneI, // zone for face
+ zoneFlip // face flip in zone
+ )
+ );
+ }
}
}
+
//
// Create 'side' faces, one per edge that is being extended.
//
@@ -1255,7 +1329,16 @@ void Foam::addPatchCellLayer::setRefinement
forAll(stringedVerts, stringedI)
{
label v = stringedVerts[stringedI];
- addVertex(meshPoints[v], newFace, newFp);
+ addVertex
+ (
+ (
+ addToMesh_
+ ? meshPoints[v]
+ : copiedPatchPoints[v]
+ ),
+ newFace,
+ newFp
+ );
}
}
else
@@ -1276,7 +1359,16 @@ void Foam::addPatchCellLayer::setRefinement
}
else
{
- addVertex(meshPoints[v], newFace, newFp);
+ addVertex
+ (
+ (
+ addToMesh_
+ ? meshPoints[v]
+ : copiedPatchPoints[v]
+ ),
+ newFace,
+ newFp
+ );
}
}
}
@@ -1316,7 +1408,16 @@ void Foam::addPatchCellLayer::setRefinement
}
else
{
- addVertex(meshPoints[v], newFace, newFp);
+ addVertex
+ (
+ (
+ addToMesh_
+ ? meshPoints[v]
+ : copiedPatchPoints[v]
+ ),
+ newFace,
+ newFp
+ );
}
}
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
index f03d23e2144fb22c5a9f093f221ed5fd86dc4a42..f4b5bb2e79325c6c9023b3b417a053e41a573650 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
@@ -26,7 +26,8 @@ Class
Foam::addPatchCellLayer
Description
- Adds layers of cells to outside of polyPatch.
+ Adds layers of cells to outside of polyPatch. 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.
@@ -164,6 +165,9 @@ class addPatchCellLayer
//- Reference to mesh
const polyMesh& mesh_;
+ //- Add layers to existing mesh or create new mesh
+ const bool addToMesh_;
+
//- 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.
@@ -253,7 +257,7 @@ public:
// Constructors
//- Construct from mesh.
- addPatchCellLayer(const polyMesh& mesh);
+ addPatchCellLayer(const polyMesh& mesh, const bool addToMesh = true);
// Member Functions
@@ -291,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.
@@ -305,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
);
@@ -329,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 f6c397640fce9eb982b3bb125df349896a8582c4..6b1874dc4b846850035734e4a009c5f0b421eb9b 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);
}
}
@@ -3176,7 +3181,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
(
autoPtr<fvMesh>& newMeshPtr,
const IOobject& io,
- const fvMesh& mesh,
+ const polyMesh& mesh,
const bool syncParallel,
const bool orderCells,
const bool orderPoints
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H
index 77801e21f0dc33b703585f15ef3ffbd65a5c70ad..e7a0ba8fb380fd8a258cebf98baa1400440f134e 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H
@@ -585,7 +585,7 @@ public:
(
autoPtr<fvMesh>& newMesh,
const IOobject& io,
- const fvMesh& mesh,
+ const polyMesh& mesh,
const bool syncParallel = true,
const bool orderCells = false,
const bool orderPoints = false
diff --git a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormatCore.C b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormatCore.C
index 887fed0fd9caa67e0ab48929a2678beeb113e0a2..a089b502e34fb121c72d0e198e6dad5c075267ae 100644
--- a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormatCore.C
+++ b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormatCore.C
@@ -48,7 +48,8 @@ int Foam::fileFormats::STLsurfaceFormatCore::detectBINARY
{
off_t dataFileSize = Foam::fileSize(filename);
- istream& is = IFstream(filename, IOstream::BINARY)().stdStream();
+ IFstream str(filename, IOstream::BINARY);
+ istream& is = str().stdStream();
// Read the STL header
char header[headerSize];