zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: snappyHexMesh: initial feature-line support

Browse Source
This commit is contained in:
mattijs
2011-04-05 11:56:57 +01:00
parent 22e752d348
commit 3e71574d8a
20 changed files with 3835 additions and 1054 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -37,6 +37,7 @@ Description
#include "autoLayerDriver.H" #include "autoLayerDriver.H"
#include "searchableSurfaces.H" #include "searchableSurfaces.H"
#include "refinementSurfaces.H" #include "refinementSurfaces.H"
#include "refinementFeatures.H"
#include "shellSurfaces.H" #include "shellSurfaces.H"
#include "decompositionMethod.H" #include "decompositionMethod.H"
#include "fvMeshDistribute.H" #include "fvMeshDistribute.H"
@ -250,6 +251,20 @@ int main(int argc, char *argv[])
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl; << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
// Read feature meshes
// ~~~~~~~~~~~~~~~~~~~
Info<< "Reading features." << endl;
refinementFeatures features
(
mesh,
refineDict.lookup("features")
);
Info<< "Read features in = "
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
// Refinement engine // Refinement engine
// ~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~
@ -265,6 +280,7 @@ int main(int argc, char *argv[])
mergeDist, // tolerance used in sorting coordinates mergeDist, // tolerance used in sorting coordinates
overwrite, // overwrite mesh files? overwrite, // overwrite mesh files?
surfaces, // for surface intersection refinement surfaces, // for surface intersection refinement
features, // for feature edges/point based refinement
shells // for volume (inside/outside) refinement shells // for volume (inside/outside) refinement
); );
Info<< "Calculated surface intersections in = " Info<< "Calculated surface intersections in = "
@ -429,13 +445,19 @@ int main(int argc, char *argv[])
// Snap parameters // Snap parameters
snapParameters snapParams(snapDict); snapParameters snapParams(snapDict);
// Temporary hack to get access to resolveFeatureAngle
scalar curvature;
{
refinementParameters refineParams(refineDict);
curvature = refineParams.curvature();
}
if (!overwrite) if (!overwrite)
{ {
const_cast<Time&>(mesh.time())++; const_cast<Time&>(mesh.time())++;
} }
snapDriver.doSnap(snapDict, motionDict, snapParams); snapDriver.doSnap(snapDict, motionDict, curvature, snapParams);
writeMesh writeMesh
( (

8
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
View File

@ -217,10 +217,10 @@ snapControls
// to surface // to surface
nSmoothPatch 3; nSmoothPatch 3;
//- Relative distance for points to be attracted by surface feature point //- Maximum relative distance for points to be attracted by surface.
// or edge. True distance is this factor times local // True distance is this factor times local maximum edge length.
// maximum edge length. // Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
tolerance 4.0; tolerance 2.0;
//- Number of mesh displacement relaxation iterations. //- Number of mesh displacement relaxation iterations.
nSolveIter 30; nSolveIter 30;

2
src/mesh/autoMesh/Make/files
View File

@ -4,6 +4,7 @@ autoHexMeshDriver = $(autoHexMesh)/autoHexMeshDriver
$(autoHexMeshDriver)/autoLayerDriver.C $(autoHexMeshDriver)/autoLayerDriver.C
$(autoHexMeshDriver)/autoLayerDriverShrink.C $(autoHexMeshDriver)/autoLayerDriverShrink.C
$(autoHexMeshDriver)/autoSnapDriver.C $(autoHexMeshDriver)/autoSnapDriver.C
$(autoHexMeshDriver)/autoSnapDriverFeature.C
$(autoHexMeshDriver)/autoRefineDriver.C $(autoHexMeshDriver)/autoRefineDriver.C
$(autoHexMeshDriver)/layerParameters/layerParameters.C $(autoHexMeshDriver)/layerParameters/layerParameters.C
@ -16,6 +17,7 @@ $(autoHexMesh)/meshRefinement/meshRefinement.C
$(autoHexMesh)/meshRefinement/meshRefinementMerge.C $(autoHexMesh)/meshRefinement/meshRefinementMerge.C
$(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C $(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
$(autoHexMesh)/meshRefinement/meshRefinementRefine.C $(autoHexMesh)/meshRefinement/meshRefinementRefine.C
$(autoHexMesh)/refinementFeatures/refinementFeatures.C
$(autoHexMesh)/refinementSurfaces/refinementSurfaces.C $(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
$(autoHexMesh)/shellSurfaces/shellSurfaces.C $(autoHexMesh)/shellSurfaces/shellSurfaces.C
$(autoHexMesh)/trackedParticle/trackedParticle.C $(autoHexMesh)/trackedParticle/trackedParticle.C

638
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View File

@ -59,622 +59,6 @@ defineTypeNameAndDebug(autoLayerDriver, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
(
const scalar minCos,
const scalar concaveCos,
const dictionary& motionDict
)
{
fvMesh& mesh = meshRefiner_.mesh();
// Patch face merging engine
combineFaces faceCombiner(mesh, true);
// Pick up all candidate cells on boundary
labelHashSet boundaryCells(mesh.nFaces()-mesh.nInternalFaces());
{
labelList patchIDs(meshRefiner_.meshedPatches());
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.instance() = mesh.time().timeName();
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
meshRefiner_.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());
}
else
{
// Delete mesh volumes.
mesh.clearOut();
}
faceCombiner.updateMesh(map);
meshRefiner_.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)
{
errorFaces.instance() = mesh.time().timeName();
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);
restoreSet.instance() = mesh.time().timeName();
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());
}
else
{
// Delete mesh volumes.
mesh.clearOut();
}
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
meshRefiner_.updateMesh
(
map,
labelList(0), // changedFaces
restoredPoints,
restoredFaces,
restoredCells
);
Info<< endl;
}
if (debug)
{
Pout<< "Writing merged-faces mesh to time "
<< meshRefiner_.timeName() << nl << endl;
mesh.write();
}
}
else
{
Info<< "No faces merged ..." << endl;
}
return nFaceSets;
}
// Remove points. pointCanBeDeleted is parallel synchronised.
Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRemovePoints
(
removePoints& pointRemover,
const boolList& pointCanBeDeleted
)
{
fvMesh& mesh = meshRefiner_.mesh();
// 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());
}
else
{
// Delete mesh volumes.
mesh.clearOut();
}
pointRemover.updateMesh(map);
meshRefiner_.updateMesh(map, labelList(0));
return map;
}
// Restore faces (which contain removed points)
Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRestorePoints
(
removePoints& pointRemover,
const labelList& facesToRestore
)
{
fvMesh& mesh = meshRefiner_.mesh();
// 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());
}
else
{
// Delete mesh volumes.
mesh.clearOut();
}
pointRemover.updateMesh(map);
meshRefiner_.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::autoLayerDriver::collectFaces
(
const labelList& candidateFaces,
const labelHashSet& set
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
// 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
);
labelList selectedFaces(findIndices(selected, true));
return selectedFaces;
}
// Pick up faces of cells of faces in set.
Foam::labelList Foam::autoLayerDriver::growFaceCellFace
(
const labelHashSet& set
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
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
);
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::autoLayerDriver::mergeEdgesUndo
(
const scalar minCos,
const dictionary& motionDict
)
{
fvMesh& mesh = meshRefiner_.mesh();
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)
{
errorFaces.instance() = mesh.time().timeName();
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 "
<< meshRefiner_.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 // For debugging: Dump displacement to .obj files
void Foam::autoLayerDriver::dumpDisplacement void Foam::autoLayerDriver::dumpDisplacement
( (
@ -793,7 +177,7 @@ void Foam::autoLayerDriver::checkMeshManifold() const
<< " points where this happens to pointSet " << " points where this happens to pointSet "
<< nonManifoldPoints.name() << endl; << nonManifoldPoints.name() << endl;
nonManifoldPoints.instance() = mesh.time().timeName(); nonManifoldPoints.instance() = meshRefiner_.timeName();
nonManifoldPoints.write(); nonManifoldPoints.write();
} }
Info<< endl; Info<< endl;
@ -1080,7 +464,6 @@ void Foam::autoLayerDriver::handleWarpedFaces
} }
//// No extrusion on cells with multiple patch faces. There ususally is a reason //// No extrusion on cells with multiple patch faces. There ususally is a reason
//// why combinePatchFaces hasn't succeeded. //// why combinePatchFaces hasn't succeeded.
//void Foam::autoLayerDriver::handleMultiplePatchFaces //void Foam::autoLayerDriver::handleMultiplePatchFaces
@ -1133,7 +516,7 @@ void Foam::autoLayerDriver::handleWarpedFaces
// //
// if (nMultiPatchCells > 0) // if (nMultiPatchCells > 0)
// { // {
// multiPatchCells.instance() = mesh.time().timeName(); // multiPatchCells.instance() = meshRefiner_.timeName();
// Info<< "Writing " << nMultiPatchCells // Info<< "Writing " << nMultiPatchCells
// << " cells with multiple (connected) patch faces to cellSet " // << " cells with multiple (connected) patch faces to cellSet "
// << multiPatchCells.objectPath() << endl; // << multiPatchCells.objectPath() << endl;
@ -2547,9 +1930,15 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
<< " (0=straight, 180=fully concave)" << nl << " (0=straight, 180=fully concave)" << nl
<< endl; << endl;
label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict); label nChanged = meshRefiner_.mergePatchFacesUndo
(
minCos,
concaveCos,
meshRefiner_.meshedPatches(),
motionDict
);
nChanged += mergeEdgesUndo(minCos, motionDict); nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
} }
@ -3191,6 +2580,7 @@ void Foam::autoLayerDriver::addLayers
"addedCells", "addedCells",
findIndices(flaggedCells, true) findIndices(flaggedCells, true)
); );
addedCellSet.instance() = meshRefiner_.timeName();
Info<< "Writing " Info<< "Writing "
<< returnReduce(addedCellSet.size(), sumOp<label>()) << returnReduce(addedCellSet.size(), sumOp<label>())
<< " added cells to cellSet " << " added cells to cellSet "
@ -3203,6 +2593,7 @@ void Foam::autoLayerDriver::addLayers
"layerFaces", "layerFaces",
findIndices(flaggedCells, true) findIndices(flaggedCells, true)
); );
layerFacesSet.instance() = meshRefiner_.timeName();
Info<< "Writing " Info<< "Writing "
<< returnReduce(layerFacesSet.size(), sumOp<label>()) << returnReduce(layerFacesSet.size(), sumOp<label>())
<< " faces inside added layer to faceSet " << " faces inside added layer to faceSet "
@ -3279,6 +2670,9 @@ void Foam::autoLayerDriver::addLayers
mesh.clearOut(); mesh.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh.setInstance(meshRefiner_.timeName());
meshRefiner_.updateMesh(map, labelList(0)); meshRefiner_.updateMesh(map, labelList(0));
@ -3345,6 +2739,7 @@ void Foam::autoLayerDriver::addLayers
// ~~~~~~~~~~ // ~~~~~~~~~~
cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true)); cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
addedCellSet.instance() = meshRefiner_.timeName();
Info<< "Writing " Info<< "Writing "
<< returnReduce(addedCellSet.size(), sumOp<label>()) << returnReduce(addedCellSet.size(), sumOp<label>())
<< " added cells to cellSet " << " added cells to cellSet "
@ -3352,6 +2747,7 @@ void Foam::autoLayerDriver::addLayers
addedCellSet.write(); addedCellSet.write();
faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true)); faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
layerFacesSet.instance() = meshRefiner_.timeName();
Info<< "Writing " Info<< "Writing "
<< returnReduce(layerFacesSet.size(), sumOp<label>()) << returnReduce(layerFacesSet.size(), sumOp<label>())
<< " faces inside added layer to faceSet " << " faces inside added layer to faceSet "

