2949 lines
80 KiB
C
2949 lines
80 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 1991-2007 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
|
|
All to do with adding cell layers
|
|
|
|
\*----------------------------------------------------------------------------*/
|
|
|
|
#include "ListOps.H"
|
|
#include "autoHexMeshDriver.H"
|
|
#include "fvMesh.H"
|
|
#include "Time.H"
|
|
#include "combineFaces.H"
|
|
#include "removePoints.H"
|
|
#include "pointFields.H"
|
|
#include "motionSmoother.H"
|
|
#include "mathematicalConstants.H"
|
|
#include "pointSet.H"
|
|
#include "faceSet.H"
|
|
#include "cellSet.H"
|
|
#include "polyTopoChange.H"
|
|
#include "mapPolyMesh.H"
|
|
#include "addPatchCellLayer.H"
|
|
#include "mapDistributePolyMesh.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
const Foam::scalar Foam::autoHexMeshDriver::defaultConcaveAngle = 90;
|
|
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
Foam::label Foam::autoHexMeshDriver::mergePatchFacesUndo
|
|
(
|
|
const scalar minCos,
|
|
const scalar concaveCos,
|
|
const dictionary& motionDict
|
|
)
|
|
{
|
|
// Patch face merging engine
|
|
combineFaces faceCombiner(mesh_, true);
|
|
|
|
// Pick up all candidate cells on boundary
|
|
labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
|
|
|
|
{
|
|
labelList patchIDs(meshRefinement::addedPatches(globalToPatch_));
|
|
|
|
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
|
|
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchI = patchIDs[i];
|
|
|
|
const polyPatch& patch = patches[patchI];
|
|
|
|
if (!patch.coupled())
|
|
{
|
|
forAll(patch, i)
|
|
{
|
|
boundaryCells.insert(mesh_.faceOwner()[patch.start()+i]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Get all sets of faces that can be merged
|
|
labelListList allFaceSets
|
|
(
|
|
faceCombiner.getMergeSets
|
|
(
|
|
minCos,
|
|
concaveCos,
|
|
boundaryCells
|
|
)
|
|
);
|
|
|
|
label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
|
|
|
|
Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
|
|
|
|
if (nFaceSets > 0)
|
|
{
|
|
if (debug_)
|
|
{
|
|
faceSet allSets(mesh_, "allFaceSets", allFaceSets.size());
|
|
forAll(allFaceSets, setI)
|
|
{
|
|
forAll(allFaceSets[setI], i)
|
|
{
|
|
allSets.insert(allFaceSets[setI][i]);
|
|
}
|
|
}
|
|
Pout<< "Writing all faces to be merged to set "
|
|
<< allSets.objectPath() << endl;
|
|
allSets.write();
|
|
}
|
|
|
|
|
|
// Topology changes container
|
|
polyTopoChange meshMod(mesh_);
|
|
|
|
// Merge all faces of a set into the first face of the set.
|
|
faceCombiner.setRefinement(allFaceSets, meshMod);
|
|
|
|
// Experimental: store data for all the points that have been deleted
|
|
meshRefinerPtr_().storeData
|
|
(
|
|
faceCombiner.savedPointLabels(), // points to store
|
|
labelList(0), // faces to store
|
|
labelList(0) // cells to store
|
|
);
|
|
|
|
// Change the mesh (no inflation)
|
|
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
|
|
|
|
// Update fields
|
|
mesh_.updateMesh(map);
|
|
|
|
// Move mesh (since morphing does not do this)
|
|
if (map().hasMotionPoints())
|
|
{
|
|
mesh_.movePoints(map().preMotionPoints());
|
|
}
|
|
|
|
faceCombiner.updateMesh(map);
|
|
|
|
meshRefinerPtr_().updateMesh(map, labelList(0));
|
|
|
|
|
|
|
|
for (label iteration = 0; iteration < 100; iteration++)
|
|
{
|
|
Info<< nl
|
|
<< "Undo iteration " << iteration << nl
|
|
<< "----------------" << endl;
|
|
|
|
|
|
// Check mesh for errors
|
|
// ~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
faceSet errorFaces
|
|
(
|
|
mesh_,
|
|
"errorFaces",
|
|
mesh_.nFaces()-mesh_.nInternalFaces()
|
|
);
|
|
bool hasErrors = motionSmoother::checkMesh
|
|
(
|
|
false, // report
|
|
mesh_,
|
|
motionDict,
|
|
errorFaces
|
|
);
|
|
//if (checkEdgeConnectivity)
|
|
//{
|
|
// Info<< "Checking edge-face connectivity (duplicate faces"
|
|
// << " or non-consecutive shared vertices)" << endl;
|
|
//
|
|
// label nOldSize = errorFaces.size();
|
|
//
|
|
// hasErrors =
|
|
// mesh_.checkFaceFaces
|
|
// (
|
|
// false,
|
|
// &errorFaces
|
|
// )
|
|
// || hasErrors;
|
|
//
|
|
// Info<< "Detected additional "
|
|
// << returnReduce(errorFaces.size()-nOldSize, sumOp<label>())
|
|
// << " faces with illegal face-face connectivity"
|
|
// << endl;
|
|
//}
|
|
|
|
if (!hasErrors)
|
|
{
|
|
break;
|
|
}
|
|
|
|
|
|
if (debug_)
|
|
{
|
|
Pout<< "Writing all faces in error to faceSet "
|
|
<< errorFaces.objectPath() << nl << endl;
|
|
errorFaces.write();
|
|
}
|
|
|
|
|
|
// Check any master cells for using any of the error faces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
DynamicList<label> mastersToRestore(allFaceSets.size());
|
|
|
|
forAll(allFaceSets, setI)
|
|
{
|
|
label masterFaceI = faceCombiner.masterFace()[setI];
|
|
|
|
if (masterFaceI != -1)
|
|
{
|
|
label masterCellII = mesh_.faceOwner()[masterFaceI];
|
|
|
|
const cell& cFaces = mesh_.cells()[masterCellII];
|
|
|
|
forAll(cFaces, i)
|
|
{
|
|
if (errorFaces.found(cFaces[i]))
|
|
{
|
|
mastersToRestore.append(masterFaceI);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
mastersToRestore.shrink();
|
|
|
|
label nRestore = returnReduce
|
|
(
|
|
mastersToRestore.size(),
|
|
sumOp<label>()
|
|
);
|
|
|
|
Info<< "Masters that need to be restored:"
|
|
<< nRestore << endl;
|
|
|
|
if (debug_)
|
|
{
|
|
faceSet restoreSet(mesh_, "mastersToRestore", mastersToRestore);
|
|
Pout<< "Writing all " << mastersToRestore.size()
|
|
<< " masterfaces to be restored to set "
|
|
<< restoreSet.objectPath() << endl;
|
|
restoreSet.write();
|
|
}
|
|
|
|
|
|
if (nRestore == 0)
|
|
{
|
|
break;
|
|
}
|
|
|
|
|
|
// Undo
|
|
// ~~~~
|
|
|
|
// Topology changes container
|
|
polyTopoChange meshMod(mesh_);
|
|
|
|
// Merge all faces of a set into the first face of the set.
|
|
// Experimental:mark all points/faces/cells that have been restored.
|
|
Map<label> restoredPoints(0);
|
|
Map<label> restoredFaces(0);
|
|
Map<label> restoredCells(0);
|
|
|
|
faceCombiner.setUnrefinement
|
|
(
|
|
mastersToRestore,
|
|
meshMod,
|
|
restoredPoints,
|
|
restoredFaces,
|
|
restoredCells
|
|
);
|
|
|
|
// Change the mesh (no inflation)
|
|
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
|
|
|
|
// Update fields
|
|
mesh_.updateMesh(map);
|
|
|
|
// Move mesh (since morphing does not do this)
|
|
if (map().hasMotionPoints())
|
|
{
|
|
mesh_.movePoints(map().preMotionPoints());
|
|
}
|
|
|
|
faceCombiner.updateMesh(map);
|
|
|
|
// Renumber restore maps
|
|
inplaceMapKey(map().reversePointMap(), restoredPoints);
|
|
inplaceMapKey(map().reverseFaceMap(), restoredFaces);
|
|
inplaceMapKey(map().reverseCellMap(), restoredCells);
|
|
|
|
// Experimental:restore all points/face/cells in maps
|
|
meshRefinerPtr_().updateMesh
|
|
(
|
|
map,
|
|
labelList(0), // changedFaces
|
|
restoredPoints,
|
|
restoredFaces,
|
|
restoredCells
|
|
);
|
|
|
|
Info<< endl;
|
|
}
|
|
|
|
if (debug_)
|
|
{
|
|
Pout<< "Writing merged-faces mesh to time "
|
|
<< mesh_.time().timeName() << nl << endl;
|
|
mesh_.write();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
Info<< "No faces merged ..." << endl;
|
|
}
|
|
|
|
return nFaceSets;
|
|
}
|
|
|
|
|
|
// Remove points. pointCanBeDeleted is parallel synchronised.
|
|
Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::doRemovePoints
|
|
(
|
|
removePoints& pointRemover,
|
|
const boolList& pointCanBeDeleted
|
|
)
|
|
{
|
|
// Topology changes container
|
|
polyTopoChange meshMod(mesh_);
|
|
|
|
pointRemover.setRefinement(pointCanBeDeleted, meshMod);
|
|
|
|
// Change the mesh (no inflation)
|
|
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
|
|
|
|
// Update fields
|
|
mesh_.updateMesh(map);
|
|
|
|
// Move mesh (since morphing does not do this)
|
|
if (map().hasMotionPoints())
|
|
{
|
|
mesh_.movePoints(map().preMotionPoints());
|
|
}
|
|
|
|
pointRemover.updateMesh(map);
|
|
meshRefinerPtr_().updateMesh(map, labelList(0));
|
|
|
|
return map;
|
|
}
|
|
|
|
|
|
// Restore faces (which contain removed points)
|
|
Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::doRestorePoints
|
|
(
|
|
removePoints& pointRemover,
|
|
const labelList& facesToRestore
|
|
)
|
|
{
|
|
// Topology changes container
|
|
polyTopoChange meshMod(mesh_);
|
|
|
|
// Determine sets of points and faces to restore
|
|
labelList localFaces, localPoints;
|
|
pointRemover.getUnrefimentSet
|
|
(
|
|
facesToRestore,
|
|
localFaces,
|
|
localPoints
|
|
);
|
|
|
|
// Undo the changes on the faces that are in error.
|
|
pointRemover.setUnrefinement
|
|
(
|
|
localFaces,
|
|
localPoints,
|
|
meshMod
|
|
);
|
|
|
|
// Change the mesh (no inflation)
|
|
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
|
|
|
|
// Update fields
|
|
mesh_.updateMesh(map);
|
|
|
|
// Move mesh (since morphing does not do this)
|
|
if (map().hasMotionPoints())
|
|
{
|
|
mesh_.movePoints(map().preMotionPoints());
|
|
}
|
|
|
|
pointRemover.updateMesh(map);
|
|
meshRefinerPtr_().updateMesh(map, labelList(0));
|
|
|
|
return map;
|
|
}
|
|
|
|
|
|
// Collect all faces that are both in candidateFaces and in the set.
|
|
// If coupled face also collects the coupled face.
|
|
Foam::labelList Foam::autoHexMeshDriver::collectFaces
|
|
(
|
|
const labelList& candidateFaces,
|
|
const labelHashSet& set
|
|
) const
|
|
{
|
|
// Has face been selected?
|
|
boolList selected(mesh_.nFaces(), false);
|
|
|
|
forAll(candidateFaces, i)
|
|
{
|
|
label faceI = candidateFaces[i];
|
|
|
|
if (set.found(faceI))
|
|
{
|
|
selected[faceI] = true;
|
|
}
|
|
}
|
|
syncTools::syncFaceList
|
|
(
|
|
mesh_,
|
|
selected,
|
|
orEqOp<bool>(), // combine operator
|
|
false // separation
|
|
);
|
|
|
|
labelList selectedFaces(findIndices(selected, true));
|
|
|
|
return selectedFaces;
|
|
}
|
|
|
|
|
|
// Pick up faces of cells of faces in set.
|
|
Foam::labelList Foam::autoHexMeshDriver::growFaceCellFace
|
|
(
|
|
const labelHashSet& set
|
|
) const
|
|
{
|
|
boolList selected(mesh_.nFaces(), false);
|
|
|
|
forAllConstIter(faceSet, set, iter)
|
|
{
|
|
label faceI = iter.key();
|
|
|
|
label own = mesh_.faceOwner()[faceI];
|
|
|
|
const cell& ownFaces = mesh_.cells()[own];
|
|
forAll(ownFaces, ownFaceI)
|
|
{
|
|
selected[ownFaces[ownFaceI]] = true;
|
|
}
|
|
|
|
if (mesh_.isInternalFace(faceI))
|
|
{
|
|
label nbr = mesh_.faceNeighbour()[faceI];
|
|
|
|
const cell& nbrFaces = mesh_.cells()[nbr];
|
|
forAll(nbrFaces, nbrFaceI)
|
|
{
|
|
selected[nbrFaces[nbrFaceI]] = true;
|
|
}
|
|
}
|
|
}
|
|
syncTools::syncFaceList
|
|
(
|
|
mesh_,
|
|
selected,
|
|
orEqOp<bool>(), // combine operator
|
|
false // separation
|
|
);
|
|
return findIndices(selected, true);
|
|
}
|
|
|
|
|
|
// Remove points not used by any face or points used by only two faces where
|
|
// the edges are in line
|
|
Foam::label Foam::autoHexMeshDriver::mergeEdgesUndo
|
|
(
|
|
const scalar minCos,
|
|
const dictionary& motionDict
|
|
)
|
|
{
|
|
Info<< nl
|
|
<< "Merging all points on surface that" << nl
|
|
<< "- are used by only two boundary faces and" << nl
|
|
<< "- make an angle with a cosine of more than " << minCos
|
|
<< "." << nl << endl;
|
|
|
|
// Point removal analysis engine with undo
|
|
removePoints pointRemover(mesh_, true);
|
|
|
|
// Count usage of points
|
|
boolList pointCanBeDeleted;
|
|
label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
|
|
|
|
if (nRemove > 0)
|
|
{
|
|
Info<< "Removing " << nRemove
|
|
<< " straight edge points ..." << nl << endl;
|
|
|
|
// Remove points
|
|
// ~~~~~~~~~~~~~
|
|
|
|
doRemovePoints(pointRemover, pointCanBeDeleted);
|
|
|
|
|
|
for (label iteration = 0; iteration < 100; iteration++)
|
|
{
|
|
Info<< nl
|
|
<< "Undo iteration " << iteration << nl
|
|
<< "----------------" << endl;
|
|
|
|
|
|
// Check mesh for errors
|
|
// ~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
faceSet errorFaces
|
|
(
|
|
mesh_,
|
|
"errorFaces",
|
|
mesh_.nFaces()-mesh_.nInternalFaces()
|
|
);
|
|
bool hasErrors = motionSmoother::checkMesh
|
|
(
|
|
false, // report
|
|
mesh_,
|
|
motionDict,
|
|
errorFaces
|
|
);
|
|
//if (checkEdgeConnectivity)
|
|
//{
|
|
// Info<< "Checking edge-face connectivity (duplicate faces"
|
|
// << " or non-consecutive shared vertices)" << endl;
|
|
//
|
|
// label nOldSize = errorFaces.size();
|
|
//
|
|
// hasErrors =
|
|
// mesh_.checkFaceFaces
|
|
// (
|
|
// false,
|
|
// &errorFaces
|
|
// )
|
|
// || hasErrors;
|
|
//
|
|
// Info<< "Detected additional "
|
|
// << returnReduce(errorFaces.size()-nOldSize, sumOp<label>())
|
|
// << " faces with illegal face-face connectivity"
|
|
// << endl;
|
|
//}
|
|
|
|
if (!hasErrors)
|
|
{
|
|
break;
|
|
}
|
|
|
|
if (debug_)
|
|
{
|
|
Pout<< "**Writing all faces in error to faceSet "
|
|
<< errorFaces.objectPath() << nl << endl;
|
|
errorFaces.write();
|
|
}
|
|
|
|
labelList masterErrorFaces
|
|
(
|
|
collectFaces
|
|
(
|
|
pointRemover.savedFaceLabels(),
|
|
errorFaces
|
|
)
|
|
);
|
|
|
|
label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
|
|
|
|
Info<< "Detected " << n
|
|
<< " error faces on boundaries that have been merged."
|
|
<< " These will be restored to their original faces." << nl
|
|
<< endl;
|
|
|
|
if (n == 0)
|
|
{
|
|
if (hasErrors)
|
|
{
|
|
Info<< "Detected "
|
|
<< returnReduce(errorFaces.size(), sumOp<label>())
|
|
<< " error faces in mesh."
|
|
<< " Restoring neighbours of faces in error." << nl
|
|
<< endl;
|
|
|
|
labelList expandedErrorFaces
|
|
(
|
|
growFaceCellFace
|
|
(
|
|
errorFaces
|
|
)
|
|
);
|
|
|
|
doRestorePoints(pointRemover, expandedErrorFaces);
|
|
}
|
|
|
|
break;
|
|
}
|
|
|
|
doRestorePoints(pointRemover, masterErrorFaces);
|
|
}
|
|
|
|
if (debug_)
|
|
{
|
|
Pout<< "Writing merged-edges mesh to time "
|
|
<< mesh_.time().timeName() << nl << endl;
|
|
mesh_.write();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
Info<< "No straight edges simplified and no points removed ..." << endl;
|
|
}
|
|
|
|
return nRemove;
|
|
}
|
|
|
|
|
|
// For debugging: Dump displacement to .obj files
|
|
void Foam::autoHexMeshDriver::dumpDisplacement
|
|
(
|
|
const fileName& prefix,
|
|
const indirectPrimitivePatch& pp,
|
|
const vectorField& patchDisp,
|
|
const List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
OFstream dispStr(prefix + "_disp.obj");
|
|
Info<< "Writing all displacements to " << dispStr.name() << nl << endl;
|
|
|
|
label vertI = 0;
|
|
|
|
forAll(patchDisp, patchPointI)
|
|
{
|
|
const point& pt = pp.localPoints()[patchPointI];
|
|
|
|
meshTools::writeOBJ(dispStr, pt); vertI++;
|
|
meshTools::writeOBJ(dispStr, pt + patchDisp[patchPointI]); vertI++;
|
|
|
|
dispStr << "l " << vertI-1 << ' ' << vertI << nl;
|
|
}
|
|
|
|
|
|
OFstream illStr(prefix + "_illegal.obj");
|
|
Info<< "Writing invalid displacements to " << illStr.name() << nl << endl;
|
|
|
|
vertI = 0;
|
|
|
|
forAll(patchDisp, patchPointI)
|
|
{
|
|
if (extrudeStatus[patchPointI] != EXTRUDE)
|
|
{
|
|
const point& pt = pp.localPoints()[patchPointI];
|
|
|
|
meshTools::writeOBJ(illStr, pt); vertI++;
|
|
meshTools::writeOBJ(illStr, pt + patchDisp[patchPointI]); vertI++;
|
|
|
|
illStr << "l " << vertI-1 << ' ' << vertI << nl;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Check that primitivePatch is not multiply connected. Collect non-manifold
|
|
// points in pointSet.
|
|
void Foam::autoHexMeshDriver::checkManifold
|
|
(
|
|
const indirectPrimitivePatch& fp,
|
|
pointSet& nonManifoldPoints
|
|
)
|
|
{
|
|
// Check for non-manifold points (surface pinched at point)
|
|
fp.checkPointManifold(false, &nonManifoldPoints);
|
|
|
|
// Check for edge-faces (surface pinched at edge)
|
|
const labelListList& edgeFaces = fp.edgeFaces();
|
|
|
|
forAll(edgeFaces, edgeI)
|
|
{
|
|
const labelList& eFaces = edgeFaces[edgeI];
|
|
|
|
if (eFaces.size() > 2)
|
|
{
|
|
const edge& e = fp.edges()[edgeI];
|
|
|
|
nonManifoldPoints.insert(fp.meshPoints()[e[0]]);
|
|
nonManifoldPoints.insert(fp.meshPoints()[e[1]]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::autoHexMeshDriver::checkMeshManifold() const
|
|
{
|
|
Info<< nl << "Checking mesh manifoldness ..." << endl;
|
|
|
|
// Get all outside faces
|
|
labelList outsideFaces(mesh_.nFaces() - mesh_.nInternalFaces());
|
|
|
|
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
|
|
{
|
|
outsideFaces[faceI - mesh_.nInternalFaces()] = faceI;
|
|
}
|
|
|
|
pointSet nonManifoldPoints
|
|
(
|
|
mesh_,
|
|
"nonManifoldPoints",
|
|
mesh_.nPoints() / 100
|
|
);
|
|
|
|
// Build primitivePatch out of faces and check it for problems.
|
|
checkManifold
|
|
(
|
|
indirectPrimitivePatch
|
|
(
|
|
IndirectList<face>(mesh_.faces(), outsideFaces),
|
|
mesh_.points()
|
|
),
|
|
nonManifoldPoints
|
|
);
|
|
|
|
label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
|
|
|
|
if (nNonManif > 0)
|
|
{
|
|
Info<< "Outside of mesh is multiply connected across edges or"
|
|
<< " points." << nl
|
|
<< "This is not a fatal error but might cause some unexpected"
|
|
<< " behaviour." << nl
|
|
<< "Writing " << nNonManif
|
|
<< " points where this happens to pointSet "
|
|
<< nonManifoldPoints.name() << endl;
|
|
|
|
nonManifoldPoints.write();
|
|
}
|
|
Info<< endl;
|
|
}
|
|
|
|
|
|
|
|
// Unset extrusion on point. Returns true if anything unset.
|
|
bool Foam::autoHexMeshDriver::unmarkExtrusion
|
|
(
|
|
const label patchPointI,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
if (extrudeStatus[patchPointI] == EXTRUDE)
|
|
{
|
|
extrudeStatus[patchPointI] = NOEXTRUDE;
|
|
patchNLayers[patchPointI] = 0;
|
|
patchDisp[patchPointI] = vector::zero;
|
|
return true;
|
|
}
|
|
else if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
|
|
{
|
|
extrudeStatus[patchPointI] = NOEXTRUDE;
|
|
patchNLayers[patchPointI] = 0;
|
|
patchDisp[patchPointI] = vector::zero;
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
|
|
|
|
// Unset extrusion on face. Returns true if anything unset.
|
|
bool Foam::autoHexMeshDriver::unmarkExtrusion
|
|
(
|
|
const face& localFace,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
bool unextruded = false;
|
|
|
|
forAll(localFace, fp)
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
localFace[fp],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
unextruded = true;
|
|
}
|
|
}
|
|
return unextruded;
|
|
}
|
|
|
|
|
|
// No extrusion at non-manifold points.
|
|
void Foam::autoHexMeshDriver::handleNonManifolds
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const labelList& meshEdges,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
Info<< nl << "Handling non-manifold points ..." << endl;
|
|
|
|
// Detect non-manifold points
|
|
Info<< nl << "Checking patch manifoldness ..." << endl;
|
|
|
|
pointSet nonManifoldPoints(mesh_, "nonManifoldPoints", pp.nPoints());
|
|
|
|
// 1. Local check
|
|
checkManifold(pp, nonManifoldPoints);
|
|
|
|
label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
|
|
|
|
Info<< "Outside of local patch is multiply connected across edges or"
|
|
<< " points at " << nNonManif << " points." << endl;
|
|
|
|
|
|
const labelList& meshPoints = pp.meshPoints();
|
|
|
|
|
|
// 2. Parallel check
|
|
// For now disable extrusion at any shared edges.
|
|
const labelHashSet sharedEdgeSet(mesh_.globalData().sharedEdgeLabels());
|
|
|
|
forAll(pp.edges(), edgeI)
|
|
{
|
|
if (sharedEdgeSet.found(meshEdges[edgeI]))
|
|
{
|
|
const edge& e = mesh_.edges()[meshEdges[edgeI]];
|
|
|
|
Pout<< "Disabling extrusion at edge "
|
|
<< mesh_.points()[e[0]] << mesh_.points()[e[1]]
|
|
<< " since it is non-manifold across coupled patches."
|
|
<< endl;
|
|
|
|
nonManifoldPoints.insert(e[0]);
|
|
nonManifoldPoints.insert(e[1]);
|
|
}
|
|
}
|
|
|
|
// 3b. extrusion can produce multiple faces between 2 cells
|
|
// across processor boundary
|
|
// This occurs when a coupled face shares more than 1 edge with a
|
|
// non-coupled boundary face.
|
|
// This is now correctly handled by addPatchCellLayer in that it
|
|
// extrudes a single face from the stringed up edges.
|
|
|
|
|
|
nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
|
|
|
|
if (nNonManif > 0)
|
|
{
|
|
Info<< "Outside of patches is multiply connected across edges or"
|
|
<< " points at " << nNonManif << " points." << nl
|
|
<< "Writing " << nNonManif
|
|
<< " points where this happens to pointSet "
|
|
<< nonManifoldPoints.name() << nl
|
|
<< "and setting layer thickness to zero on these points."
|
|
<< endl;
|
|
|
|
nonManifoldPoints.write();
|
|
|
|
forAll(meshPoints, patchPointI)
|
|
{
|
|
if (nonManifoldPoints.found(meshPoints[patchPointI]))
|
|
{
|
|
unmarkExtrusion
|
|
(
|
|
patchPointI,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Set displacement to zero for all " << nNonManif
|
|
<< " non-manifold points" << endl;
|
|
}
|
|
|
|
|
|
// Parallel feature edge detection. Assumes non-manifold edges already handled.
|
|
void Foam::autoHexMeshDriver::handleFeatureAngle
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const labelList& meshEdges,
|
|
const scalar minCos,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
Info<< nl << "Handling feature edges ..." << endl;
|
|
|
|
if (minCos < 1-SMALL)
|
|
{
|
|
// Normal component of normals of connected faces.
|
|
vectorField edgeNormal(mesh_.nEdges(), wallPoint::greatPoint);
|
|
|
|
const labelListList& edgeFaces = pp.edgeFaces();
|
|
|
|
forAll(edgeFaces, edgeI)
|
|
{
|
|
const labelList& eFaces = pp.edgeFaces()[edgeI];
|
|
|
|
label meshEdgeI = meshEdges[edgeI];
|
|
|
|
forAll(eFaces, i)
|
|
{
|
|
nomalsCombine()
|
|
(
|
|
edgeNormal[meshEdgeI],
|
|
pp.faceNormals()[eFaces[i]]
|
|
);
|
|
}
|
|
}
|
|
|
|
syncTools::syncEdgeList
|
|
(
|
|
mesh_,
|
|
edgeNormal,
|
|
nomalsCombine(),
|
|
wallPoint::greatPoint, // null value
|
|
false // no separation
|
|
);
|
|
|
|
label vertI = 0;
|
|
autoPtr<OFstream> str;
|
|
if (debug_)
|
|
{
|
|
str.reset(new OFstream(mesh_.time().path()/"featureEdges.obj"));
|
|
Info<< "Writing feature edges to " << str().name() << endl;
|
|
}
|
|
|
|
label nFeats = 0;
|
|
|
|
// Now on coupled edges the edgeNormal will have been truncated and
|
|
// only be still be the old value where two faces have the same normal
|
|
forAll(edgeFaces, edgeI)
|
|
{
|
|
const labelList& eFaces = pp.edgeFaces()[edgeI];
|
|
|
|
label meshEdgeI = meshEdges[edgeI];
|
|
|
|
const vector& n = edgeNormal[meshEdgeI];
|
|
|
|
if (n != wallPoint::greatPoint)
|
|
{
|
|
scalar cos = n & pp.faceNormals()[eFaces[0]];
|
|
|
|
if (cos < minCos)
|
|
{
|
|
const edge& e = pp.edges()[edgeI];
|
|
|
|
unmarkExtrusion
|
|
(
|
|
e[0],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
unmarkExtrusion
|
|
(
|
|
e[1],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
nFeats++;
|
|
|
|
if (str.valid())
|
|
{
|
|
meshTools::writeOBJ(str(), pp.localPoints()[e[0]]);
|
|
vertI++;
|
|
meshTools::writeOBJ(str(), pp.localPoints()[e[1]]);
|
|
vertI++;
|
|
str()<< "l " << vertI-1 << ' ' << vertI << nl;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Set displacement to zero for points on "
|
|
<< returnReduce(nFeats, sumOp<label>())
|
|
<< " feature edges" << endl;
|
|
}
|
|
}
|
|
|
|
|
|
// No extrusion on cells with warped faces. Calculates the thickness of the
|
|
// layer and compares it to the space the warped face takes up. Disables
|
|
// extrusion if layer thickness is more than faceRatio of the thickness of
|
|
// the face.
|
|
void Foam::autoHexMeshDriver::handleWarpedFaces
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const scalar faceRatio,
|
|
const scalar edge0Len,
|
|
const labelList& cellLevel,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
Info<< nl << "Handling cells with warped patch faces ..." << nl;
|
|
|
|
const pointField& points = mesh_.points();
|
|
|
|
label nWarpedFaces = 0;
|
|
|
|
forAll(pp, i)
|
|
{
|
|
const face& f = pp[i];
|
|
|
|
if (f.size() > 3)
|
|
{
|
|
label faceI = pp.addressing()[i];
|
|
|
|
label ownLevel = cellLevel[mesh_.faceOwner()[faceI]];
|
|
scalar edgeLen = edge0Len/(1<<ownLevel);
|
|
|
|
// Normal distance to face centre plane
|
|
const point& fc = mesh_.faceCentres()[faceI];
|
|
const vector& fn = pp.faceNormals()[i];
|
|
|
|
scalarField vProj(f.size());
|
|
|
|
forAll(f, fp)
|
|
{
|
|
vector n = points[f[fp]] - fc;
|
|
vProj[fp] = (n & fn);
|
|
}
|
|
|
|
// Get normal 'span' of face
|
|
scalar minVal = min(vProj);
|
|
scalar maxVal = max(vProj);
|
|
|
|
if ((maxVal - minVal) > faceRatio * edgeLen)
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
pp.localFaces()[i],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nWarpedFaces++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Set displacement to zero on "
|
|
<< returnReduce(nWarpedFaces, sumOp<label>())
|
|
<< " warped faces since layer would be > " << faceRatio
|
|
<< " of the size of the bounding box." << endl;
|
|
}
|
|
|
|
|
|
|
|
//// No extrusion on cells with multiple patch faces. There ususally is a reason
|
|
//// why combinePatchFaces hasn't succeeded.
|
|
//void Foam::autoHexMeshDriver::handleMultiplePatchFaces
|
|
//(
|
|
// const indirectPrimitivePatch& pp,
|
|
// pointField& patchDisp,
|
|
// labelList& patchNLayers,
|
|
// List<extrudeMode>& extrudeStatus
|
|
//) const
|
|
//{
|
|
// Info<< nl << "Handling cells with multiple patch faces ..." << nl;
|
|
//
|
|
// const labelListList& pointFaces = pp.pointFaces();
|
|
//
|
|
// // Cells that should not get an extrusion layer
|
|
// cellSet multiPatchCells(mesh_, "multiPatchCells", pp.size());
|
|
//
|
|
// // Detect points that use multiple faces on same cell.
|
|
// forAll(pointFaces, patchPointI)
|
|
// {
|
|
// const labelList& pFaces = pointFaces[patchPointI];
|
|
//
|
|
// labelHashSet pointCells(pFaces.size());
|
|
//
|
|
// forAll(pFaces, i)
|
|
// {
|
|
// label cellI = mesh_.faceOwner()[pp.addressing()[pFaces[i]]];
|
|
//
|
|
// if (!pointCells.insert(cellI))
|
|
// {
|
|
// // Second or more occurrence of cell so cell has two or more
|
|
// // pp faces connected to this point.
|
|
// multiPatchCells.insert(cellI);
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// label nMultiPatchCells = returnReduce
|
|
// (
|
|
// multiPatchCells.size(),
|
|
// sumOp<label>()
|
|
// );
|
|
//
|
|
// Info<< "Detected " << nMultiPatchCells
|
|
// << " cells with multiple (connected) patch faces." << endl;
|
|
//
|
|
// label nChanged = 0;
|
|
//
|
|
// if (nMultiPatchCells > 0)
|
|
// {
|
|
// Info<< "Writing " << nMultiPatchCells
|
|
// << " cells with multiple (connected) patch faces to cellSet "
|
|
// << multiPatchCells.objectPath() << endl;
|
|
// multiPatchCells.write();
|
|
//
|
|
//
|
|
// // Go through all points and remove extrusion on any cell in
|
|
// // multiPatchCells
|
|
// // (has to be done in separate loop since having one point on
|
|
// // multipatches has to reset extrusion on all points of cell)
|
|
//
|
|
// forAll(pointFaces, patchPointI)
|
|
// {
|
|
// if (extrudeStatus[patchPointI] != NOEXTRUDE)
|
|
// {
|
|
// const labelList& pFaces = pointFaces[patchPointI];
|
|
//
|
|
// forAll(pFaces, i)
|
|
// {
|
|
// label cellI =
|
|
// mesh_.faceOwner()[pp.addressing()[pFaces[i]]];
|
|
//
|
|
// if (multiPatchCells.found(cellI))
|
|
// {
|
|
// if
|
|
// (
|
|
// unmarkExtrusion
|
|
// (
|
|
// patchPointI,
|
|
// patchDisp,
|
|
// patchNLayers,
|
|
// extrudeStatus
|
|
// )
|
|
// )
|
|
// {
|
|
// nChanged++;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// reduce(nChanged, sumOp<label>());
|
|
// }
|
|
//
|
|
// Info<< "Prevented extrusion on " << nChanged
|
|
// << " points due to multiple patch faces." << nl << endl;
|
|
//}
|
|
|
|
|
|
// No extrusion on faces with differing number of layers for points
|
|
void Foam::autoHexMeshDriver::setNumLayers
|
|
(
|
|
const labelList& patchIDs,
|
|
const indirectPrimitivePatch& pp,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
Info<< nl << "Handling points with inconsistent layer specification ..."
|
|
<< endl;
|
|
|
|
const labelList& nSurfLayers = surfaces().numLayers();
|
|
|
|
// Build map from patch to layers
|
|
Map<label> patchToNLayers(nSurfLayers.size());
|
|
forAll(nSurfLayers, region)
|
|
{
|
|
if (globalToPatch_[region] != -1)
|
|
{
|
|
patchToNLayers.insert(globalToPatch_[region], nSurfLayers[region]);
|
|
}
|
|
}
|
|
|
|
// Get for every point (really only nessecary on patch external points)
|
|
// the max and min of any patch faces using it.
|
|
labelList maxLayers(patchNLayers.size(), labelMin);
|
|
labelList minLayers(patchNLayers.size(), labelMax);
|
|
|
|
forAll(patchIDs, i)
|
|
{
|
|
label patchI = patchIDs[i];
|
|
|
|
const labelList& meshPoints = mesh_.boundaryMesh()[patchI].meshPoints();
|
|
|
|
label wantedLayers = patchToNLayers[patchI];
|
|
|
|
forAll(meshPoints, patchPointI)
|
|
{
|
|
label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
|
|
|
|
maxLayers[ppPointI] = max(wantedLayers, maxLayers[ppPointI]);
|
|
minLayers[ppPointI] = min(wantedLayers, minLayers[ppPointI]);
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh_,
|
|
pp.meshPoints(),
|
|
maxLayers,
|
|
maxEqOp<label>(),
|
|
labelMin, // null value
|
|
false // no separation
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh_,
|
|
pp.meshPoints(),
|
|
minLayers,
|
|
minEqOp<label>(),
|
|
labelMax, // null value
|
|
false // no separation
|
|
);
|
|
|
|
// Unmark any point with different min and max
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
label nConflicts = 0;
|
|
|
|
forAll(maxLayers, i)
|
|
{
|
|
if (maxLayers[i] == labelMin || minLayers[i] == labelMax)
|
|
{
|
|
FatalErrorIn("setNumLayers(..)")
|
|
<< "Patchpoint:" << i << " coord:" << pp.localPoints()[i]
|
|
<< " maxLayers:" << maxLayers
|
|
<< " minLayers:" << minLayers
|
|
<< abort(FatalError);
|
|
}
|
|
else if (maxLayers[i] == minLayers[i])
|
|
{
|
|
// Ok setting.
|
|
patchNLayers[i] = maxLayers[i];
|
|
}
|
|
else
|
|
{
|
|
// Inconsistent num layers between patch faces using point
|
|
//if
|
|
//(
|
|
// unmarkExtrusion
|
|
// (
|
|
// i,
|
|
// patchDisp,
|
|
// patchNLayers,
|
|
// extrudeStatus
|
|
// )
|
|
//)
|
|
//{
|
|
// nConflicts++;
|
|
//}
|
|
patchNLayers[i] = maxLayers[i];
|
|
}
|
|
}
|
|
|
|
reduce(nConflicts, sumOp<label>());
|
|
|
|
Info<< "Set displacement to zero for " << nConflicts
|
|
<< " points due to points being on multiple regions"
|
|
<< " with inconsistent nLayers specification." << endl;
|
|
}
|
|
|
|
|
|
// Grow no-extrusion layer
|
|
void Foam::autoHexMeshDriver::growNoExtrusion
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
|
|
|
|
List<extrudeMode> grownExtrudeStatus(extrudeStatus);
|
|
|
|
const faceList& localFaces = pp.localFaces();
|
|
|
|
label nGrown = 0;
|
|
|
|
forAll(localFaces, faceI)
|
|
{
|
|
const face& f = localFaces[faceI];
|
|
|
|
bool hasSqueeze = false;
|
|
forAll(f, fp)
|
|
{
|
|
if (extrudeStatus[f[fp]] == NOEXTRUDE)
|
|
{
|
|
hasSqueeze = true;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (hasSqueeze)
|
|
{
|
|
// Squeeze all points of face
|
|
forAll(f, fp)
|
|
{
|
|
if
|
|
(
|
|
extrudeStatus[f[fp]] == NOEXTRUDE
|
|
&& grownExtrudeStatus[f[fp]] != NOEXTRUDE
|
|
)
|
|
{
|
|
grownExtrudeStatus[f[fp]] = NOEXTRUDE;
|
|
nGrown++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
extrudeStatus.transfer(grownExtrudeStatus);
|
|
|
|
forAll(extrudeStatus, patchPointI)
|
|
{
|
|
if (extrudeStatus[patchPointI] == NOEXTRUDE)
|
|
{
|
|
patchDisp[patchPointI] = vector::zero;
|
|
patchNLayers[patchPointI] = 0;
|
|
}
|
|
}
|
|
|
|
reduce(nGrown, sumOp<label>());
|
|
|
|
Info<< "Set displacement to zero for an additional " << nGrown
|
|
<< " points." << endl;
|
|
}
|
|
|
|
|
|
|
|
// Calculate pointwise wanted and minimum thickness.
|
|
// thickness: wanted thickness per point
|
|
void Foam::autoHexMeshDriver::calculateLayerThickness
|
|
(
|
|
const scalar expansionRatio,
|
|
const scalar finalLayerRatio,
|
|
const scalar relMinThickness,
|
|
const indirectPrimitivePatch& pp,
|
|
const labelList& cellLevel,
|
|
const labelList& patchNLayers,
|
|
const scalar edge0Len,
|
|
scalarField& thickness,
|
|
scalarField& minThickness
|
|
) const
|
|
{
|
|
if (relMinThickness < 0 || relMinThickness > 2)
|
|
{
|
|
FatalErrorIn("calculateLayerThickness(..)")
|
|
<< "Thickness should be factor of local undistorted cell size."
|
|
<< " Valid values are [0..2]." << nl
|
|
<< " minThickness:" << relMinThickness
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
thickness.setSize(pp.nPoints());
|
|
minThickness.setSize(pp.nPoints());
|
|
|
|
// Per point the max cell level of connected cells
|
|
labelList maxPointLevel(pp.nPoints(), labelMin);
|
|
|
|
forAll(pp, i)
|
|
{
|
|
label ownLevel = cellLevel[mesh_.faceOwner()[pp.addressing()[i]]];
|
|
|
|
const face& f = pp.localFaces()[i];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh_,
|
|
pp.meshPoints(),
|
|
maxPointLevel,
|
|
maxEqOp<label>(),
|
|
labelMin, // null value
|
|
false // no separation
|
|
);
|
|
|
|
|
|
forAll(maxPointLevel, pointI)
|
|
{
|
|
// Find undistorted edge size for this level.
|
|
scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
|
|
|
|
// Calculate layer thickness based on expansion ratio
|
|
// and final layer height
|
|
if (expansionRatio == 1)
|
|
{
|
|
thickness[pointI] = finalLayerRatio*patchNLayers[pointI]*edgeLen;
|
|
minThickness[pointI] = relMinThickness*edgeLen;
|
|
}
|
|
else
|
|
{
|
|
scalar invExpansion = 1.0 / expansionRatio;
|
|
label nLay = patchNLayers[pointI];
|
|
thickness[pointI] = finalLayerRatio*edgeLen
|
|
* (1.0 - pow(invExpansion, nLay))
|
|
/ (1.0 - invExpansion);
|
|
minThickness[pointI] = relMinThickness*edgeLen;
|
|
}
|
|
}
|
|
|
|
Info<< "calculateLayerThickness : min:" << gMin(thickness)
|
|
<< " max:" << gMax(thickness) << endl;
|
|
}
|
|
|
|
|
|
// Synchronize displacement among coupled patches.
|
|
void Foam::autoHexMeshDriver::syncPatchDisplacement
|
|
(
|
|
const motionSmoother& meshMover,
|
|
const scalarField& minThickness,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
const labelList& meshPoints = meshMover.patch().meshPoints();
|
|
|
|
label nChangedTotal = 0;
|
|
|
|
while (true)
|
|
{
|
|
label nChanged = 0;
|
|
|
|
// Sync displacement (by taking min)
|
|
syncTools::syncPointList
|
|
(
|
|
mesh_,
|
|
meshPoints,
|
|
patchDisp,
|
|
minEqOp<vector>(),
|
|
wallPoint::greatPoint, // null value
|
|
false // no separation
|
|
);
|
|
|
|
// Unmark if displacement too small
|
|
forAll(patchDisp, i)
|
|
{
|
|
if (mag(patchDisp[i]) < minThickness[i])
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
i,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
|
|
labelList syncPatchNLayers(patchNLayers);
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh_,
|
|
meshPoints,
|
|
syncPatchNLayers,
|
|
minEqOp<label>(),
|
|
labelMax, // null value
|
|
false // no separation
|
|
);
|
|
|
|
// Reset if differs
|
|
forAll(syncPatchNLayers, i)
|
|
{
|
|
if (syncPatchNLayers[i] != patchNLayers[i])
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
i,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh_,
|
|
meshPoints,
|
|
syncPatchNLayers,
|
|
maxEqOp<label>(),
|
|
labelMin, // null value
|
|
false // no separation
|
|
);
|
|
|
|
// Reset if differs
|
|
forAll(syncPatchNLayers, i)
|
|
{
|
|
if (syncPatchNLayers[i] != patchNLayers[i])
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
i,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
nChangedTotal += nChanged;
|
|
|
|
if (!returnReduce(nChanged, sumOp<label>()))
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
Info<< "Prevented extrusion on "
|
|
<< returnReduce(nChangedTotal, sumOp<label>())
|
|
<< " coupled patch points during syncPatchDisplacement." << endl;
|
|
}
|
|
|
|
|
|
// Calculate displacement vector for all patch points. Uses pointNormal.
|
|
// Checks that displaced patch point would be visible from all centres
|
|
// of the faces using it.
|
|
// extrudeStatus is both input and output and gives the status of each
|
|
// patch point.
|
|
void Foam::autoHexMeshDriver::getPatchDisplacement
|
|
(
|
|
const motionSmoother& meshMover,
|
|
const scalarField& thickness,
|
|
const scalarField& minThickness,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
Info<< nl << "Determining displacement for added points"
|
|
<< " according to pointNormal ..." << endl;
|
|
|
|
const indirectPrimitivePatch& pp = meshMover.patch();
|
|
const vectorField& faceNormals = pp.faceNormals();
|
|
const labelListList& pointFaces = pp.pointFaces();
|
|
const pointField& localPoints = pp.localPoints();
|
|
const labelList& meshPoints = pp.meshPoints();
|
|
|
|
// Determine pointNormal
|
|
// ~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
pointField pointNormals(pp.nPoints(), vector::zero);
|
|
{
|
|
labelList nPointFaces(pp.nPoints(), 0);
|
|
|
|
forAll(faceNormals, faceI)
|
|
{
|
|
const face& f = pp.localFaces()[faceI];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
pointNormals[f[fp]] += faceNormals[faceI];
|
|
nPointFaces[f[fp]] ++;
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh_,
|
|
meshPoints,
|
|
pointNormals,
|
|
plusEqOp<vector>(),
|
|
vector::zero, // null value
|
|
false // no separation
|
|
);
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh_,
|
|
meshPoints,
|
|
nPointFaces,
|
|
plusEqOp<label>(),
|
|
0, // null value
|
|
false // no separation
|
|
);
|
|
|
|
forAll(pointNormals, i)
|
|
{
|
|
pointNormals[i] /= nPointFaces[i];
|
|
}
|
|
}
|
|
|
|
|
|
// Determine local length scale on patch
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Start off from same thickness everywhere (except where no extrusion)
|
|
patchDisp = thickness*pointNormals;
|
|
|
|
// Check if no extrude possible.
|
|
forAll(pointNormals, patchPointI)
|
|
{
|
|
label meshPointI = pp.meshPoints()[patchPointI];
|
|
|
|
if (extrudeStatus[patchPointI] == NOEXTRUDE)
|
|
{
|
|
// Do not use unmarkExtrusion; forcibly set to zero extrusion.
|
|
patchNLayers[patchPointI] = 0;
|
|
patchDisp[patchPointI] = vector::zero;
|
|
}
|
|
else
|
|
{
|
|
// Get normal
|
|
const vector& n = pointNormals[patchPointI];
|
|
|
|
if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointI]))
|
|
{
|
|
Pout<< "No valid normal for point " << meshPointI
|
|
<< ' ' << pp.points()[meshPointI]
|
|
<< "; setting displacement to " << patchDisp[patchPointI]
|
|
<< endl;
|
|
|
|
extrudeStatus[patchPointI] = EXTRUDEREMOVE;
|
|
}
|
|
}
|
|
}
|
|
|
|
// At illegal points make displacement average of new neighbour positions
|
|
forAll(extrudeStatus, patchPointI)
|
|
{
|
|
if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
|
|
{
|
|
point avg(vector::zero);
|
|
label nPoints = 0;
|
|
|
|
const labelList& pEdges = pp.pointEdges()[patchPointI];
|
|
|
|
forAll(pEdges, i)
|
|
{
|
|
label edgeI = pEdges[i];
|
|
|
|
label otherPointI = pp.edges()[edgeI].otherVertex(patchPointI);
|
|
|
|
if (extrudeStatus[otherPointI] != NOEXTRUDE)
|
|
{
|
|
avg += localPoints[otherPointI] + patchDisp[otherPointI];
|
|
nPoints++;
|
|
}
|
|
}
|
|
|
|
if (nPoints > 0)
|
|
{
|
|
Pout<< "Displacement at illegal point "
|
|
<< localPoints[patchPointI]
|
|
<< " set to " << (avg / nPoints - localPoints[patchPointI])
|
|
<< endl;
|
|
|
|
patchDisp[patchPointI] =
|
|
avg / nPoints
|
|
- localPoints[patchPointI];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Make sure displacement is equal on both sides of coupled patches.
|
|
syncPatchDisplacement
|
|
(
|
|
meshMover,
|
|
minThickness,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
Info<< endl;
|
|
}
|
|
|
|
|
|
// Truncates displacement
|
|
// - for all patchFaces in the faceset displacement gets set to zero
|
|
// - all displacement < minThickness gets set to zero
|
|
Foam::label Foam::autoHexMeshDriver::truncateDisplacement
|
|
(
|
|
const motionSmoother& meshMover,
|
|
const scalarField& minThickness,
|
|
const faceSet& illegalPatchFaces,
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
) const
|
|
{
|
|
const polyMesh& mesh = meshMover.mesh();
|
|
const indirectPrimitivePatch& pp = meshMover.patch();
|
|
|
|
label nChanged = 0;
|
|
|
|
const Map<label>& meshPointMap = pp.meshPointMap();
|
|
|
|
forAllConstIter(faceSet, illegalPatchFaces, iter)
|
|
{
|
|
label faceI = iter.key();
|
|
|
|
if (mesh.isInternalFace(faceI))
|
|
{
|
|
FatalErrorIn("truncateDisplacement(..)")
|
|
<< "Faceset " << illegalPatchFaces.name()
|
|
<< " contains internal face " << faceI << nl
|
|
<< "It should only contain patch faces" << abort(FatalError);
|
|
}
|
|
|
|
const face& f = mesh.faces()[faceI];
|
|
|
|
|
|
forAll(f, fp)
|
|
{
|
|
if (meshPointMap.found(f[fp]))
|
|
{
|
|
label patchPointI = meshPointMap[f[fp]];
|
|
|
|
if (extrudeStatus[patchPointI] != NOEXTRUDE)
|
|
{
|
|
unmarkExtrusion
|
|
(
|
|
patchPointI,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
forAll(patchDisp, patchPointI)
|
|
{
|
|
if (mag(patchDisp[patchPointI]) < minThickness[patchPointI])
|
|
{
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
patchPointI,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
else if (extrudeStatus[patchPointI] == NOEXTRUDE)
|
|
{
|
|
// Make sure displacement is 0. Should already be so but ...
|
|
patchDisp[patchPointI] = vector::zero;
|
|
patchNLayers[patchPointI] = 0;
|
|
}
|
|
}
|
|
|
|
|
|
const faceList& localFaces = pp.localFaces();
|
|
|
|
while (true)
|
|
{
|
|
syncPatchDisplacement
|
|
(
|
|
meshMover,
|
|
minThickness,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Make sure that a face doesn't have two non-consecutive areas
|
|
// not extruded (e.g. quad where vertex 0 and 2 are not extruded
|
|
// but 1 and 3 are) since this gives topological errors.
|
|
|
|
label nPinched = 0;
|
|
|
|
forAll(localFaces, i)
|
|
{
|
|
const face& localF = localFaces[i];
|
|
|
|
// Count number of transitions from unsnapped to snapped.
|
|
label nTrans = 0;
|
|
|
|
extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
|
|
|
|
forAll(localF, fp)
|
|
{
|
|
extrudeMode fpMode = extrudeStatus[localF[fp]];
|
|
|
|
if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
|
|
{
|
|
nTrans++;
|
|
}
|
|
prevMode = fpMode;
|
|
}
|
|
|
|
if (nTrans > 1)
|
|
{
|
|
// Multiple pinches. Reset whole face as unextruded.
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
localF,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nPinched++;
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
|
|
reduce(nPinched, sumOp<label>());
|
|
|
|
Info<< "truncateDisplacement : Unextruded " << nPinched
|
|
<< " faces due to non-consecutive vertices being extruded." << endl;
|
|
|
|
|
|
// Make sure that a face has consistent number of layers for all
|
|
// its vertices.
|
|
|
|
label nDiffering = 0;
|
|
|
|
//forAll(localFaces, i)
|
|
//{
|
|
// const face& localF = localFaces[i];
|
|
//
|
|
// label numLayers = -1;
|
|
//
|
|
// forAll(localF, fp)
|
|
// {
|
|
// if (patchNLayers[localF[fp]] > 0)
|
|
// {
|
|
// if (numLayers == -1)
|
|
// {
|
|
// numLayers = patchNLayers[localF[fp]];
|
|
// }
|
|
// else if (numLayers != patchNLayers[localF[fp]])
|
|
// {
|
|
// // Differing number of layers
|
|
// if
|
|
// (
|
|
// unmarkExtrusion
|
|
// (
|
|
// localF,
|
|
// patchDisp,
|
|
// patchNLayers,
|
|
// extrudeStatus
|
|
// )
|
|
// )
|
|
// {
|
|
// nDiffering++;
|
|
// nChanged++;
|
|
// }
|
|
// break;
|
|
// }
|
|
// }
|
|
// }
|
|
//}
|
|
//
|
|
//reduce(nDiffering, sumOp<label>());
|
|
//
|
|
//Info<< "truncateDisplacement : Unextruded " << nDiffering
|
|
// << " faces due to having differing number of layers." << endl;
|
|
|
|
if (nPinched+nDiffering == 0)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
return nChanged;
|
|
}
|
|
|
|
|
|
// Setup layer information (at points and faces) to modify mesh topology in
|
|
// regions where layer mesh terminates.
|
|
void Foam::autoHexMeshDriver::setupLayerInfoTruncation
|
|
(
|
|
const motionSmoother& meshMover,
|
|
const labelList& patchNLayers,
|
|
const List<extrudeMode>& extrudeStatus,
|
|
const label nBufferCellsNoExtrude,
|
|
labelList& nPatchPointLayers,
|
|
labelList& nPatchFaceLayers
|
|
) const
|
|
{
|
|
Info<< nl << "Setting up information for layer truncation ..." << endl;
|
|
|
|
const indirectPrimitivePatch& pp = meshMover.patch();
|
|
const polyMesh& mesh = meshMover.mesh();
|
|
|
|
if (nBufferCellsNoExtrude < 0)
|
|
{
|
|
Info<< nl << "Performing no layer truncation."
|
|
<< " nBufferCellsNoExtrude set to less than 0 ..." << endl;
|
|
|
|
// Face layers if any point get extruded
|
|
forAll(pp.localFaces(), patchFaceI)
|
|
{
|
|
const face& f = pp.localFaces()[patchFaceI];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
if (patchNLayers[f[fp]] > 0)
|
|
{
|
|
nPatchFaceLayers[patchFaceI] = patchNLayers[f[fp]];
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
nPatchPointLayers = patchNLayers;
|
|
}
|
|
else
|
|
{
|
|
// Determine max point layers per face.
|
|
labelList maxLevel(pp.size(), 0);
|
|
|
|
forAll(pp.localFaces(), patchFaceI)
|
|
{
|
|
const face& f = pp.localFaces()[patchFaceI];
|
|
|
|
// find patch faces where layer terminates (i.e contains extrude
|
|
// and noextrude points).
|
|
|
|
bool noExtrude = false;
|
|
label mLevel = 0;
|
|
|
|
forAll(f, fp)
|
|
{
|
|
if (extrudeStatus[f[fp]] == NOEXTRUDE)
|
|
{
|
|
noExtrude = true;
|
|
}
|
|
mLevel = max(mLevel, patchNLayers[f[fp]]);
|
|
}
|
|
|
|
if (mLevel > 0)
|
|
{
|
|
// So one of the points is extruded. Check if all are extruded
|
|
// or is a mix.
|
|
|
|
if (noExtrude)
|
|
{
|
|
nPatchFaceLayers[patchFaceI] = 1;
|
|
maxLevel[patchFaceI] = mLevel;
|
|
}
|
|
else
|
|
{
|
|
maxLevel[patchFaceI] = mLevel;
|
|
}
|
|
}
|
|
}
|
|
|
|
// We have the seed faces (faces with nPatchFaceLayers != maxLevel)
|
|
// Now do a meshwave across the patch where we pick up neighbours
|
|
// of seed faces.
|
|
// Note: quite inefficient. Could probably be coded better.
|
|
|
|
const labelListList& pointFaces = pp.pointFaces();
|
|
|
|
label nLevels = gMax(patchNLayers);
|
|
|
|
// flag neighbouring patch faces with number of layers to grow
|
|
for (label ilevel = 1; ilevel < nLevels; ilevel++)
|
|
{
|
|
label nBuffer;
|
|
|
|
if (ilevel == 1)
|
|
{
|
|
nBuffer = nBufferCellsNoExtrude - 1;
|
|
}
|
|
else
|
|
{
|
|
nBuffer = nBufferCellsNoExtrude;
|
|
}
|
|
|
|
for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
|
|
{
|
|
labelList tempCounter(nPatchFaceLayers);
|
|
|
|
boolList foundNeighbour(pp.nPoints(), false);
|
|
|
|
forAll(pp.meshPoints(), patchPointI)
|
|
{
|
|
forAll(pointFaces[patchPointI], pointFaceI)
|
|
{
|
|
label faceI = pointFaces[patchPointI][pointFaceI];
|
|
|
|
if
|
|
(
|
|
nPatchFaceLayers[faceI] != -1
|
|
&& maxLevel[faceI] > 0
|
|
)
|
|
{
|
|
foundNeighbour[patchPointI] = true;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
foundNeighbour,
|
|
orEqOp<bool>(),
|
|
false, // null value
|
|
false // no separation
|
|
);
|
|
|
|
forAll(pp.meshPoints(), patchPointI)
|
|
{
|
|
if (foundNeighbour[patchPointI])
|
|
{
|
|
forAll(pointFaces[patchPointI], pointFaceI)
|
|
{
|
|
label faceI = pointFaces[patchPointI][pointFaceI];
|
|
if
|
|
(
|
|
nPatchFaceLayers[faceI] == -1
|
|
&& maxLevel[faceI] > 0
|
|
&& ilevel < maxLevel[faceI]
|
|
)
|
|
{
|
|
tempCounter[faceI] = ilevel;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
nPatchFaceLayers = tempCounter;
|
|
}
|
|
}
|
|
|
|
forAll(pp.localFaces(), patchFaceI)
|
|
{
|
|
if (nPatchFaceLayers[patchFaceI] == -1)
|
|
{
|
|
nPatchFaceLayers[patchFaceI] = maxLevel[patchFaceI];
|
|
}
|
|
}
|
|
|
|
forAll(pp.meshPoints(), patchPointI)
|
|
{
|
|
if (extrudeStatus[patchPointI] != NOEXTRUDE)
|
|
{
|
|
forAll(pointFaces[patchPointI], pointFaceI)
|
|
{
|
|
label face = pointFaces[patchPointI][pointFaceI];
|
|
nPatchPointLayers[patchPointI] = max
|
|
(
|
|
nPatchPointLayers[patchPointI],
|
|
nPatchFaceLayers[face]
|
|
);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
nPatchPointLayers[patchPointI] = 0;
|
|
}
|
|
}
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
nPatchPointLayers,
|
|
maxEqOp<label>(),
|
|
0, // null value
|
|
false // no separation
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// Does any of the cells use a face from faces?
|
|
bool Foam::autoHexMeshDriver::cellsUseFace
|
|
(
|
|
const polyMesh& mesh,
|
|
const labelList& cellLabels,
|
|
const labelHashSet& faces
|
|
)
|
|
{
|
|
forAll(cellLabels, i)
|
|
{
|
|
const cell& cFaces = mesh.cells()[cellLabels[i]];
|
|
|
|
forAll(cFaces, cFaceI)
|
|
{
|
|
if (faces.found(cFaces[cFaceI]))
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
|
|
// Checks the newly added cells and locally unmarks points so they
|
|
// will not get extruded next time round. Returns global number of unmarked
|
|
// points (0 if all was fine)
|
|
Foam::label Foam::autoHexMeshDriver::checkAndUnmark
|
|
(
|
|
const addPatchCellLayer& addLayer,
|
|
const dictionary& motionDict,
|
|
const indirectPrimitivePatch& pp,
|
|
const fvMesh& newMesh,
|
|
|
|
pointField& patchDisp,
|
|
labelList& patchNLayers,
|
|
List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
// Check the resulting mesh for errors
|
|
Info<< nl << "Checking mesh with layer ..." << endl;
|
|
faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
|
|
motionSmoother::checkMesh(false, newMesh, motionDict, wrongFaces);
|
|
Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
|
|
<< " illegal faces"
|
|
<< " (concave, zero area or negative cell pyramid volume)"
|
|
<< endl;
|
|
|
|
// Undo local extrusion if
|
|
// - any of the added cells in error
|
|
|
|
label nChanged = 0;
|
|
|
|
// Get all cells in the layer.
|
|
labelListList addedCells
|
|
(
|
|
addPatchCellLayer::addedCells
|
|
(
|
|
newMesh,
|
|
addLayer.layerFaces()
|
|
)
|
|
);
|
|
|
|
// Check if any of the faces in error uses any face of an added cell
|
|
forAll(addedCells, oldPatchFaceI)
|
|
{
|
|
// Get the cells (in newMesh labels) per old patch face (in mesh
|
|
// labels)
|
|
const labelList& fCells = addedCells[oldPatchFaceI];
|
|
|
|
if (cellsUseFace(newMesh, fCells, wrongFaces))
|
|
{
|
|
// Unmark points on old mesh
|
|
if
|
|
(
|
|
unmarkExtrusion
|
|
(
|
|
pp.localFaces()[oldPatchFaceI],
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
)
|
|
)
|
|
{
|
|
nChanged++;
|
|
}
|
|
}
|
|
}
|
|
|
|
return returnReduce(nChanged, sumOp<label>());
|
|
}
|
|
|
|
|
|
//- Count global number of extruded faces
|
|
Foam::label Foam::autoHexMeshDriver::countExtrusion
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const List<extrudeMode>& extrudeStatus
|
|
)
|
|
{
|
|
// Count number of extruded patch faces
|
|
label nExtruded = 0;
|
|
{
|
|
const faceList& localFaces = pp.localFaces();
|
|
|
|
forAll(localFaces, i)
|
|
{
|
|
const face& localFace = localFaces[i];
|
|
|
|
forAll(localFace, fp)
|
|
{
|
|
if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
|
|
{
|
|
nExtruded++;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return returnReduce(nExtruded, sumOp<label>());
|
|
}
|
|
|
|
|
|
// Collect layer faces and layer cells into bools for ease of handling
|
|
void Foam::autoHexMeshDriver::getLayerCellsFaces
|
|
(
|
|
const polyMesh& mesh,
|
|
const addPatchCellLayer& addLayer,
|
|
boolList& flaggedCells,
|
|
boolList& flaggedFaces
|
|
)
|
|
{
|
|
flaggedCells.setSize(mesh.nCells());
|
|
flaggedCells = false;
|
|
flaggedFaces.setSize(mesh.nFaces());
|
|
flaggedFaces = false;
|
|
|
|
// Mark all faces in the layer
|
|
const labelListList& layerFaces = addLayer.layerFaces();
|
|
|
|
// Mark all cells in the layer.
|
|
labelListList addedCells(addPatchCellLayer::addedCells(mesh, layerFaces));
|
|
|
|
forAll(addedCells, oldPatchFaceI)
|
|
{
|
|
const labelList& added = addedCells[oldPatchFaceI];
|
|
|
|
forAll(added, i)
|
|
{
|
|
flaggedCells[added[i]] = true;
|
|
}
|
|
}
|
|
|
|
forAll(layerFaces, oldPatchFaceI)
|
|
{
|
|
const labelList& layer = layerFaces[oldPatchFaceI];
|
|
|
|
if (layer.size() > 0)
|
|
{
|
|
for (label i = 1; i < layer.size()-1; i++)
|
|
{
|
|
flaggedFaces[layer[i]] = true;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
void Foam::autoHexMeshDriver::mergePatchFacesUndo
|
|
(
|
|
const dictionary& shrinkDict,
|
|
const dictionary& motionDict
|
|
)
|
|
{
|
|
const scalar featureAngle(readScalar(shrinkDict.lookup("featureAngle")));
|
|
scalar minCos = Foam::cos(featureAngle*mathematicalConstant::pi/180.0);
|
|
|
|
scalar concaveAngle = defaultConcaveAngle;
|
|
|
|
if (shrinkDict.found("concaveAngle"))
|
|
{
|
|
concaveAngle = readScalar(shrinkDict.lookup("concaveAngle"));
|
|
}
|
|
scalar concaveCos = Foam::cos(concaveAngle*mathematicalConstant::pi/180.0);
|
|
|
|
Info<< nl
|
|
<< "Merging all faces of a cell" << nl
|
|
<< "---------------------------" << nl
|
|
<< " - which are on the same patch" << nl
|
|
<< " - which make an angle < " << featureAngle << " degrees"
|
|
<< nl
|
|
<< " (cos:" << minCos << ')' << nl
|
|
<< " - as long as the resulting face doesn't become concave"
|
|
<< " by more than "
|
|
<< concaveAngle << " degrees (0=straight, 180=fully concave)" << nl
|
|
<< endl;
|
|
|
|
label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
|
|
|
|
nChanged += mergeEdgesUndo(minCos, motionDict);
|
|
}
|
|
|
|
|
|
void Foam::autoHexMeshDriver::addLayers
|
|
(
|
|
const dictionary& shrinkDict,
|
|
const dictionary& motionDict,
|
|
const label nAllowableErrors,
|
|
motionSmoother& meshMover
|
|
)
|
|
{
|
|
// Read some more dictionary settings
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Min thickness per cell
|
|
const scalar relMinThickness(readScalar(shrinkDict.lookup("minThickness")));
|
|
|
|
// Warped faces recoginition
|
|
const scalar maxFaceThicknessRatio
|
|
(
|
|
readScalar(shrinkDict.lookup("maxFaceThicknessRatio"))
|
|
);
|
|
const scalar featureAngle(readScalar(shrinkDict.lookup("featureAngle")));
|
|
|
|
// Feature angle
|
|
const scalar minCos = Foam::cos(featureAngle*mathematicalConstant::pi/180.);
|
|
|
|
// Number of layers for to grow non-extrusion region by
|
|
const label nGrow(readLabel(shrinkDict.lookup("nGrow")));
|
|
|
|
//const label nSmoothDisp(readLabel(shrinkDict.lookup("nSmoothDispl")));
|
|
// Snapping iterations
|
|
const label nSnap(readLabel(shrinkDict.lookup("nSnap")));
|
|
|
|
// Mesh termination and medial axis smoothing quantities
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
const scalar minCosLayerTermination =
|
|
Foam::cos(0.5*featureAngle*mathematicalConstant::pi/180.);
|
|
const scalar expansionRatio
|
|
(
|
|
readScalar(shrinkDict.lookup("expansionRatio"))
|
|
);
|
|
const scalar finalLayerRatio
|
|
(
|
|
readScalar(shrinkDict.lookup("finalLayerRatio"))
|
|
);
|
|
const label nBufferCellsNoExtrude
|
|
(
|
|
readLabel(shrinkDict.lookup("nBufferCellsNoExtrude"))
|
|
);
|
|
const label nSmoothSurfaceNormals
|
|
(
|
|
readLabel(shrinkDict.lookup("nSmoothSurfaceNormals"))
|
|
);
|
|
const label nSmoothNormals(readLabel(shrinkDict.lookup("nSmoothNormals")));
|
|
const label nSmoothThickness
|
|
(
|
|
readLabel(shrinkDict.lookup("nSmoothThickness"))
|
|
);
|
|
const scalar maxThicknessToMedialRatio
|
|
(
|
|
readScalar(shrinkDict.lookup("maxThicknessToMedialRatio"))
|
|
);
|
|
|
|
// Medial axis setup
|
|
|
|
const scalar minMedianAxisAngle
|
|
(
|
|
readScalar(shrinkDict.lookup("minMedianAxisAngle"))
|
|
);
|
|
const scalar minMedianAxisAngleCos =
|
|
Foam::cos(minMedianAxisAngle*mathematicalConstant::pi/180.);
|
|
|
|
|
|
// Precalculate mesh edge labels for patch edges
|
|
const indirectPrimitivePatch& pp = meshMover.patch();
|
|
const labelList& meshPoints = pp.meshPoints();
|
|
|
|
labelList meshEdges(pp.nEdges());
|
|
forAll(pp.edges(), edgeI)
|
|
{
|
|
const edge& ppEdge = pp.edges()[edgeI];
|
|
label v0 = meshPoints[ppEdge[0]];
|
|
label v1 = meshPoints[ppEdge[1]];
|
|
meshEdges[edgeI] = meshTools::findEdge
|
|
(
|
|
mesh_.edges(),
|
|
mesh_.pointEdges()[v0],
|
|
v0,
|
|
v1
|
|
);
|
|
}
|
|
|
|
// Displacement for all pp.localPoints.
|
|
vectorField patchDisp(pp.nPoints(), vector::one);
|
|
|
|
// Number of layers for all pp.localPoints. Note: only valid if
|
|
// extrudeStatus = EXTRUDE.
|
|
labelList patchNLayers(pp.nPoints(), 0);
|
|
|
|
// Whether to add edge for all pp.localPoints.
|
|
List<extrudeMode> extrudeStatus(pp.nPoints(), EXTRUDE);
|
|
|
|
|
|
// Get number of layer per point from number of layers per patch
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
setNumLayers
|
|
(
|
|
meshMover.adaptPatchIDs(),
|
|
pp,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Disable extrusion on non-manifold points
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
handleNonManifolds
|
|
(
|
|
pp,
|
|
meshEdges,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Disable extrusion on feature angles
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
handleFeatureAngle
|
|
(
|
|
pp,
|
|
meshEdges,
|
|
minCos,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Disable extrusion on warped faces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Undistorted edge length
|
|
const scalar edge0Len = meshRefinerPtr_().meshCutter().level0EdgeLength();
|
|
const labelList& cellLevel = meshRefinerPtr_().meshCutter().cellLevel();
|
|
|
|
handleWarpedFaces
|
|
(
|
|
pp,
|
|
maxFaceThicknessRatio,
|
|
edge0Len,
|
|
cellLevel,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
//// Disable extrusion on cells with multiple patch faces
|
|
//// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
//
|
|
//handleMultiplePatchFaces
|
|
//(
|
|
// pp,
|
|
//
|
|
// patchDisp,
|
|
// patchNLayers,
|
|
// extrudeStatus
|
|
//);
|
|
|
|
|
|
// Grow out region of non-extrusion
|
|
for (label i = 0; i < nGrow; i++)
|
|
{
|
|
growNoExtrusion
|
|
(
|
|
pp,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
}
|
|
|
|
// Determine point-wise layer thickness
|
|
scalarField thickness(pp.nPoints());
|
|
scalarField minThickness(pp.nPoints());
|
|
calculateLayerThickness
|
|
(
|
|
expansionRatio,
|
|
finalLayerRatio,
|
|
relMinThickness,
|
|
pp,
|
|
cellLevel,
|
|
patchNLayers,
|
|
edge0Len,
|
|
|
|
thickness,
|
|
minThickness
|
|
);
|
|
|
|
// Calculate wall to medial axis distance for smoothing displacement
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
pointScalarField pointMedialDist
|
|
(
|
|
IOobject
|
|
(
|
|
"pointMedialDist",
|
|
mesh_.time().timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
meshMover.pMesh(),
|
|
dimensionedScalar("pointMedialDist", dimless, 0.0)
|
|
);
|
|
|
|
pointVectorField dispVec
|
|
(
|
|
IOobject
|
|
(
|
|
"dispVec",
|
|
mesh_.time().timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
meshMover.pMesh(),
|
|
dimensionedVector("dispVec", dimless, vector::zero)
|
|
);
|
|
|
|
pointScalarField medialRatio
|
|
(
|
|
IOobject
|
|
(
|
|
"medialRatio",
|
|
mesh_.time().timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
meshMover.pMesh(),
|
|
dimensionedScalar("medialRatio", dimless, 0.0)
|
|
);
|
|
|
|
// Setup information for medial axis smoothing. Calculates medial axis
|
|
// and a smoothed displacement direction.
|
|
// - pointMedialDist : distance to medial axis
|
|
// - dispVec : normalised direction of nearest displacement
|
|
// - medialRatio : ratio of medial distance to wall distance.
|
|
// (1 at wall, 0 at medial axis)
|
|
medialAxisSmoothingInfo
|
|
(
|
|
meshMover,
|
|
nSmoothNormals,
|
|
nSmoothSurfaceNormals,
|
|
minMedianAxisAngleCos,
|
|
|
|
dispVec,
|
|
medialRatio,
|
|
pointMedialDist
|
|
);
|
|
|
|
|
|
|
|
// Saved old points
|
|
pointField oldPoints(mesh_.points());
|
|
|
|
// Last set of topology changes. (changing mesh clears out polyTopoChange)
|
|
polyTopoChange savedMeshMod(mesh_.boundaryMesh().size());
|
|
|
|
boolList flaggedCells;
|
|
boolList flaggedFaces;
|
|
|
|
while (true)
|
|
{
|
|
// Make sure displacement is equal on both sides of coupled patches.
|
|
syncPatchDisplacement
|
|
(
|
|
meshMover,
|
|
minThickness,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Displacement acc. to pointnormals
|
|
getPatchDisplacement
|
|
(
|
|
meshMover,
|
|
thickness,
|
|
minThickness,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
// Shrink mesh by displacement value first.
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
{
|
|
pointField oldPatchPos(pp.localPoints());
|
|
|
|
//// Laplacian displacement shrinking.
|
|
//shrinkMeshDistance
|
|
//(
|
|
// debug,
|
|
// meshMover,
|
|
// -patchDisp, // Shrink in opposite direction of addedPoints
|
|
// nSmoothDisp,
|
|
// nSnap
|
|
//);
|
|
|
|
// Medial axis based shrinking
|
|
shrinkMeshMedialDistance
|
|
(
|
|
meshMover,
|
|
|
|
nSmoothThickness,
|
|
maxThicknessToMedialRatio,
|
|
nAllowableErrors,
|
|
nSnap,
|
|
minCosLayerTermination,
|
|
|
|
thickness,
|
|
minThickness,
|
|
|
|
dispVec,
|
|
medialRatio,
|
|
pointMedialDist,
|
|
|
|
extrudeStatus,
|
|
patchDisp,
|
|
patchNLayers
|
|
);
|
|
|
|
// Update patchDisp (since not all might have been honoured)
|
|
patchDisp = oldPatchPos - pp.localPoints();
|
|
}
|
|
|
|
// Truncate displacements that are too small (this will do internal
|
|
// ones, coupled ones have already been truncated by syncPatch)
|
|
faceSet dummySet(mesh_, "wrongPatchFaces", 0);
|
|
truncateDisplacement
|
|
(
|
|
meshMover,
|
|
minThickness,
|
|
dummySet,
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
|
|
// Dump to .obj file for debugging.
|
|
if (debug_)
|
|
{
|
|
dumpDisplacement
|
|
(
|
|
mesh_.time().path()/"layer",
|
|
pp,
|
|
patchDisp,
|
|
extrudeStatus
|
|
);
|
|
|
|
const_cast<Time&>(mesh_.time())++;
|
|
Info<< "Writing shrunk mesh to " << mesh_.time().timeName() << endl;
|
|
mesh_.write();
|
|
}
|
|
|
|
|
|
// Mesh topo change engine
|
|
polyTopoChange meshMod(mesh_);
|
|
|
|
// Grow layer of cells on to patch. Handles zero sized displacement.
|
|
addPatchCellLayer addLayer(mesh_);
|
|
|
|
// Determine per point/per face number of layers to extrude. Also
|
|
// handles the slow termination of layers when going switching layers
|
|
|
|
labelList nPatchPointLayers(pp.nPoints(),-1);
|
|
labelList nPatchFaceLayers(pp.localFaces().size(),-1);
|
|
setupLayerInfoTruncation
|
|
(
|
|
meshMover,
|
|
patchNLayers,
|
|
extrudeStatus,
|
|
nBufferCellsNoExtrude,
|
|
nPatchPointLayers,
|
|
nPatchFaceLayers
|
|
);
|
|
|
|
// Calculate displacement for first layer for addPatchLayer
|
|
vectorField firstDisp(patchNLayers.size(), vector::zero);
|
|
|
|
forAll(patchNLayers, i)
|
|
{
|
|
if (patchNLayers[i] > 0)
|
|
{
|
|
if (expansionRatio == 1.0)
|
|
{
|
|
firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
|
|
}
|
|
else
|
|
{
|
|
label nLay = nPatchPointLayers[i];
|
|
scalar h =
|
|
pow(expansionRatio,nLay - 1)
|
|
* (mag(patchDisp[i])*(1.0 - expansionRatio))
|
|
/ (1.0 - pow(expansionRatio, nLay));
|
|
firstDisp[i] = h/mag(patchDisp[i])*patchDisp[i];
|
|
}
|
|
}
|
|
}
|
|
|
|
scalar invExpansionRatio = 1.0 / expansionRatio;
|
|
|
|
// Add topo regardless of whether extrudeStatus is extruderemove.
|
|
// Not add layer if patchDisp is zero.
|
|
addLayer.setRefinement
|
|
(
|
|
invExpansionRatio,
|
|
pp,
|
|
nPatchFaceLayers, // layers per face
|
|
nPatchPointLayers, // layers per point
|
|
firstDisp, // thickness of first layer
|
|
meshMod
|
|
);
|
|
|
|
if (debug_)
|
|
{
|
|
const_cast<Time&>(mesh_.time())++;
|
|
}
|
|
|
|
// Store mesh changes for if mesh is correct.
|
|
savedMeshMod = meshMod;
|
|
|
|
|
|
// With the stored topo changes we create a new mesh so we can
|
|
// undo if neccesary.
|
|
|
|
autoPtr<fvMesh> newMeshPtr;
|
|
autoPtr<mapPolyMesh> map = meshMod.makeMesh
|
|
(
|
|
newMeshPtr,
|
|
IOobject
|
|
(
|
|
//mesh.name()+"_layer",
|
|
mesh_.name(),
|
|
static_cast<polyMesh&>(mesh_).instance(),
|
|
mesh_.db(),
|
|
static_cast<polyMesh&>(mesh_).readOpt(),
|
|
static_cast<polyMesh&>(mesh_).writeOpt()
|
|
), // io params from original mesh but new name
|
|
mesh_, // original mesh
|
|
true // parallel sync
|
|
);
|
|
fvMesh& newMesh = newMeshPtr();
|
|
|
|
//?neccesary? Update fields
|
|
newMesh.updateMesh(map);
|
|
|
|
// Update numbering on addLayer:
|
|
// - cell/point labels to be newMesh.
|
|
// - patchFaces to remain in oldMesh order.
|
|
addLayer.updateMesh
|
|
(
|
|
map,
|
|
identity(pp.size()),
|
|
identity(pp.nPoints())
|
|
);
|
|
|
|
// Collect layer faces and cells for outside loop.
|
|
getLayerCellsFaces
|
|
(
|
|
newMesh,
|
|
addLayer,
|
|
flaggedCells,
|
|
flaggedFaces
|
|
);
|
|
|
|
|
|
if (debug_)
|
|
{
|
|
Info<< "Writing layer mesh to " << mesh_.time().timeName() << endl;
|
|
newMesh.write();
|
|
cellSet addedCellSet
|
|
(
|
|
newMesh,
|
|
"addedCells",
|
|
findIndices(flaggedCells, true)
|
|
);
|
|
Info<< "Writing "
|
|
<< returnReduce(addedCellSet.size(), sumOp<label>())
|
|
<< " added cells to cellSet "
|
|
<< addedCellSet.objectPath()
|
|
<< endl;
|
|
addedCellSet.write();
|
|
|
|
faceSet layerFacesSet
|
|
(
|
|
newMesh,
|
|
"layerFaces",
|
|
findIndices(flaggedCells, true)
|
|
);
|
|
Info<< "Writing "
|
|
<< returnReduce(layerFacesSet.size(), sumOp<label>())
|
|
<< " faces inside added layer to faceSet "
|
|
<< layerFacesSet.objectPath()
|
|
<< endl;
|
|
layerFacesSet.write();
|
|
}
|
|
|
|
// Unset the extrusion at the pp.
|
|
label nTotChanged = checkAndUnmark
|
|
(
|
|
addLayer,
|
|
motionDict,
|
|
pp,
|
|
newMesh,
|
|
|
|
patchDisp,
|
|
patchNLayers,
|
|
extrudeStatus
|
|
);
|
|
|
|
Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
|
|
<< " out of " << returnReduce(pp.size(), sumOp<label>())
|
|
<< " faces. Removed extrusion at " << nTotChanged << " faces."
|
|
<< endl;
|
|
|
|
if (nTotChanged == 0)
|
|
{
|
|
break;
|
|
}
|
|
|
|
// Reset mesh points and start again
|
|
meshMover.movePoints(oldPoints);
|
|
meshMover.correct();
|
|
}
|
|
|
|
|
|
// At this point we have a (shrunk) mesh and a set of topology changes
|
|
// which will make a valid mesh with layer. Apply these changes to the
|
|
// current mesh.
|
|
|
|
// Apply the stored topo changes to the current mesh.
|
|
autoPtr<mapPolyMesh> map = savedMeshMod.changeMesh(mesh_, false);
|
|
|
|
// Update fields
|
|
mesh_.updateMesh(map);
|
|
|
|
// Move mesh (since morphing does not do this)
|
|
if (map().hasMotionPoints())
|
|
{
|
|
mesh_.movePoints(map().preMotionPoints());
|
|
}
|
|
|
|
meshRefinerPtr_().updateMesh(map, labelList(0));
|
|
|
|
|
|
// Do final balancing
|
|
// ~~~~~~~~~~~~~~~~~~
|
|
|
|
if (Pstream::parRun())
|
|
{
|
|
// Balance. No restriction on face zones and baffles.
|
|
autoPtr<mapDistributePolyMesh> map = balance(false, false);
|
|
|
|
// Re-distribute flag of layer faces and cells
|
|
map().distributeCellData(flaggedCells);
|
|
map().distributeFaceData(flaggedFaces);
|
|
}
|
|
|
|
|
|
// Write mesh
|
|
// ~~~~~~~~~~
|
|
|
|
//writeMesh("Layer mesh");
|
|
cellSet addedCellSet(mesh_, "addedCells", findIndices(flaggedCells, true));
|
|
Info<< "Writing "
|
|
<< returnReduce(addedCellSet.size(), sumOp<label>())
|
|
<< " added cells to cellSet "
|
|
<< addedCellSet.objectPath()
|
|
<< endl;
|
|
addedCellSet.write();
|
|
|
|
faceSet layerFacesSet(mesh_, "layerFaces", findIndices(flaggedFaces, true));
|
|
Info<< "Writing "
|
|
<< returnReduce(layerFacesSet.size(), sumOp<label>())
|
|
<< " faces inside added layer to faceSet "
|
|
<< layerFacesSet.objectPath()
|
|
<< endl;
|
|
layerFacesSet.write();
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|