fa93ce8cd7
- markup codingStyleGuide.org examples so they actually indent correctly - use 'Info<<' as per codingStyleGuide instead of 'Info <<'
703 lines
19 KiB
C
703 lines
19 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / 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
|
|
|
|
Description
|
|
Extrude mesh from existing patch (by default outwards facing normals;
|
|
optional flips faces) or from patch read from file.
|
|
|
|
Note: Merges close points so be careful.
|
|
|
|
Type of extrusion prescribed by run-time selectable model.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "argList.H"
|
|
#include "Time.H"
|
|
#include "dimensionedTypes.H"
|
|
#include "IFstream.H"
|
|
#include "polyTopoChange.H"
|
|
#include "polyTopoChanger.H"
|
|
#include "edgeCollapser.H"
|
|
#include "globalMeshData.H"
|
|
#include "perfectInterface.H"
|
|
#include "addPatchCellLayer.H"
|
|
#include "fvMesh.H"
|
|
#include "MeshedSurfaces.H"
|
|
|
|
#include "extrudedMesh.H"
|
|
#include "extrudeModel.H"
|
|
|
|
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 wordList& names)
|
|
{
|
|
label n = 0;
|
|
|
|
forAll(names, i)
|
|
{
|
|
const polyPatch& pp = patches[findPatchID(patches, names[i])];
|
|
|
|
n += pp.size();
|
|
}
|
|
labelList faceLabels(n);
|
|
n = 0;
|
|
forAll(names, i)
|
|
{
|
|
const polyPatch& pp = patches[findPatchID(patches, names[i])];
|
|
|
|
forAll(pp, j)
|
|
{
|
|
faceLabels[n++] = pp.start()+j;
|
|
}
|
|
}
|
|
|
|
return faceLabels;
|
|
}
|
|
|
|
|
|
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"
|
|
|
|
// 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
|
|
(
|
|
IOobject
|
|
(
|
|
"extrudeProperties",
|
|
runTimeExtruded.constant(),
|
|
regionDir,
|
|
runTimeExtruded,
|
|
IOobject::MUST_READ
|
|
)
|
|
);
|
|
|
|
// Point generator
|
|
autoPtr<extrudeModel> model(extrudeModel::New(dict));
|
|
|
|
// Whether to flip normals
|
|
const Switch flipNormals(dict.lookup("flipNormals"));
|
|
|
|
// What to extrude
|
|
const ExtrudeMode mode = ExtrudeModeNames.read
|
|
(
|
|
dict.lookup("constructFrom")
|
|
);
|
|
|
|
|
|
// Generated mesh (one of either)
|
|
autoPtr<fvMesh> meshFromMesh;
|
|
autoPtr<polyMesh> 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)
|
|
{
|
|
if (flipNormals)
|
|
{
|
|
FatalErrorIn(args.executable())
|
|
<< "Flipping normals not supported for extrusions from patch."
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
fileName sourceCasePath(dict.lookup("sourceCase"));
|
|
sourceCasePath.expand();
|
|
fileName sourceRootDir = sourceCasePath.path();
|
|
fileName sourceCaseDir = sourceCasePath.name();
|
|
wordList sourcePatches;
|
|
dict.lookup("sourcePatches") >> sourcePatches;
|
|
|
|
if (sourcePatches.size() == 1)
|
|
{
|
|
frontPatchName = sourcePatches[0];
|
|
}
|
|
|
|
Info<< "Extruding patches " << sourcePatches
|
|
<< " on mesh " << sourceCasePath << nl
|
|
<< endl;
|
|
|
|
Time runTime
|
|
(
|
|
Time::controlDictName,
|
|
sourceRootDir,
|
|
sourceCaseDir
|
|
);
|
|
#include "createMesh.H"
|
|
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
|
|
|
|
// Topo change container. Either copy an existing mesh or start
|
|
// with empty storage (number of patches only needed for checking)
|
|
autoPtr<polyTopoChange> meshMod
|
|
(
|
|
(
|
|
mode == MESH
|
|
? new polyTopoChange(mesh)
|
|
: new polyTopoChange(patches.size())
|
|
)
|
|
);
|
|
|
|
// Extrusion engine. Either adding to existing mesh or
|
|
// creating separate mesh.
|
|
addPatchCellLayer layerExtrude(mesh, (mode == MESH));
|
|
|
|
indirectPrimitivePatch extrudePatch
|
|
(
|
|
IndirectList<face>
|
|
(
|
|
mesh.faces(),
|
|
patchFaces(patches, sourcePatches)
|
|
),
|
|
mesh.points()
|
|
);
|
|
|
|
// Only used for addPatchCellLayer into new mesh
|
|
labelList exposedPatchIDs;
|
|
if (mode == PATCH)
|
|
{
|
|
dict.lookup("exposedPatchName") >> backPatchName;
|
|
exposedPatchIDs.setSize
|
|
(
|
|
extrudePatch.size(),
|
|
findPatchID(patches, backPatchName)
|
|
);
|
|
}
|
|
|
|
|
|
pointField layer0Points(extrudePatch.nPoints());
|
|
pointField displacement(extrudePatch.nPoints());
|
|
forAll(displacement, pointI)
|
|
{
|
|
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].first();
|
|
frontPatchFaces[i] = layerFaces[i].last();
|
|
}
|
|
|
|
// 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
|
|
{
|
|
// Read from surface
|
|
fileName surfName(dict.lookup("surface"));
|
|
|
|
Info<< "Extruding surfaceMesh read from file " << surfName << nl
|
|
<< endl;
|
|
|
|
MeshedSurface<face> fMesh(surfName);
|
|
|
|
if (flipNormals)
|
|
{
|
|
Info<< "Flipping faces." << nl << endl;
|
|
faceList& faces = const_cast<faceList&>(fMesh.faces());
|
|
forAll(faces, i)
|
|
{
|
|
faces[i] = fMesh[i].reverseFace();
|
|
}
|
|
}
|
|
|
|
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(),
|
|
wordList(1, frontPatchName)
|
|
);
|
|
backPatchName = "otherSide";
|
|
backPatchFaces = patchFaces
|
|
(
|
|
meshFromSurface().boundaryMesh(),
|
|
wordList(1, backPatchName)
|
|
);
|
|
}
|
|
|
|
|
|
polyMesh& mesh =
|
|
(
|
|
meshFromMesh.valid()
|
|
? meshFromMesh()
|
|
: meshFromSurface()
|
|
);
|
|
|
|
|
|
const boundBox& bb = mesh.globalData().bb();
|
|
const vector span = bb.span();
|
|
const scalar mergeDim = 1E-4 * bb.minDim();
|
|
|
|
Info<< "Mesh bounding box : " << bb << nl
|
|
<< " with span : " << span << nl
|
|
<< "Merge distance : " << mergeDim << nl
|
|
<< endl;
|
|
|
|
|
|
// Collapse edges
|
|
// ~~~~~~~~~~~~~~
|
|
|
|
{
|
|
Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl;
|
|
|
|
// Edge collapsing engine
|
|
edgeCollapser collapser(mesh);
|
|
|
|
const edgeList& edges = mesh.edges();
|
|
const pointField& points = mesh.points();
|
|
|
|
forAll(edges, edgeI)
|
|
{
|
|
const edge& e = edges[edgeI];
|
|
|
|
scalar d = e.mag(points);
|
|
|
|
if (d < mergeDim)
|
|
{
|
|
Info<< "Merging edge " << e << " since length " << d
|
|
<< " << " << mergeDim << nl;
|
|
|
|
// Collapse edge to e[0]
|
|
collapser.collapseEdge(edgeI, e[0]);
|
|
}
|
|
}
|
|
|
|
// Topo change container
|
|
polyTopoChange meshMod(mesh);
|
|
// Put all modifications into meshMod
|
|
bool anyChange = collapser.setRefinement(meshMod);
|
|
|
|
if (anyChange)
|
|
{
|
|
// Construct new mesh from polyTopoChange.
|
|
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
|
|
|
|
// Update fields
|
|
mesh.updateMesh(map);
|
|
|
|
// Update stored data
|
|
updateFaceLabels(map(), frontPatchFaces);
|
|
updateFaceLabels(map(), backPatchFaces);
|
|
|
|
// Move mesh (if inflation used)
|
|
if (map().hasMotionPoints())
|
|
{
|
|
mesh.movePoints(map().preMotionPoints());
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Merging front and back patch faces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
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 "
|
|
<< 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);
|
|
|
|
const word cutZoneName("originalCutFaceZone");
|
|
|
|
List<faceZone*> fz
|
|
(
|
|
1,
|
|
new faceZone
|
|
(
|
|
cutZoneName,
|
|
frontPatchFaces,
|
|
boolList(frontPatchFaces.size(), false),
|
|
0,
|
|
mesh.faceZones()
|
|
)
|
|
);
|
|
|
|
mesh.addZones(List<pointZone*>(0), fz, List<cellZone*>(0));
|
|
|
|
// Add the perfect interface mesh modifier
|
|
perfectInterface perfectStitcher
|
|
(
|
|
"couple",
|
|
0,
|
|
stitcher,
|
|
cutZoneName,
|
|
word::null, // dummy patch name
|
|
word::null // dummy patch name
|
|
);
|
|
|
|
// Topo change container
|
|
polyTopoChange meshMod(mesh);
|
|
|
|
perfectStitcher.setRefinement
|
|
(
|
|
indirectPrimitivePatch
|
|
(
|
|
IndirectList<face>
|
|
(
|
|
mesh.faces(),
|
|
frontPatchFaces
|
|
),
|
|
mesh.points()
|
|
),
|
|
indirectPrimitivePatch
|
|
(
|
|
IndirectList<face>
|
|
(
|
|
mesh.faces(),
|
|
backPatchFaces
|
|
),
|
|
mesh.points()
|
|
),
|
|
meshMod
|
|
);
|
|
|
|
// Construct new mesh from polyTopoChange.
|
|
autoPtr<mapPolyMesh> 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"
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
Info<< "End\n" << endl;
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|