40
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
View File

@ -105,46 +105,6 @@ class autoLayerDriver
// Private Member Functions // Private Member Functions
// Face merging
//- Merge patch faces. Undo until no checkMesh errors.
label mergePatchFacesUndo
(
const scalar minCos,
const scalar concaveCos,
const dictionary&
);
//- Remove points.
autoPtr<mapPolyMesh> doRemovePoints
(
removePoints& pointRemover,
const boolList& pointCanBeDeleted
);
//- Restore faces (which contain removed points)
autoPtr<mapPolyMesh> doRestorePoints
(
removePoints& pointRemover,
const labelList& facesToRestore
);
//- Return candidateFaces that are also in set.
labelList collectFaces
(
const labelList& candidateFaces,
const labelHashSet& set
) const;
//- Pick up faces of cells of faces in set.
labelList growFaceCellFace(const labelHashSet&) const;
//- Remove points not used by any face or points used by only
// two faces where the edges are in line
label mergeEdgesUndo(const scalar minCos, const dictionary&);
// Layers // Layers
//- For debugging: Dump displacement to .obj files //- For debugging: Dump displacement to .obj files

101
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -30,8 +30,8 @@ License
#include "cellSet.H" #include "cellSet.H"
#include "syncTools.H" #include "syncTools.H"
#include "refinementParameters.H" #include "refinementParameters.H"
#include "featureEdgeMesh.H"
#include "refinementSurfaces.H" #include "refinementSurfaces.H"
#include "refinementFeatures.H"
#include "shellSurfaces.H" #include "shellSurfaces.H"
#include "mapDistributePolyMesh.H" #include "mapDistributePolyMesh.H"
#include "unitConversion.H" #include "unitConversion.H"
@ -48,62 +48,6 @@ defineTypeNameAndDebug(autoRefineDriver, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Read explicit feature edges
Foam::label Foam::autoRefineDriver::readFeatureEdges
(
const PtrList<dictionary>& featDicts,
PtrList<featureEdgeMesh>& featureMeshes,
labelList& featureLevels
) const
{
Info<< "Reading external feature lines." << endl;
const fvMesh& mesh = meshRefiner_.mesh();
featureMeshes.setSize(featDicts.size());
featureLevels.setSize(featDicts.size());
forAll(featDicts, i)
{
const dictionary& dict = featDicts[i];
fileName featFileName(dict.lookup("file"));
featureMeshes.set
(
i,
new featureEdgeMesh
(
IOobject
(
featFileName, // name
//mesh.time().findInstance("triSurface", featFileName),
// // instance
mesh.time().constant(), // instance
"triSurface", // local
mesh.time(), // registry
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
)
);
featureMeshes[i].mergePoints(meshRefiner_.mergeDistance());
featureLevels[i] = readLabel(dict.lookup("level"));
Info<< "Refinement level " << featureLevels[i]
<< " for all cells crossed by feature " << featFileName
<< " (" << featureMeshes[i].points().size() << " points, "
<< featureMeshes[i].edges().size() << " edges)." << endl;
}
Info<< "Read feature lines in = "
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
return featureMeshes.size();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -128,23 +72,15 @@ Foam::autoRefineDriver::autoRefineDriver
Foam::label Foam::autoRefineDriver::featureEdgeRefine Foam::label Foam::autoRefineDriver::featureEdgeRefine
( (
const refinementParameters& refineParams, const refinementParameters& refineParams,
const PtrList<dictionary>& featDicts,
const label maxIter, const label maxIter,
const label minRefine const label minRefine
) )
{ {
const fvMesh& mesh = meshRefiner_.mesh(); const fvMesh& mesh = meshRefiner_.mesh();
// Read explicit feature edges
PtrList<featureEdgeMesh> featureMeshes;
// Per feature the refinement level
labelList featureLevels;
readFeatureEdges(featDicts, featureMeshes, featureLevels);
label iter = 0; label iter = 0;
if (featureMeshes.size() && maxIter > 0) if (meshRefiner_.features().size() && maxIter > 0)
{ {
for (; iter < maxIter; iter++) for (; iter < maxIter; iter++)
{ {
@ -160,9 +96,6 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
refineParams.keepPoints()[0], // For now only use one. refineParams.keepPoints()[0], // For now only use one.
refineParams.curvature(), refineParams.curvature(),
featureMeshes,
featureLevels,
true, // featureRefinement true, // featureRefinement
false, // internalRefinement false, // internalRefinement
false, // surfaceRefinement false, // surfaceRefinement
@ -266,8 +199,6 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
// Only look at surface intersections (minLevel and surface curvature), // Only look at surface intersections (minLevel and surface curvature),
// do not do internal refinement (refinementShells) // do not do internal refinement (refinementShells)
const PtrList<featureEdgeMesh> dummyFeatures;
labelList candidateCells labelList candidateCells
( (
meshRefiner_.refineCandidates meshRefiner_.refineCandidates
@ -275,9 +206,6 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
refineParams.keepPoints()[0], refineParams.keepPoints()[0],
refineParams.curvature(), refineParams.curvature(),
dummyFeatures, // dummy featureMeshes;
labelList(0), // dummy featureLevels;
false, // featureRefinement false, // featureRefinement
false, // internalRefinement false, // internalRefinement
true, // surfaceRefinement true, // surfaceRefinement
@ -432,8 +360,6 @@ Foam::label Foam::autoRefineDriver::shellRefine
<< "----------------------------" << nl << "----------------------------" << nl
<< endl; << endl;
const PtrList<featureEdgeMesh> dummyFeatures;
labelList candidateCells labelList candidateCells
( (
meshRefiner_.refineCandidates meshRefiner_.refineCandidates
@ -441,9 +367,6 @@ Foam::label Foam::autoRefineDriver::shellRefine
refineParams.keepPoints()[0], refineParams.keepPoints()[0],
refineParams.curvature(), refineParams.curvature(),
dummyFeatures, // dummy featureMeshes;
labelList(0), // dummy featureLevels;
false, // featureRefinement false, // featureRefinement
true, // internalRefinement true, // internalRefinement
false, // surfaceRefinement false, // surfaceRefinement
@ -459,7 +382,7 @@ Foam::label Foam::autoRefineDriver::shellRefine
<< " cells to cellSet candidateCellsFromShells." << endl; << " cells to cellSet candidateCellsFromShells." << endl;
cellSet c(mesh, "candidateCellsFromShells", candidateCells); cellSet c(mesh, "candidateCellsFromShells", candidateCells);
c.instance() = mesh.time().timeName(); c.instance() = meshRefiner_.timeName();
c.write(); c.write();
} }
@ -677,6 +600,7 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
handleSnapProblems, // remove perp edge connected cells handleSnapProblems, // remove perp edge connected cells
perpAngle, // perp angle perpAngle, // perp angle
false, // merge free standing baffles? false, // merge free standing baffles?
//true, // merge free standing baffles?
motionDict, motionDict,
const_cast<Time&>(mesh.time()), const_cast<Time&>(mesh.time()),
globalToPatch_, globalToPatch_,
@ -735,7 +659,8 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
void Foam::autoRefineDriver::mergePatchFaces void Foam::autoRefineDriver::mergePatchFaces
( (
const refinementParameters& refineParams const refinementParameters& refineParams,
const dictionary& motionDict
) )
{ {
const fvMesh& mesh = meshRefiner_.mesh(); const fvMesh& mesh = meshRefiner_.mesh();
@ -750,11 +675,12 @@ void Foam::autoRefineDriver::mergePatchFaces
const_cast<Time&>(mesh.time())++; const_cast<Time&>(mesh.time())++;
} }
meshRefiner_.mergePatchFaces meshRefiner_.mergePatchFacesUndo
( (
Foam::cos(degToRad(45.0)), Foam::cos(degToRad(45.0)),
Foam::cos(degToRad(45.0)), Foam::cos(degToRad(45.0)),
meshRefiner_.meshedPatches() meshRefiner_.meshedPatches(),
motionDict
); );
if (debug) if (debug)
@ -762,7 +688,7 @@ void Foam::autoRefineDriver::mergePatchFaces
meshRefiner_.checkData(); meshRefiner_.checkData();
} }
meshRefiner_.mergeEdges(Foam::cos(degToRad(45.0))); meshRefiner_.mergeEdgesUndo(Foam::cos(degToRad(45.0)), motionDict);
if (debug) if (debug)
{ {
@ -789,13 +715,10 @@ void Foam::autoRefineDriver::doRefine
// Check that all the keep points are inside the mesh. // Check that all the keep points are inside the mesh.
refineParams.findCells(mesh); refineParams.findCells(mesh);
PtrList<dictionary> featDicts(refineDict.lookup("features"));
// Refine around feature edges // Refine around feature edges
featureEdgeRefine featureEdgeRefine
( (
refineParams, refineParams,
featDicts,
100, // maxIter 100, // maxIter
0 // min cells to refine 0 // min cells to refine
); );
@ -833,7 +756,7 @@ void Foam::autoRefineDriver::doRefine
// Do something about cells with refined faces on the boundary // Do something about cells with refined faces on the boundary
if (prepareForSnapping) if (prepareForSnapping)
{ {
mergePatchFaces(refineParams); mergePatchFaces(refineParams, motionDict);
} }

15
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -42,7 +42,6 @@ namespace Foam
{ {
// Forward declaration of classes // Forward declaration of classes
class featureEdgeMesh;
class refinementParameters; class refinementParameters;
class meshRefinement; class meshRefinement;
class decompositionMethod; class decompositionMethod;
@ -71,19 +70,10 @@ class autoRefineDriver
// Private Member Functions // Private Member Functions
//- Read explicit feature edges
label readFeatureEdges
(
const PtrList<dictionary>& featDicts,
PtrList<featureEdgeMesh>& featureMeshes,
labelList& featureLevel
) const;
//- Refine all cells pierced by explicit feature edges //- Refine all cells pierced by explicit feature edges
label featureEdgeRefine label featureEdgeRefine
( (
const refinementParameters& refineParams, const refinementParameters& refineParams,
const PtrList<dictionary>& featDicts,
const label maxIter, const label maxIter,
const label minRefine const label minRefine
); );
@ -130,7 +120,8 @@ class autoRefineDriver
//- Merge refined boundary faces (from exposing coarser cell) //- Merge refined boundary faces (from exposing coarser cell)
void mergePatchFaces void mergePatchFaces
( (
const refinementParameters& refineParams const refinementParameters& refineParams,
const dictionary& motionDict
); );

182
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -27,7 +27,6 @@ Description
\*----------------------------------------------------------------------------*/ \*----------------------------------------------------------------------------*/
#include "autoSnapDriver.H" #include "autoSnapDriver.H"
#include "Time.H"
#include "motionSmoother.H" #include "motionSmoother.H"
#include "polyTopoChange.H" #include "polyTopoChange.H"
#include "OFstream.H" #include "OFstream.H"
@ -36,13 +35,12 @@ Description
#include "Time.H" #include "Time.H"
#include "OFstream.H" #include "OFstream.H"
#include "mapPolyMesh.H" #include "mapPolyMesh.H"
#include "motionSmoother.H"
#include "pointEdgePoint.H" #include "pointEdgePoint.H"
#include "PointEdgeWave.H" #include "PointEdgeWave.H"
#include "mergePoints.H" #include "mergePoints.H"
#include "snapParameters.H" #include "snapParameters.H"
#include "refinementSurfaces.H" #include "refinementSurfaces.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -509,9 +507,7 @@ void Foam::autoSnapDriver::dumpMove
) )
{ {
// Dump direction of growth into file // Dump direction of growth into file
Pout<< nl << "Dumping move direction to " << fName << nl Pout<< nl << "Dumping move direction to " << fName << endl;
<< "View this Lightwave-OBJ file with e.g. javaview" << nl
<< endl;
OFstream nearestStream(fName); OFstream nearestStream(fName);
@ -652,8 +648,7 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
-GREAT // null value -GREAT // null value
); );
tmp<scalarField> tfld = snapParams.snapTol()*maxEdgeLen; return snapParams.snapTol()*maxEdgeLen;
return tfld();
} }
@ -688,6 +683,7 @@ void Foam::autoSnapDriver::preSmoothPatch
// The current mesh is the starting mesh to smooth from. // The current mesh is the starting mesh to smooth from.
meshMover.setDisplacement(patchDisp); meshMover.setDisplacement(patchDisp);
meshMover.correct(); meshMover.correct();
scalar oldErrorReduction = -1; scalar oldErrorReduction = -1;
@ -828,7 +824,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
( (
unzonedSurfaces, unzonedSurfaces,
localPoints, localPoints,
sqr(4*snapDist), // sqr of attract distance sqr(snapDist), // sqr of attract distance
hitSurface, hitSurface,
hitInfo hitInfo
); );
@ -880,7 +876,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
( (
labelList(1, zoneSurfI), labelList(1, zoneSurfI),
pointField(localPoints, zonePointIndices), pointField(localPoints, zonePointIndices),
sqr(4*scalarField(minSnapDist, zonePointIndices)), sqr(scalarField(minSnapDist, zonePointIndices)),
hitSurface, hitSurface,
hitInfo hitInfo
); );
@ -956,8 +952,8 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
mesh, mesh,
pp.meshPoints(), pp.meshPoints(),
patchDisp, patchDisp,
minMagEqOp(), // combine op minMagSqrEqOp<point>(), // combine op
vector(GREAT, GREAT, GREAT) // null value vector(GREAT, GREAT, GREAT) // null value (note: cannot use VGREAT)
); );
@ -990,13 +986,12 @@ void Foam::autoSnapDriver::smoothDisplacement
) const ) const
{ {
const fvMesh& mesh = meshRefiner_.mesh(); const fvMesh& mesh = meshRefiner_.mesh();
const pointMesh& pMesh = meshMover.pMesh();
const indirectPrimitivePatch& pp = meshMover.patch(); const indirectPrimitivePatch& pp = meshMover.patch();
Info<< "Smoothing displacement ..." << endl; Info<< "Smoothing displacement ..." << endl;
// Set edge diffusivity as inverse of distance to patch // Set edge diffusivity as inverse of distance to patch
scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL)); scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + SMALL));
//scalarField edgeGamma(mesh.nEdges(), 1.0); //scalarField edgeGamma(mesh.nEdges(), 1.0);
//scalarField edgeGamma(wallGamma(mesh, pp, 10, 1)); //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
@ -1010,7 +1005,6 @@ void Foam::autoSnapDriver::smoothDisplacement
Info<< "Iteration " << iter << endl; Info<< "Iteration " << iter << endl;
} }
pointVectorField oldDisp(disp); pointVectorField oldDisp(disp);
meshMover.smooth(oldDisp, edgeGamma, disp); meshMover.smooth(oldDisp, edgeGamma, disp);
} }
Info<< "Displacement smoothed in = " Info<< "Displacement smoothed in = "
@ -1045,7 +1039,7 @@ void Foam::autoSnapDriver::smoothDisplacement
} }
void Foam::autoSnapDriver::scaleMesh bool Foam::autoSnapDriver::scaleMesh
( (
const snapParameters& snapParams, const snapParameters& snapParams,
const label nInitErrors, const label nInitErrors,
@ -1061,6 +1055,8 @@ void Foam::autoSnapDriver::scaleMesh
scalar oldErrorReduction = -1; scalar oldErrorReduction = -1;
bool meshOk = false;
Info<< "Moving mesh ..." << endl; Info<< "Moving mesh ..." << endl;
for (label iter = 0; iter < 2*snapParams.nSnap(); iter++) for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
{ {
@ -1072,10 +1068,11 @@ void Foam::autoSnapDriver::scaleMesh
oldErrorReduction = meshMover.setErrorReduction(0.0); oldErrorReduction = meshMover.setErrorReduction(0.0);
} }
if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors)) meshOk = meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors);
if (meshOk)
{ {
Info<< "Successfully moved mesh" << endl; Info<< "Successfully moved mesh" << endl;
break; break;
} }
if (debug) if (debug)
@ -1098,6 +1095,8 @@ void Foam::autoSnapDriver::scaleMesh
} }
Info<< "Moved mesh in = " Info<< "Moved mesh in = "
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl; << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
return meshOk;
} }
@ -1191,7 +1190,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
( (
unzonedSurfaces, unzonedSurfaces,
localFaceCentres, localFaceCentres,
sqr(4*faceSnapDist), // sqr of attract distance sqr(faceSnapDist), // sqr of attract distance
hitSurface, hitSurface,
hitRegion hitRegion
); );
@ -1260,6 +1259,7 @@ void Foam::autoSnapDriver::doSnap
( (
const dictionary& snapDict, const dictionary& snapDict,
const dictionary& motionDict, const dictionary& motionDict,
const scalar featureCos,
const snapParameters& snapParams const snapParameters& snapParams
) )
{ {
@ -1278,6 +1278,21 @@ void Foam::autoSnapDriver::doSnap
List<labelPair> baffles; List<labelPair> baffles;
meshRefiner_.createZoneBaffles(globalToPatch_, baffles); meshRefiner_.createZoneBaffles(globalToPatch_, baffles);
bool doFeatures = false;
label nFeatIter = 1;
if (snapParams.nFeatureSnap() > 0)
{
doFeatures = true;
nFeatIter = snapParams.nFeatureSnap();
Info<< "Snapping to features in " << nFeatIter
<< " iterations ..." << endl;
}
bool meshOk = false;
{ {
autoPtr<indirectPrimitivePatch> ppPtr autoPtr<indirectPrimitivePatch> ppPtr
( (
@ -1287,10 +1302,9 @@ void Foam::autoSnapDriver::doSnap
adaptPatchIDs adaptPatchIDs
) )
); );
indirectPrimitivePatch& pp = ppPtr();
// Distance to attract to nearest feature on surface // Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(snapParams, pp)); const scalarField snapDist(calcSnapDistance(snapParams, ppPtr()));
// Construct iterative mesh mover. // Construct iterative mesh mover.
@ -1302,7 +1316,7 @@ void Foam::autoSnapDriver::doSnap
motionSmoother meshMover motionSmoother meshMover
( (
mesh, mesh,
pp, ppPtr(),
adaptPatchIDs, adaptPatchIDs,
meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs), meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
motionDict motionDict
@ -1330,19 +1344,97 @@ void Foam::autoSnapDriver::doSnap
// Pre-smooth patch vertices (so before determining nearest) // Pre-smooth patch vertices (so before determining nearest)
preSmoothPatch(snapParams, nInitErrors, baffles, meshMover); preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
if (debug)
{
Pout<< "Writing patch smoothed mesh to time "
<< meshRefiner_.timeName() << '.' << endl;
meshRefiner_.write
(
debug,
mesh.time().path()/meshRefiner_.timeName()
);
Pout<< "Dumped mesh in = "
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
}
for (label iter = 0; iter < nFeatIter; iter++)
{
Info<< nl
<< "Morph iteration " << iter << nl
<< "-----------------" << endl;
// Calculate displacement at every patch point. Insert into // Calculate displacement at every patch point. Insert into
// meshMover. // meshMover.
calcNearestSurface(snapDist, meshMover); vectorField disp = calcNearestSurface(snapDist, meshMover);
//// Get smoothly varying internal displacement field. // Override displacement with feature edge attempt
//- 2009-12-16 : was not found to be beneficial. Keeping internal if (doFeatures)
// fields fixed slightly increases skewness (on boundary) {
// but lowers non-orthogonality quite a bit (e.g. 65->59 degrees). disp = calcNearestSurfaceFeature
// Maybe if better smoother? (
//smoothDisplacement(snapParams, meshMover); iter,
featureCos,
scalar(iter+1)/nFeatIter,
snapDist,
disp,
meshMover
);
}
if (debug&meshRefinement::OBJINTERSECTIONS)
{
dumpMove
(
mesh.time().path()
/ "patchDisplacement_" + name(iter) + ".obj",
ppPtr().localPoints(),
ppPtr().localPoints() + disp
);
}
// Get smoothly varying internal displacement field.
smoothDisplacement(snapParams, meshMover);
// Apply internal displacement to mesh. // Apply internal displacement to mesh.
scaleMesh(snapParams, nInitErrors, baffles, meshMover); meshOk = scaleMesh
(
snapParams,
nInitErrors,
baffles,
meshMover
);
if (!meshOk)
{
Info<< "Did not succesfully snap mesh. Giving up."
<< nl << endl;
// Use current mesh as base mesh
meshMover.correct();
break;
}
if (debug)
{
const_cast<Time&>(mesh.time())++;
Pout<< "Writing scaled mesh to time "
<< meshRefiner_.timeName() << endl;
meshRefiner_.write
(
debug,
mesh.time().path()/meshRefiner_.timeName()
);
Pout<< "Writing displacement field ..." << endl;
meshMover.displacement().write();
tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
magDisp().write();
}
// Use current mesh as base mesh
meshMover.correct();
}
} }
// Merge any introduced baffles. // Merge any introduced baffles.
@ -1350,6 +1442,34 @@ void Foam::autoSnapDriver::doSnap
// Repatch faces according to nearest. // Repatch faces according to nearest.
repatchToSurface(snapParams, adaptPatchIDs); repatchToSurface(snapParams, adaptPatchIDs);
// Repatching might have caused faces to be on same patch and hence
// mergeable so try again to merge coplanar faces
label nChanged = meshRefiner_.mergePatchFacesUndo
(
featureCos, // minCos
featureCos, // concaveCos
meshRefiner_.meshedPatches(),
motionDict
);
nChanged += meshRefiner_.mergeEdgesUndo
(
featureCos,
motionDict
);
if (nChanged > 0 && debug)
{
const_cast<Time&>(mesh.time())++;
Pout<< "Writing patchFace merged mesh to time "
<< meshRefiner_.timeName() << endl;
meshRefiner_.write
(
debug,
mesh.time().path()/meshRefiner_.timeName()
);
}
} }

197
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -29,6 +29,7 @@ Description
SourceFiles SourceFiles
autoSnapDriver.C autoSnapDriver.C
autoSnapDriverFeature.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
@ -45,6 +46,7 @@ namespace Foam
// Forward declaration of classes // Forward declaration of classes
class motionSmoother; class motionSmoother;
class snapParameters; class snapParameters;
class pointConstraint;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class autoSnapDriver Declaration Class autoSnapDriver Declaration
@ -52,23 +54,6 @@ class snapParameters;
class autoSnapDriver class autoSnapDriver
{ {
// Private classes
//- Combine operator class for equalizing displacements.
class minMagEqOp
{
public:
void operator()(vector& x, const vector& y) const
{
if (magSqr(y) < magSqr(x))
{
x = y;
}
}
};
// Private data // Private data
//- Mesh+surface //- Mesh+surface
@ -124,6 +109,174 @@ class autoSnapDriver
const vectorField& const vectorField&
); );
// Feature line snapping
void smoothAndConstrain
(
const indirectPrimitivePatch& pp,
const List<pointConstraint>& constraints,
vectorField& disp
) const;
void calcNearest
(
const pointField& points,
vectorField& disp,
vectorField& surfaceNormal
) const;
void calcNearestFace
(
const label iter,
const indirectPrimitivePatch& pp,
vectorField& faceDisp,
vectorField& faceSurfaceNormal,
vectorField& faceRotation
) const;
void interpolateFaceToPoint
(
const label iter,
const indirectPrimitivePatch& pp,
const vectorField& faceSurfaceNormal,
const vectorField& faceDisp,
const vectorField& faceRotation,
vectorField& patchDisp,
vectorField& patchRotationDisp
) const;
void correctAttraction
(
const DynamicList<point>& surfacePoints,
const DynamicList<label>& surfaceCount,
const point& edgePt,
const vector& edgeNormal, // normalised normal
const point& pt,
vector& edgeOffset // offset from pt to point on edge
) const;
void binFeatureFace
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch& pp,
const scalar snapDist,
const point& fc,
const vector& faceSurfaceNormal,
const vector& faceDisp,
DynamicList<point>& surfacePoints,
DynamicList<vector>& surfaceNormals,
DynamicList<label>& surfaceCount
) const;
void binFeatureFaces
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const label pointI,
const List<List<point> >& pointFaceNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
DynamicList<point>& surfacePoints,
DynamicList<vector>& surfaceNormals,
DynamicList<label>& surfaceCount
) const;
void featureAttractionUsingReconstruction
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const List<List<point> >& pointFaceNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
void determineAllFeatures
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch&,
const scalarField&,
const List<List<point> >& pointFaceNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
List<labelList>& pointAttractor,
List<List<pointConstraint> >& pointConstraints,
// Feature-edge to pp point
List<List<DynamicList<point> > >& edgeAttractors,
List<List<DynamicList<pointConstraint> > >& edgeConstraints,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
void determineFeatures
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch&,
const scalarField&,
const List<List<point> >& pointFaceNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
List<labelList>& pointAttractor,
List<List<pointConstraint> >& pointConstraints,
// Feature-edge to pp point
List<List<DynamicList<point> > >& edgeAttractors,
List<List<DynamicList<pointConstraint> > >& edgeConstraints,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
void featureAttractionUsingFeatureEdges
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const List<List<point> >& pointFaceNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
void preventFaceSqueeze
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
vectorField calcNearestSurfaceFeature
(
const label iter,
const scalar featureCos,
const scalar featureAttract,
const scalarField& snapDist,
const vectorField& nearestDisp,
motionSmoother& meshMover
) const;
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
autoSnapDriver(const autoSnapDriver&); autoSnapDriver(const autoSnapDriver&);
@ -188,6 +341,7 @@ public:
motionSmoother& meshMover motionSmoother& meshMover
) const; ) const;
//- Smooth the displacement field to the internal. //- Smooth the displacement field to the internal.
void smoothDisplacement void smoothDisplacement
( (
@ -196,8 +350,9 @@ public:
) const; ) const;
//- Do the hard work: move the mesh according to displacement, //- Do the hard work: move the mesh according to displacement,
// locally relax the displacement. // locally relax the displacement. Return true if ended up with
void scaleMesh // correct mesh, false if not.
bool scaleMesh
( (
const snapParameters& snapParams, const snapParameters& snapParams,
const label nInitErrors, const label nInitErrors,
@ -212,11 +367,11 @@ public:
const labelList& adaptPatchIDs const labelList& adaptPatchIDs
); );
void doSnap void doSnap
( (
const dictionary& snapDict, const dictionary& snapDict,
const dictionary& motionDict, const dictionary& motionDict,
const scalar featureCos,
const snapParameters& snapParams const snapParameters& snapParams
); );

2146
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C Normal file
View File

File diff suppressed because it is too large Load Diff

15
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -27,23 +27,14 @@ License
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from dictionary
Foam::snapParameters::snapParameters(const dictionary& dict, const label dummy)
:
nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
snapTol_(readScalar(dict.lookup("snapTol"))),
nSmoothDispl_(readLabel(dict.lookup("nSmoothDispl"))),
nSnap_(readLabel(dict.lookup("nSnap")))
{}
// Construct from dictionary // Construct from dictionary
Foam::snapParameters::snapParameters(const dictionary& dict) Foam::snapParameters::snapParameters(const dictionary& dict)
: :
nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))), nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
snapTol_(readScalar(dict.lookup("tolerance"))), snapTol_(readScalar(dict.lookup("tolerance"))),
nSmoothDispl_(readLabel(dict.lookup("nSolveIter"))), nSmoothDispl_(readLabel(dict.lookup("nSolveIter"))),
nSnap_(readLabel(dict.lookup("nRelaxIter"))) nSnap_(readLabel(dict.lookup("nRelaxIter"))),
nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1))
{} {}

15
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -61,6 +61,8 @@ class snapParameters
const label nSnap_; const label nSnap_;
const label nFeatureSnap_;
// Private Member Functions // Private Member Functions
@ -75,10 +77,7 @@ public:
// Constructors // Constructors
//- Construct from dictionary - old syntax //- Construct from dictionary
snapParameters(const dictionary& dict, const label dummy);
//- Construct from dictionary - new syntax
snapParameters(const dictionary& dict); snapParameters(const dictionary& dict);
@ -114,6 +113,12 @@ public:
{ {
return nSnap_; return nSnap_;
} }
label nFeatureSnap() const
{
return nFeatureSnap_;
}
}; };

55
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
View File

@ -31,6 +31,7 @@ License
#include "Time.H" #include "Time.H"
#include "refinementHistory.H" #include "refinementHistory.H"
#include "refinementSurfaces.H" #include "refinementSurfaces.H"
#include "refinementFeatures.H"
#include "decompositionMethod.H" #include "decompositionMethod.H"
#include "regionSplit.H" #include "regionSplit.H"
#include "fvMeshDistribute.H" #include "fvMeshDistribute.H"
@ -497,6 +498,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
mesh_.clearOut(); mesh_.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
setInstance(mesh_.facesInstance());
// Update local mesh data // Update local mesh data
cellRemover.updateMesh(map); cellRemover.updateMesh(map);
@ -857,6 +862,7 @@ Foam::meshRefinement::meshRefinement
const scalar mergeDistance, const scalar mergeDistance,
const bool overwrite, const bool overwrite,
const refinementSurfaces& surfaces, const refinementSurfaces& surfaces,
const refinementFeatures& features,
const shellSurfaces& shells const shellSurfaces& shells
) )
: :
@ -865,6 +871,7 @@ Foam::meshRefinement::meshRefinement
overwrite_(overwrite), overwrite_(overwrite),
oldInstance_(mesh.pointsInstance()), oldInstance_(mesh.pointsInstance()),
surfaces_(surfaces), surfaces_(surfaces),
features_(features),
shells_(shells), shells_(shells),
meshCutter_ meshCutter_
( (
@ -1271,6 +1278,10 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
// Update numbering of meshRefiner // Update numbering of meshRefiner
distribute(map); distribute(map);
// Set correct instance (for if overwrite)
mesh_.setInstance(timeName());
setInstance(mesh_.facesInstance());
} }
return map; return map;
} }
@ -1449,6 +1460,8 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
} }
} }
// Note: time().timeName() instead of meshRefinement::timeName() since
// postprocessable field.
tmp<pointVectorField> tfld tmp<pointVectorField> tfld
( (
new pointVectorField new pointVectorField
@ -1456,7 +1469,7 @@ Foam::tmp<Foam::pointVectorField> Foam::meshRefinement::makeDisplacementField
IOobject IOobject
( (
"pointDisplacement", "pointDisplacement",
mesh.time().timeName(), //timeName(), mesh.time().timeName(),
mesh, mesh,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE IOobject::AUTO_WRITE
@ -1930,10 +1943,6 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
pointMap.clear(); pointMap.clear();
} }
} }
// If necessary reset the instance
mesh_.setInstance(timeName());
setInstance(mesh_.facesInstance());
} }
@ -2059,10 +2068,6 @@ void Foam::meshRefinement::updateMesh
data.transfer(newFaceData); data.transfer(newFaceData);
} }
} }
// If necessary reset the instance
mesh_.setInstance(timeName());
setInstance(mesh_.facesInstance());
} }
@ -2118,12 +2123,32 @@ void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
<< " points(local):" << mesh_.nPoints() << " points(local):" << mesh_.nPoints()
<< endl; << endl;
} }
else
{ {
PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
label nMasterFaces = 0;
forAll(isMasterFace, i)
{
if (isMasterFace[i])
{
nMasterFaces++;
}
}
PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh_));
label nMasterPoints = 0;
forAll(isMasterPoint, i)
{
if (isMasterPoint[i])
{
nMasterPoints++;
}
}
Info<< msg.c_str() Info<< msg.c_str()
<< " : cells:" << pData.nTotalCells() << " : cells:" << pData.nTotalCells()
<< " faces:" << pData.nTotalFaces() << " faces:" << returnReduce(nMasterFaces, sumOp<label>())
<< " points:" << pData.nTotalPoints() << " points:" << returnReduce(nMasterPoints, sumOp<label>())
<< endl; << endl;
} }
@ -2168,12 +2193,14 @@ Foam::word Foam::meshRefinement::timeName() const
void Foam::meshRefinement::dumpRefinementLevel() const void Foam::meshRefinement::dumpRefinementLevel() const
{ {
// Note: use time().timeName(), not meshRefinement::timeName()
// so as to dump the fields to 0, not to constant.
volScalarField volRefLevel volScalarField volRefLevel
( (
IOobject IOobject
( (
"cellLevel", "cellLevel",
mesh_.time().timeName(),// Dump to current time, not to mesh inst mesh_.time().timeName(),
mesh_, mesh_,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::AUTO_WRITE, IOobject::AUTO_WRITE,
@ -2201,7 +2228,7 @@ void Foam::meshRefinement::dumpRefinementLevel() const
IOobject IOobject
( (
"pointLevel", "pointLevel",
mesh_.time().timeName(),// Dump to current time, not to mesh inst mesh_.time().timeName(),
mesh_, mesh_,
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE, IOobject::NO_WRITE,

64
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
View File

@ -62,13 +62,14 @@ class fvMesh;
class mapDistributePolyMesh; class mapDistributePolyMesh;
class decompositionMethod; class decompositionMethod;
class refinementSurfaces; class refinementSurfaces;
class refinementFeatures;
class shellSurfaces; class shellSurfaces;
class removeCells; class removeCells;
class featureEdgeMesh;
class fvMeshDistribute; class fvMeshDistribute;
class searchableSurface; class searchableSurface;
class regionSplit; class regionSplit;
class globalIndex; class globalIndex;
class removePoints;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class meshRefinement Declaration Class meshRefinement Declaration
@ -118,6 +119,9 @@ private:
//- All surface-intersection interaction //- All surface-intersection interaction
const refinementSurfaces& surfaces_; const refinementSurfaces& surfaces_;
//- All feature-edge interaction
const refinementFeatures& features_;
//- All shell-refinement interaction //- All shell-refinement interaction
const shellSurfaces& shells_; const shellSurfaces& shells_;
@ -260,8 +264,6 @@ private:
label markFeatureRefinement label markFeatureRefinement
( (
const point& keepPoint, const point& keepPoint,
const PtrList<featureEdgeMesh>& featureMeshes,
const labelList& featureLevels,
const label nAllowRefine, const label nAllowRefine,
labelList& refineCell, labelList& refineCell,
@ -415,12 +417,6 @@ private:
const labelList& globalToPatch const labelList& globalToPatch
) const; ) const;
////- Initial test of marking faces using geometric information.
//labelList markFacesOnProblemCellsGeometric
//(
// const dictionary& motionDict
//) const;
// Baffle merging // Baffle merging
@ -501,6 +497,7 @@ public:
const scalar mergeDistance, const scalar mergeDistance,
const bool overwrite, const bool overwrite,
const refinementSurfaces&, const refinementSurfaces&,
const refinementFeatures&,
const shellSurfaces& const shellSurfaces&
); );
@ -542,6 +539,12 @@ public:
return surfaces_; return surfaces_;
} }
//- reference to feature edge mesh
const refinementFeatures& features() const
{
return features_;
}
//- reference to refinement shells (regions) //- reference to refinement shells (regions)
const shellSurfaces& shells() const const shellSurfaces& shells() const
{ {
@ -648,9 +651,6 @@ public:
const point& keepPoint, const point& keepPoint,
const scalar curvature, const scalar curvature,
const PtrList<featureEdgeMesh>& featureMeshes,
const labelList& featureLevels,
const bool featureRefinement, const bool featureRefinement,
const bool internalRefinement, const bool internalRefinement,
const bool surfaceRefinement, const bool surfaceRefinement,
@ -803,19 +803,45 @@ public:
const Map<label>& cellsToRestore const Map<label>& cellsToRestore
); );
// Merging coplanar faces and edges
//- Merge faces on the same patch (usually from exposing refinement) label mergePatchFacesUndo
// Returns global number of faces merged.
label mergePatchFaces
( (
const scalar minCos, const scalar minCos,
const scalar concaveCos, const scalar concaveCos,
const labelList& patchIDs const labelList& patchIDs,
const dictionary& motionDict
); );
//- Remove points not used by any face or points used autoPtr<mapPolyMesh> doRemovePoints
// by only two faces where the edges are in line (
autoPtr<mapPolyMesh> mergeEdges(const scalar minCos); removePoints& pointRemover,
const boolList& pointCanBeDeleted
);
autoPtr<mapPolyMesh> doRestorePoints
(
removePoints& pointRemover,
const labelList& facesToRestore
);
labelList collectFaces
(
const labelList& candidateFaces,
const labelHashSet& set
) const;
// Pick up faces of cells of faces in set.
labelList growFaceCellFace
(
const labelHashSet& set
) const;
label mergeEdgesUndo
(
const scalar minCos,
const dictionary& motionDict
);
// Debug/IO // Debug/IO

34
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
View File

@ -524,6 +524,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
mesh_.clearOut(); mesh_.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
//- Redo the intersections on the newly create baffle faces. Note that //- Redo the intersections on the newly create baffle faces. Note that
// this changes also the cell centre positions. // this changes also the cell centre positions.
faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles); faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
@ -644,7 +648,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::getDuplicateFaces
} }
Pout<< "Writing duplicate faces (baffles) to faceSet " Pout<< "Writing duplicate faces (baffles) to faceSet "
<< duplicateFaceSet.name() << nl << endl; << duplicateFaceSet.name() << nl << endl;
duplicateFaceSet.instance() = mesh_.time().timeName(); duplicateFaceSet.instance() = timeName();
duplicateFaceSet.write(); duplicateFaceSet.write();
} }
@ -724,9 +728,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
if (debug) if (debug)
{ {
const_cast<Time&>(mesh_.time())++; const_cast<Time&>(mesh_.time())++;
Pout<< "Writing baffled mesh to time " Pout<< "Writing zone-baffled mesh to time " << timeName()
<< mesh_.time().timeName() << endl; << endl;
mesh_.write(); write(debug, mesh_.time().path()/"baffles");
} }
} }
Info<< "Created " << nZoneFaces << " baffles in = " Info<< "Created " << nZoneFaces << " baffles in = "
@ -746,14 +750,6 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
const List<labelPair>& couples const List<labelPair>& couples
) const ) const
{ {
// Construct addressing engine for all duplicate faces (only one
// for each pair)
// Duplicate faces in mesh labels (first face of each pair only)
// (reused later on to mark off filtered couples. see below)
labelList duplicateFaces(couples.size());
// All duplicate faces on edge of the patch are to be merged. // All duplicate faces on edge of the patch are to be merged.
// So we count for all edges of duplicate faces how many duplicate // So we count for all edges of duplicate faces how many duplicate
// faces use them. // faces use them.
@ -824,7 +820,8 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
// Baffles which are not next to other boundaries and baffles will have // Baffles which are not next to other boundaries and baffles will have
// value 2*1000000+2*1 // nBafflesPerEdge value 2*1000000+2*1 (from 2 boundary faces which are
// both baffle faces)
List<labelPair> filteredCouples(couples.size()); List<labelPair> filteredCouples(couples.size());
label filterI = 0; label filterI = 0;
@ -967,6 +964,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
mesh_.clearOut(); mesh_.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
// Update intersections. Recalculate intersections on merged faces since // Update intersections. Recalculate intersections on merged faces since
// this seems to give problems? Note: should not be nessecary since // this seems to give problems? Note: should not be nessecary since
// baffles preserve intersections from when they were created. // baffles preserve intersections from when they were created.
@ -1789,7 +1789,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
problemTopo.insert(faceI); problemTopo.insert(faceI);
} }
} }
problemTopo.instance() = mesh_.time().timeName(); problemTopo.instance() = timeName();
Pout<< "Dumping " << problemTopo.size() Pout<< "Dumping " << problemTopo.size()
<< " problem faces to " << problemTopo.objectPath() << endl; << " problem faces to " << problemTopo.objectPath() << endl;
problemTopo.write(); problemTopo.write();
@ -2246,6 +2246,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
mesh_.clearOut(); mesh_.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
// Update intersections. Is mapping only (no faces created, positions stay // Update intersections. Is mapping only (no faces created, positions stay
// same) so no need to recalculate intersections. // same) so no need to recalculate intersections.
updateMesh(map, labelList(0)); updateMesh(map, labelList(0));
@ -2782,6 +2785,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
mesh_.clearOut(); mesh_.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
// Print some stats (note: zones are synchronised) // Print some stats (note: zones are synchronised)
if (mesh_.cellZones().size() > 0) if (mesh_.cellZones().size() > 0)
{ {

849
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -27,25 +27,240 @@ License
#include "combineFaces.H" #include "combineFaces.H"
#include "polyTopoChange.H" #include "polyTopoChange.H"
#include "removePoints.H" #include "removePoints.H"
#include "faceSet.H"
#include "Time.H"
#include "motionSmoother.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Merge faces that are in-line. //// Merge faces that are in-line.
Foam::label Foam::meshRefinement::mergePatchFaces //Foam::label Foam::meshRefinement::mergePatchFaces
//(
// const scalar minCos,
// const scalar concaveCos,
// const labelList& patchIDs
//)
//{
// // Patch face merging engine
// combineFaces faceCombiner(mesh_);
//
// const polyBoundaryMesh& patches = mesh_.boundaryMesh();
//
// // Pick up all candidate cells on boundary
// labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
//
// 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 mergeSets
// (
// faceCombiner.getMergeSets
// (
// minCos,
// concaveCos,
// boundaryCells
// )
// );
//
// label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>());
//
// Info<< "mergePatchFaces : Merging " << nFaceSets
// << " sets of faces." << endl;
//
// if (nFaceSets > 0)
// {
// // Topology changes container
// polyTopoChange meshMod(mesh_);
//
// // Merge all faces of a set into the first face of the set. Remove
// // unused points.
// faceCombiner.setRefinement(mergeSets, 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());
// }
// else
// {
// // Delete mesh volumes. No other way to do this?
// mesh_.clearOut();
// }
//
//
// // Reset the instance for if in overwrite mode
// mesh_.setInstance(timeName());
//
// faceCombiner.updateMesh(map);
//
// // Get the kept faces that need to be recalculated.
// // Merging two boundary faces might shift the cell centre
// // (unless the faces are absolutely planar)
// labelHashSet retestFaces(6*mergeSets.size());
//
// forAll(mergeSets, setI)
// {
// label oldMasterI = mergeSets[setI][0];
//
// label faceI = map().reverseFaceMap()[oldMasterI];
//
// // faceI is always uncoupled boundary face
// const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
//
// forAll(cFaces, i)
// {
// retestFaces.insert(cFaces[i]);
// }
// }
// updateMesh(map, retestFaces.toc());
// }
//
//
// return nFaceSets;
//}
//
//
//// Remove points not used by any face or points used by only two faces where
//// the edges are in line
//Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
//(
// const scalar minCos
//)
//{
// // Point removal analysis engine
// removePoints pointRemover(mesh_);
//
// // Count usage of points
// boolList pointCanBeDeleted;
// label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
//
// Info<< "Removing " << nRemove
// << " straight edge points." << endl;
//
// autoPtr<mapPolyMesh> map;
//
// if (nRemove > 0)
// {
// // Save my local faces that will change. These changed faces might
// // cause a shift in the cell centre which needs to be retested.
// // Have to do this before changing mesh since point will be removed.
// labelHashSet retestOldFaces(nRemove / Pstream::nProcs());
//
// {
// const faceList& faces = mesh_.faces();
//
// forAll(faces, faceI)
// {
// const face& f = faces[faceI];
//
// forAll(f, fp)
// {
// if (pointCanBeDeleted[f[fp]])
// {
// retestOldFaces.insert(faceI);
// break;
// }
// }
// }
// }
//
// // Topology changes container
// polyTopoChange meshMod(mesh_);
//
// pointRemover.setRefinement(pointCanBeDeleted, meshMod);
//
// // Change the mesh (no inflation)
// 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());
// }
// else
// {
// // Delete mesh volumes. No other way to do this?
// mesh_.clearOut();
// }
//
// // Reset the instance for if in overwrite mode
// mesh_.setInstance(timeName());
//
// pointRemover.updateMesh(map);
//
// // Get the kept faces that need to be recalculated.
// labelHashSet retestFaces(6*retestOldFaces.size());
//
// const cellList& cells = mesh_.cells();
//
// forAllConstIter(labelHashSet, retestOldFaces, iter)
// {
// label faceI = map().reverseFaceMap()[iter.key()];
//
// const cell& ownFaces = cells[mesh_.faceOwner()[faceI]];
//
// forAll(ownFaces, i)
// {
// retestFaces.insert(ownFaces[i]);
// }
//
// if (mesh_.isInternalFace(faceI))
// {
// const cell& neiFaces = cells[mesh_.faceNeighbour()[faceI]];
//
// forAll(neiFaces, i)
// {
// retestFaces.insert(neiFaces[i]);
// }
// }
// }
// updateMesh(map, retestFaces.toc());
// }
//
// return map;
//}
Foam::label Foam::meshRefinement::mergePatchFacesUndo
( (
const scalar minCos, const scalar minCos,
const scalar concaveCos, const scalar concaveCos,
const labelList& patchIDs const labelList& patchIDs,
const dictionary& motionDict
) )
{ {
// Patch face merging engine // Patch face merging engine
combineFaces faceCombiner(mesh_); combineFaces faceCombiner(mesh_, true);
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
// Pick up all candidate cells on boundary // Pick up all candidate cells on boundary
labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces()); labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
forAll(patchIDs, i) forAll(patchIDs, i)
{ {
label patchI = patchIDs[i]; label patchI = patchIDs[i];
@ -60,9 +275,10 @@ Foam::label Foam::meshRefinement::mergePatchFaces
} }
} }
} }
}
// Get all sets of faces that can be merged // Get all sets of faces that can be merged
labelListList mergeSets labelListList allFaceSets
( (
faceCombiner.getMergeSets faceCombiner.getMergeSets
( (
@ -72,19 +288,44 @@ Foam::label Foam::meshRefinement::mergePatchFaces
) )
); );
label nFaceSets = returnReduce(mergeSets.size(), sumOp<label>()); label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
Info<< "mergePatchFaces : Merging " << nFaceSets Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
<< " sets of faces." << endl;
if (nFaceSets > 0) 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.instance() = timeName();
allSets.write();
const_cast<Time&>(mesh_.time())++;
}
// Topology changes container // Topology changes container
polyTopoChange meshMod(mesh_); polyTopoChange meshMod(mesh_);
// Merge all faces of a set into the first face of the set. Remove // Merge all faces of a set into the first face of the set.
// unused points. faceCombiner.setRefinement(allFaceSets, meshMod);
faceCombiner.setRefinement(mergeSets, meshMod);
// Experimental: store data for all the points that have been deleted
storeData
(
faceCombiner.savedPointLabels(), // points to store
labelList(0), // faces to store
labelList(0) // cells to store
);
// Change the mesh (no inflation) // Change the mesh (no inflation)
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true); autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
@ -99,90 +340,169 @@ Foam::label Foam::meshRefinement::mergePatchFaces
} }
else else
{ {
// Delete mesh volumes. No other way to do this? // Delete mesh volumes.
mesh_.clearOut(); mesh_.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
faceCombiner.updateMesh(map); faceCombiner.updateMesh(map);
// Get the kept faces that need to be recalculated. updateMesh(map, labelList(0));
// Merging two boundary faces might shift the cell centre
// (unless the faces are absolutely planar)
labelHashSet retestFaces(6*mergeSets.size());
forAll(mergeSets, setI) if (debug)
{ {
label oldMasterI = mergeSets[setI][0]; Pout<< "Writing initial merged-faces mesh to time "
<< timeName() << nl << endl;
write();
}
label faceI = map().reverseFaceMap()[oldMasterI]; for (label iteration = 0; iteration < 100; iteration++)
{
Info<< nl
<< "Undo iteration " << iteration << nl
<< "----------------" << endl;
// faceI is always uncoupled boundary face
const cell& cFaces = mesh_.cells()[mesh_.faceOwner()[faceI]]; // 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)
{
errorFaces.instance() = timeName();
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) forAll(cFaces, i)
{ {
retestFaces.insert(cFaces[i]); if (errorFaces.found(cFaces[i]))
}
}
updateMesh(map, retestFaces.toc());
}
return nFaceSets;
}
// Remove points not used by any face or points used by only two faces where
// the edges are in line
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
(
const scalar minCos
)
{ {
// Point removal analysis engine mastersToRestore.append(masterFaceI);
removePoints pointRemover(mesh_);
// Count usage of points
boolList pointCanBeDeleted;
label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
Info<< "Removing " << nRemove
<< " straight edge points." << endl;
autoPtr<mapPolyMesh> map;
if (nRemove > 0)
{
// Save my local faces that will change. These changed faces might
// cause a shift in the cell centre which needs to be retested.
// Have to do this before changing mesh since point will be removed.
labelHashSet retestOldFaces(nRemove / Pstream::nProcs());
{
const faceList& faces = mesh_.faces();
forAll(faces, faceI)
{
const face& f = faces[faceI];
forAll(f, fp)
{
if (pointCanBeDeleted[f[fp]])
{
retestOldFaces.insert(faceI);
break; 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);
restoreSet.instance() = timeName();
Pout<< "Writing all " << mastersToRestore.size()
<< " masterfaces to be restored to set "
<< restoreSet.objectPath() << endl;
restoreSet.write();
}
if (nRestore == 0)
{
break;
}
// Undo
// ~~~~
if (debug)
{
const_cast<Time&>(mesh_.time())++;
}
// Topology changes container // Topology changes container
polyTopoChange meshMod(mesh_); polyTopoChange meshMod(mesh_);
pointRemover.setRefinement(pointCanBeDeleted, meshMod); // 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) // Change the mesh (no inflation)
map = meshMod.changeMesh(mesh_, false, true); autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
// Update fields // Update fields
mesh_.updateMesh(map); mesh_.updateMesh(map);
@ -194,43 +514,364 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
} }
else else
{ {
// Delete mesh volumes. No other way to do this? // Delete mesh volumes.
mesh_.clearOut(); mesh_.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
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
updateMesh
(
map,
labelList(0), // changedFaces
restoredPoints,
restoredFaces,
restoredCells
);
if (debug)
{
Pout<< "Writing merged-faces mesh to time "
<< timeName() << nl << endl;
write();
}
Info<< endl;
}
}
else
{
Info<< "No faces merged ..." << endl;
}
return nFaceSets;
}
// Remove points. pointCanBeDeleted is parallel synchronised.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::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());
}
else
{
// Delete mesh volumes.
mesh_.clearOut();
}
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
pointRemover.updateMesh(map); pointRemover.updateMesh(map);
updateMesh(map, labelList(0));
// Get the kept faces that need to be recalculated.
labelHashSet retestFaces(6*retestOldFaces.size());
const cellList& cells = mesh_.cells();
forAllConstIter(labelHashSet, retestOldFaces, iter)
{
label faceI = map().reverseFaceMap()[iter.key()];
const cell& ownFaces = cells[mesh_.faceOwner()[faceI]];
forAll(ownFaces, i)
{
retestFaces.insert(ownFaces[i]);
}
if (mesh_.isInternalFace(faceI))
{
const cell& neiFaces = cells[mesh_.faceNeighbour()[faceI]];
forAll(neiFaces, i)
{
retestFaces.insert(neiFaces[i]);
}
}
}
updateMesh(map, retestFaces.toc());
}
return map; return map;
} }
// Restore faces (which contain removed points)
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::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());
}
else
{
// Delete mesh volumes.
mesh_.clearOut();
}
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
pointRemover.updateMesh(map);
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::meshRefinement::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
);
labelList selectedFaces(findIndices(selected, true));
return selectedFaces;
}
// Pick up faces of cells of faces in set.
Foam::labelList Foam::meshRefinement::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
);
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::meshRefinement::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)
{
errorFaces.instance() = timeName();
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)
{
const_cast<Time&>(mesh_.time())++;
Pout<< "Writing merged-edges mesh to time "
<< timeName() << nl << endl;
write();
}
}
else
{
Info<< "No straight edges simplified and no points removed ..." << endl;
}
return nRemove;
}
// ************************************************************************* // // ************************************************************************* //

12
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd. \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -202,7 +202,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
if (debug) if (debug)
{ {
faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces); faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
fSet.instance() = mesh_.time().timeName(); fSet.instance() = timeName();
Pout<< "Writing " << fSet.size() Pout<< "Writing " << fSet.size()
<< " with problematic topology to faceSet " << " with problematic topology to faceSet "
<< fSet.objectPath() << endl; << fSet.objectPath() << endl;
@ -265,7 +265,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
if (debug) if (debug)
{ {
perpFaces.instance() = mesh_.time().timeName(); perpFaces.instance() = timeName();
Pout<< "Writing " << perpFaces.size() Pout<< "Writing " << perpFaces.size()
<< " faces that are perpendicular to the surface to set " << " faces that are perpendicular to the surface to set "
<< perpFaces.objectPath() << endl; << perpFaces.objectPath() << endl;
@ -427,6 +427,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
// Count of faces marked for baffling // Count of faces marked for baffling
label nBaffleFaces = 0; label nBaffleFaces = 0;
PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
// Count of faces not baffled since would not cause a collapse // Count of faces not baffled since would not cause a collapse
label nPrevented = 0; label nPrevented = 0;
@ -483,7 +484,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
if (debug) if (debug)
{ {
cellSet problemCellSet(mesh_, "problemCells", problemCells.toc()); cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
problemCellSet.instance() = mesh_.time().timeName(); problemCellSet.instance() = timeName();
Pout<< "Writing " << problemCellSet.size() Pout<< "Writing " << problemCellSet.size()
<< " cells that are edge connected to coarser cell to set " << " cells that are edge connected to coarser cell to set "
<< problemCellSet.objectPath() << endl; << problemCellSet.objectPath() << endl;
@ -918,7 +919,10 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
else else
{ {
facePatch[faceI] = getBafflePatch(facePatch, faceI); facePatch[faceI] = getBafflePatch(facePatch, faceI);
if (isMasterFace[faceI])
{
nBaffleFaces++; nBaffleFaces++;
}
// Do NOT update boundary data since this would grow // Do NOT update boundary data since this would grow
// blocked faces across gaps. // blocked faces across gaps.

56
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -28,6 +28,7 @@ License
#include "syncTools.H" #include "syncTools.H"
#include "Time.H" #include "Time.H"
#include "refinementSurfaces.H" #include "refinementSurfaces.H"
#include "refinementFeatures.H"
#include "shellSurfaces.H" #include "shellSurfaces.H"
#include "faceSet.H" #include "faceSet.H"
#include "decompositionMethod.H" #include "decompositionMethod.H"
@ -53,6 +54,11 @@ Foam::labelList Foam::meshRefinement::getChangedFaces
const polyMesh& mesh = map.mesh(); const polyMesh& mesh = map.mesh();
labelList changedFaces; labelList changedFaces;
// For reporting: number of masterFaces changed
label nMasterChanged = 0;
PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
{ {
// Mark any face on a cell which has been added or changed // Mark any face on a cell which has been added or changed
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -182,36 +188,26 @@ Foam::labelList Foam::meshRefinement::getChangedFaces
// Now we have in changedFace marked all affected faces. Pack. // Now we have in changedFace marked all affected faces. Pack.
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nChanged = 0; changedFaces = findIndices(changedFace, true);
// Count changed master faces.
nMasterChanged = 0;
forAll(changedFace, faceI) forAll(changedFace, faceI)
{ {
if (changedFace[faceI]) if (changedFace[faceI] && isMasterFace[faceI])
{ {
nChanged++; nMasterChanged++;
} }
} }
changedFaces.setSize(nChanged);
nChanged = 0;
forAll(changedFace, faceI)
{
if (changedFace[faceI])
{
changedFaces[nChanged++] = faceI;
} }
}
}
label nChangedFaces = changedFaces.size();
reduce(nChangedFaces, sumOp<label>());
if (debug) if (debug)
{ {
Pout<< "getChangedFaces : Detected " Pout<< "getChangedFaces : Detected "
<< " local:" << changedFaces.size() << " local:" << changedFaces.size()
<< " global:" << nChangedFaces << " global:" << returnReduce(nMasterChanged, sumOp<label>())
<< " changed faces out of " << mesh.globalData().nTotalFaces() << " changed faces out of " << mesh.globalData().nTotalFaces()
<< endl; << endl;
@ -252,8 +248,6 @@ bool Foam::meshRefinement::markForRefine
Foam::label Foam::meshRefinement::markFeatureRefinement Foam::label Foam::meshRefinement::markFeatureRefinement
( (
const point& keepPoint, const point& keepPoint,
const PtrList<featureEdgeMesh>& featureMeshes,
const labelList& featureLevels,
const label nAllowRefine, const label nAllowRefine,
labelList& refineCell, labelList& refineCell,
@ -287,9 +281,11 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
if (cellI != -1) if (cellI != -1)
{ {
forAll(featureMeshes, featI) forAll(features_, featI)
{ {
const featureEdgeMesh& featureMesh = featureMeshes[featI]; const featureEdgeMesh& featureMesh = features_[featI];
const label featureLevel = features_.levels()[featI];
const labelListList& pointEdges = featureMesh.pointEdges(); const labelListList& pointEdges = featureMesh.pointEdges();
forAll(pointEdges, pointI) forAll(pointEdges, pointI)
@ -315,7 +311,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
tetFaceI, tetFaceI,
tetPtI, tetPtI,
featureMesh.points()[pointI], // endpos featureMesh.points()[pointI], // endpos
featureLevels[featI], // level featureLevel, // level
featI, // featureMesh featI, // featureMesh
pointI // end point pointI // end point
) )
@ -339,11 +335,11 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
maxFeatureLevel = -1; maxFeatureLevel = -1;
// Whether edge has been visited. // Whether edge has been visited.
List<PackedBoolList> featureEdgeVisited(featureMeshes.size()); List<PackedBoolList> featureEdgeVisited(features_.size());
forAll(featureMeshes, featI) forAll(features_, featI)
{ {
featureEdgeVisited[featI].setSize(featureMeshes[featI].edges().size()); featureEdgeVisited[featI].setSize(features_[featI].edges().size());
featureEdgeVisited[featI] = 0u; featureEdgeVisited[featI] = 0u;
} }
@ -359,7 +355,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
label featI = tp.i(); label featI = tp.i();
label pointI = tp.j(); label pointI = tp.j();
const featureEdgeMesh& featureMesh = featureMeshes[featI]; const featureEdgeMesh& featureMesh = features_[featI];
const labelList& pEdges = featureMesh.pointEdges()[pointI]; const labelList& pEdges = featureMesh.pointEdges()[pointI];
// Particle now at pointI. Check connected edges to see which one // Particle now at pointI. Check connected edges to see which one
@ -1068,9 +1064,6 @@ Foam::labelList Foam::meshRefinement::refineCandidates
const point& keepPoint, const point& keepPoint,
const scalar curvature, const scalar curvature,
const PtrList<featureEdgeMesh>& featureMeshes,
const labelList& featureLevels,
const bool featureRefinement, const bool featureRefinement,
const bool internalRefinement, const bool internalRefinement,
const bool surfaceRefinement, const bool surfaceRefinement,
@ -1134,8 +1127,6 @@ Foam::labelList Foam::meshRefinement::refineCandidates
label nFeatures = markFeatureRefinement label nFeatures = markFeatureRefinement
( (
keepPoint, keepPoint,
featureMeshes,
featureLevels,
nAllowRefine, nAllowRefine,
refineCell, refineCell,
@ -1251,6 +1242,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::refine
mesh_.clearOut(); mesh_.clearOut();
} }
// Reset the instance for if in overwrite mode
mesh_.setInstance(timeName());
// Update intersection info // Update intersection info
updateMesh(map, getChangedFaces(map, cellsToRefine)); updateMesh(map, getChangedFaces(map, cellsToRefine));

233
src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.C Normal file
View File

@ -0,0 +1,233 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "refinementFeatures.H"
#include "Time.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::refinementFeatures::refinementFeatures
(
const objectRegistry& io,
const PtrList<dictionary>& featDicts
)
:
PtrList<featureEdgeMesh>(featDicts.size()),
levels_(featDicts.size()),
edgeTrees_(featDicts.size()),
pointTrees_(featDicts.size())
{
// Read features
forAll(featDicts, i)
{
const dictionary& dict = featDicts[i];
fileName featFileName(dict.lookup("file"));
set
(
i,
new featureEdgeMesh
(
IOobject
(
featFileName, // name
io.time().constant(), // instance
"triSurface", // local
io.time(), // registry
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
)
);
const featureEdgeMesh& eMesh = operator[](i);
//eMesh.mergePoints(meshRefiner_.mergeDistance());
levels_[i] = readLabel(dict.lookup("level"));
Info<< "Refinement level " << levels_[i]
<< " for all cells crossed by feature " << featFileName
<< " (" << eMesh.points().size() << " points, "
<< eMesh.edges().size() << " edges)." << endl;
}
// Search engines
forAll(*this, i)
{
const featureEdgeMesh& eMesh = operator[](i);
const pointField& points = eMesh.points();
const edgeList& edges = eMesh.edges();
// Calculate bb of all points
const treeBoundBox bb(points);
edgeTrees_.set
(
i,
new indexedOctree<treeDataEdge>
(
treeDataEdge
(
false, // do not cache bb
edges,
points,
identity(edges.size())
),
bb, // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
// Detect feature points from edges.
const labelListList& pointEdges = eMesh.pointEdges();
DynamicList<label> featurePoints;
forAll(pointEdges, pointI)
{
if (pointEdges[pointI].size() > 2)
{
featurePoints.append(pointI);
}
}
Info<< "Detected " << featurePoints.size()
<< " featurePoints out of " << points.size() << endl;
pointTrees_.set
(
i,
new indexedOctree<treeDataPoint>
(
treeDataPoint(points, featurePoints),
bb, // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::refinementFeatures::findNearestEdge
(
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearFeature,
List<pointIndexHit>& nearInfo
) const
{
nearFeature.setSize(samples.size());
nearFeature = -1;
nearInfo.setSize(samples.size());
forAll(edgeTrees_, featI)
{
const indexedOctree<treeDataEdge>& tree = edgeTrees_[featI];
forAll(samples, sampleI)
{
const point& sample = samples[sampleI];
scalar distSqr;
if (nearInfo[sampleI].hit())
{
distSqr = magSqr(nearInfo[sampleI].hitPoint()-sample);
}
else
{
distSqr = nearestDistSqr[sampleI];
}
pointIndexHit info = tree.findNearest(sample, distSqr);
if (info.hit())
{
nearInfo[sampleI] = info;
nearFeature[sampleI] = featI;
}
}
}
}
void Foam::refinementFeatures::findNearestPoint
(
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearFeature,
labelList& nearIndex
) const
{
nearFeature.setSize(samples.size());
nearFeature = -1;
nearIndex.setSize(samples.size());
nearIndex = -1;
forAll(pointTrees_, featI)
{
const indexedOctree<treeDataPoint>& tree = pointTrees_[featI];
forAll(samples, sampleI)
{
const point& sample = samples[sampleI];
scalar distSqr;
if (nearFeature[sampleI] != -1)
{
label nearFeatI = nearFeature[sampleI];
const indexedOctree<treeDataPoint>& nearTree =
pointTrees_[nearFeatI];
label featPointI =
nearTree.shapes().pointLabels()[nearIndex[sampleI]];
const point& featPt =
operator[](nearFeatI).points()[featPointI];
distSqr = magSqr(featPt-sample);
}
else
{
distSqr = nearestDistSqr[sampleI];
}
pointIndexHit info = tree.findNearest(sample, distSqr);
if (info.hit())
{
nearFeature[sampleI] = featI;
nearIndex[sampleI] = info.index();
}
}
}
}
// ************************************************************************* //

139
src/mesh/autoMesh/autoHexMesh/refinementFeatures/refinementFeatures.H Normal file
View File

@ -0,0 +1,139 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2011 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 3 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, see <http://www.gnu.org/licenses/>.
Class
Foam::refinementFeatures
Description
Encapsulates queries for features.
SourceFiles
refinementFeatures.C
\*---------------------------------------------------------------------------*/
#ifndef refinementFeatures_H
#define refinementFeatures_H
#include "featureEdgeMesh.H"
#include "indexedOctree.H"
#include "treeDataEdge.H"
#include "treeDataPoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class refinementFeatures Declaration
\*---------------------------------------------------------------------------*/
class refinementFeatures
:
public PtrList<featureEdgeMesh>
{
private:
// Private data
//- Refinement levels
labelList levels_;
//- Edge
PtrList<indexedOctree<treeDataEdge> > edgeTrees_;
//- Features points
PtrList<indexedOctree<treeDataPoint> > pointTrees_;
// Private Member Functions
public:
// Constructors
//- Construct from components
refinementFeatures
(
const objectRegistry& io,
const PtrList<dictionary>& featDicts
);
// Member Functions
// Access
const labelList& levels() const
{
return levels_;
}
const PtrList<indexedOctree<treeDataEdge> >& edgeTrees() const
{
return edgeTrees_;
}
const PtrList<indexedOctree<treeDataPoint> >& pointTrees() const
{
return pointTrees_;
}
// Query
//- Find nearest point on nearest feature edge
void findNearestEdge
(
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearFeature,
List<pointIndexHit>& nearInfo
) const;
//- Find nearest feature point. Is an index into feature points
// which itself is an index into the edgeMesh points.
// So the point index is
// pointTrees()[nearFeature].shapes().pointLabels()[nearIndex]
// Wip.
void findNearestPoint
(
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearFeature,
labelList& nearIndex
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //