ENH: snappyHexMesh: added proximity refinement

This commit is contained in:
mattijs
2013-05-27 11:15:39 +01:00
parent 35a944beea
commit 4b9f2ae6fc
7 changed files with 1147 additions and 27 deletions

View File

@ -46,9 +46,6 @@ defineTypeNameAndDebug(autoRefineDriver, 0);
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
@ -97,12 +94,14 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
(
refineParams.keepPoints()[0], // For now only use one.
refineParams.curvature(),
refineParams.planarAngle(),
true, // featureRefinement
false, // featureDistanceRefinement
false, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement
false, // gapRefinement
refineParams.maxGlobalCells(),
refineParams.maxLocalCells()
)
@ -208,12 +207,14 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
(
refineParams.keepPoints()[0],
refineParams.curvature(),
refineParams.planarAngle(),
false, // featureRefinement
false, // featureDistanceRefinement
false, // internalRefinement
true, // surfaceRefinement
true, // curvatureRefinement
false, // gapRefinement
refineParams.maxGlobalCells(),
refineParams.maxLocalCells()
)
@ -294,6 +295,371 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
}
Foam::label Foam::autoRefineDriver::gapOnlyRefine
(
const refinementParameters& refineParams,
const label maxIter
)
{
const fvMesh& mesh = meshRefiner_.mesh();
// Determine the maximum refinement level over all surfaces. This
// determines the minumum number of surface refinement iterations.
label maxIncrement = 0;
const labelList& maxLevel = meshRefiner_.surfaces().maxLevel();
const labelList& gapLevel = meshRefiner_.surfaces().gapLevel();
forAll(maxLevel, i)
{
maxIncrement = max(maxIncrement, gapLevel[i]-maxLevel[i]);
}
label iter = 0;
if (maxIncrement == 0)
{
return iter;
}
for (iter = 0; iter < maxIter; iter++)
{
Info<< nl
<< "Gap refinement iteration " << iter << nl
<< "--------------------------" << nl
<< endl;
// Determine cells to refine
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Only look at surface intersections (minLevel and surface curvature),
// do not do internal refinement (refinementShells)
labelList candidateCells
(
meshRefiner_.refineCandidates
(
refineParams.keepPoints()[0],
refineParams.curvature(),
refineParams.planarAngle(),
false, // featureRefinement
false, // featureDistanceRefinement
false, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement
true, // gapRefinement
refineParams.maxGlobalCells(),
refineParams.maxLocalCells()
)
);
if (debug&meshRefinement::MESH)
{
Pout<< "Dumping " << candidateCells.size()
<< " cells to cellSet candidateCellsFromGap." << endl;
cellSet c(mesh, "candidateCellsFromGap", candidateCells);
c.instance() = meshRefiner_.timeName();
c.write();
}
// Grow by one layer to make sure we're covering the gap
{
boolList isCandidateCell(mesh.nCells(), false);
forAll(candidateCells, i)
{
isCandidateCell[candidateCells[i]] = true;
}
for (label i=0; i<1; i++)
{
boolList newIsCandidateCell(isCandidateCell);
// Internal faces
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
label own = mesh.faceOwner()[faceI];
label nei = mesh.faceNeighbour()[faceI];
if (isCandidateCell[own] != isCandidateCell[nei])
{
newIsCandidateCell[own] = true;
newIsCandidateCell[nei] = true;
}
}
// Get coupled boundary condition values
boolList neiIsCandidateCell;
syncTools::swapBoundaryCellList
(
mesh,
isCandidateCell,
neiIsCandidateCell
);
// Boundary faces
for
(
label faceI = mesh.nInternalFaces();
faceI < mesh.nFaces();
faceI++
)
{
label own = mesh.faceOwner()[faceI];
label bFaceI = faceI-mesh.nInternalFaces();
if (isCandidateCell[own] != neiIsCandidateCell[bFaceI])
{
newIsCandidateCell[own] = true;
}
}
isCandidateCell.transfer(newIsCandidateCell);
}
label n = 0;
forAll(isCandidateCell, cellI)
{
if (isCandidateCell[cellI])
{
n++;
}
}
candidateCells.setSize(n);
n = 0;
forAll(isCandidateCell, cellI)
{
if (isCandidateCell[cellI])
{
candidateCells[n++] = cellI;
}
}
}
if (debug&meshRefinement::MESH)
{
Pout<< "Dumping " << candidateCells.size()
<< " cells to cellSet candidateCellsFromGapPlusBuffer." << endl;
cellSet c(mesh, "candidateCellsFromGapPlusBuffer", candidateCells);
c.instance() = meshRefiner_.timeName();
c.write();
}
labelList cellsToRefine
(
meshRefiner_.meshCutter().consistentRefinement
(
candidateCells,
true
)
);
Info<< "Determined cells to refine in = "
<< mesh.time().cpuTimeIncrement() << " s" << endl;
label nCellsToRefine = cellsToRefine.size();
reduce(nCellsToRefine, sumOp<label>());
Info<< "Selected for refinement : " << nCellsToRefine
<< " cells (out of " << mesh.globalData().nTotalCells()
<< ')' << endl;
// Stop when no cells to refine or have done minimum necessary
// iterations and not enough cells to refine.
if
(
nCellsToRefine == 0
|| (
iter >= maxIncrement
&& nCellsToRefine <= refineParams.minRefineCells()
)
)
{
Info<< "Stopping refining since too few cells selected."
<< nl << endl;
break;
}
if (debug)
{
const_cast<Time&>(mesh.time())++;
}
if
(
returnReduce
(
(mesh.nCells() >= refineParams.maxLocalCells()),
orOp<bool>()
)
)
{
meshRefiner_.balanceAndRefine
(
"gap refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
else
{
meshRefiner_.refineAndBalance
(
"gap refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
}
return iter;
}
Foam::label Foam::autoRefineDriver::danglingCellRefine
(
const refinementParameters& refineParams,
const label nFaces,
const label maxIter
)
{
const fvMesh& mesh = meshRefiner_.mesh();
label iter;
for (iter = 0; iter < maxIter; iter++)
{
Info<< nl
<< "Dangling coarse cells refinement iteration " << iter << nl
<< "--------------------------------------------" << nl
<< endl;
// Determine cells to refine
// ~~~~~~~~~~~~~~~~~~~~~~~~~
const cellList& cells = mesh.cells();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
labelList candidateCells;
{
cellSet candidateCellSet(mesh, "candidateCells", cells.size()/1000);
forAll(cells, cellI)
{
const cell& cFaces = cells[cellI];
label nIntFaces = 0;
forAll(cFaces, i)
{
label bFaceI = cFaces[i]-mesh.nInternalFaces();
if (bFaceI < 0)
{
nIntFaces++;
}
else
{
label patchI = pbm.patchID()[bFaceI];
if (pbm[patchI].coupled())
{
nIntFaces++;
}
}
}
if (nIntFaces == nFaces)
{
candidateCellSet.insert(cellI);
}
}
if (debug&meshRefinement::MESH)
{
Pout<< "Dumping " << candidateCellSet.size()
<< " cells to cellSet candidateCellSet." << endl;
candidateCellSet.instance() = meshRefiner_.timeName();
candidateCellSet.write();
}
candidateCells = candidateCellSet.toc();
}
labelList cellsToRefine
(
meshRefiner_.meshCutter().consistentRefinement
(
candidateCells,
true
)
);
Info<< "Determined cells to refine in = "
<< mesh.time().cpuTimeIncrement() << " s" << endl;
label nCellsToRefine = cellsToRefine.size();
reduce(nCellsToRefine, sumOp<label>());
Info<< "Selected for refinement : " << nCellsToRefine
<< " cells (out of " << mesh.globalData().nTotalCells()
<< ')' << endl;
// Stop when no cells to refine. No checking of minRefineCells since
// too few cells
if (nCellsToRefine == 0)
{
Info<< "Stopping refining since too few cells selected."
<< nl << endl;
break;
}
if (debug)
{
const_cast<Time&>(mesh.time())++;
}
if
(
returnReduce
(
(mesh.nCells() >= refineParams.maxLocalCells()),
orOp<bool>()
)
)
{
meshRefiner_.balanceAndRefine
(
"coarse cell refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
else
{
meshRefiner_.refineAndBalance
(
"coarse cell refinement iteration " + name(iter),
decomposer_,
distributor_,
cellsToRefine,
refineParams.maxLoadUnbalance()
);
}
}
return iter;
}
void Foam::autoRefineDriver::removeInsideCells
(
const refinementParameters& refineParams,
@ -371,12 +737,14 @@ Foam::label Foam::autoRefineDriver::shellRefine
(
refineParams.keepPoints()[0],
refineParams.curvature(),
refineParams.planarAngle(),
false, // featureRefinement
true, // featureDistanceRefinement
true, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement
false, // gapRefinement
refineParams.maxGlobalCells(),
refineParams.maxLocalCells()
)
@ -653,6 +1021,16 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
// Debug:test all is still synced across proc patches
meshRefiner_.checkData();
}
// Remove any now dangling parts
meshRefiner_.splitMeshRegions(refineParams.keepPoints()[0]);
if (debug)
{
// Debug:test all is still synced across proc patches
meshRefiner_.checkData();
}
Info<< "Merged free-standing baffles in = "
<< mesh.time().cpuTimeIncrement() << " s." << endl;
}
@ -732,6 +1110,12 @@ void Foam::autoRefineDriver::doRefine
100 // maxIter
);
gapOnlyRefine
(
refineParams,
100 // maxIter
);
// Remove cells (a certain distance) beyond surface intersections
removeInsideCells
(
@ -746,6 +1130,20 @@ void Foam::autoRefineDriver::doRefine
100 // maxIter
);
// Refine any hexes with 5 or 6 faces refined to make smooth edges
danglingCellRefine
(
refineParams,
21, // 1 coarse face + 5 refined faces
100 // maxIter
);
danglingCellRefine
(
refineParams,
24, // 0 coarse faces + 6 refined faces
100 // maxIter
);
// Introduce baffles at surface intersections. Remove sections unreachable
// from keepPoint.
baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -88,6 +88,21 @@ class autoRefineDriver
const label maxIter
);
//- Refine all cells in small gaps
label gapOnlyRefine
(
const refinementParameters& refineParams,
const label maxIter
);
//- Refine cells with almost all sides refined
label danglingCellRefine
(
const refinementParameters& refineParams,
const label nFaces,
const label maxIter
);
//- Remove all cells within intersected region
void removeInsideCells
(

View File

@ -335,6 +335,48 @@ private:
label& nRefine
) const;
//- Is local topology a small gap?
bool isGap
(
const scalar,
const vector&,
const vector&,
const vector&,
const vector&
) const;
//- Mark cell if local a gap topology or
bool checkProximity
(
const scalar planarCos,
const label nAllowRefine,
const label surfaceLevel,
const vector& surfaceLocation,
const vector& surfaceNormal,
const label cellI,
label& cellMaxLevel,
vector& cellMaxLocation,
vector& cellMaxNormal,
labelList& refineCell,
label& nRefine
) const;
//- Mark cells for surface proximity based refinement.
label markProximityRefinement
(
const scalar curvature,
const label nAllowRefine,
const labelList& neiLevel,
const pointField& neiCc,
labelList& refineCell,
label& nRefine
) const;
// Baffle handling
//- Get faces to repatch. Returns map from face to patch.
@ -676,12 +718,14 @@ public:
(
const point& keepPoint,
const scalar curvature,
const scalar planarAngle,
const bool featureRefinement,
const bool featureDistanceRefinement,
const bool internalRefinement,
const bool surfaceRefinement,
const bool curvatureRefinement,
const bool gapRefinement,
const label maxGlobalCells,
const label maxLocalCells
) const;

View File

@ -1965,7 +1965,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
const bool mergeFreeStanding,
const scalar freeStandingAngle,
const scalar planarAngle,
const dictionary& motionDict,
Time& runTime,
const labelList& globalToMasterPatch,
@ -2172,7 +2172,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
identity(mesh_.nFaces()-mesh_.nInternalFaces())
+mesh_.nInternalFaces()
),
freeStandingAngle
planarAngle
)
);

View File

@ -38,6 +38,7 @@ License
#include "featureEdgeMesh.H"
#include "Cloud.H"
//#include "globalIndex.H"
//#include "OBJstream.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -1022,6 +1023,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
end,
minLevel, // max level of surface has to be bigger
// than min level of neighbouring cells
surfaces_.maxLevel(),
surfaceNormal,
surfaceLevel
);
@ -1210,6 +1214,473 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
}
// Mark small gaps
bool Foam::meshRefinement::isGap
(
const scalar planarCos,
const vector& point0,
const vector& normal0,
const vector& point1,
const vector& normal1
) const
{
////- hits differ and angles are oppositeish
//return
// (mag(point0-point1) > mergeDistance())
// && ((normal0 & normal1) < (-1+planarCos));
//- hits differ and angles are oppositeish and
// hits have a normal distance
vector d = point1-point0;
scalar magD = mag(d);
if (magD > mergeDistance())
{
scalar cosAngle = (normal0 & normal1);
vector avg = vector::zero;
if (cosAngle < (-1+planarCos))
{
// Opposite normals
avg = 0.5*(normal0-normal1);
}
else if (cosAngle > (1-planarCos))
{
avg = 0.5*(normal0+normal1);
}
if (avg != vector::zero)
{
avg /= mag(avg);
d /= magD;
// Check average normal with respect to intersection locations
if (mag(avg&d) > Foam::cos(degToRad(45)))
{
return true;
}
else
{
return false;
}
}
else
{
return false;
}
}
else
{
return false;
}
}
bool Foam::meshRefinement::checkProximity
(
const scalar planarCos,
const label nAllowRefine,
const label surfaceLevel, // current intersection max level
const vector& surfaceLocation, // current intersection location
const vector& surfaceNormal, // current intersection normal
const label cellI,
label& cellMaxLevel, // cached max surface level for this cell
vector& cellMaxLocation, // cached surface normal for this cell
vector& cellMaxNormal, // cached surface normal for this cell
labelList& refineCell,
label& nRefine
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
// Test if surface applicable
if (surfaceLevel > cellLevel[cellI])
{
if (cellMaxLevel == -1)
{
// First visit of cell. Store
cellMaxLevel = surfaceLevel;
cellMaxLocation = surfaceLocation;
cellMaxNormal = surfaceNormal;
}
else
{
// Second or more visit.
// Check if
// - different location
// - opposite surface
bool closeSurfaces = isGap
(
planarCos,
cellMaxLocation,
cellMaxNormal,
surfaceLocation,
surfaceNormal
);
// Set normal to that of highest surface. Not really necessary
// over here but we reuse cellMax info when doing coupled faces.
if (surfaceLevel > cellMaxLevel)
{
cellMaxLevel = surfaceLevel;
cellMaxLocation = surfaceLocation;
cellMaxNormal = surfaceNormal;
}
if (closeSurfaces)
{
//Pout<< "Found gap:" << nl
// << " location:" << surfaceLocation
// << "\tnormal :" << surfaceNormal << nl
/// << " location:" << cellMaxLocation
// << "\tnormal :" << cellMaxNormal << nl
// << "\tcos :" << (surfaceNormal&cellMaxNormal) << nl
// << endl;
return markForRefine
(
surfaceLevel, // mark with any non-neg number.
nAllowRefine,
refineCell[cellI],
nRefine
);
}
}
}
// Did not reach refinement limit.
return true;
}
Foam::label Foam::meshRefinement::markProximityRefinement
(
const scalar planarCos,
const label nAllowRefine,
const labelList& neiLevel,
const pointField& neiCc,
labelList& refineCell,
label& nRefine
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
const pointField& cellCentres = mesh_.cellCentres();
label oldNRefine = nRefine;
// 1. local test: any cell on more than one surface gets refined
// (if its current level is < max of the surface max level)
// 2. 'global' test: any cell on only one surface with a neighbour
// on a different surface gets refined (if its current level etc.)
// Collect candidate faces (i.e. intersecting any surface and
// owner/neighbour not yet refined.
labelList testFaces(getRefineCandidateFaces(refineCell));
// Collect segments
pointField start(testFaces.size());
pointField end(testFaces.size());
labelList minLevel(testFaces.size());
forAll(testFaces, i)
{
label faceI = testFaces[i];
label own = mesh_.faceOwner()[faceI];
if (mesh_.isInternalFace(faceI))
{
label nei = mesh_.faceNeighbour()[faceI];
start[i] = cellCentres[own];
end[i] = cellCentres[nei];
minLevel[i] = min(cellLevel[own], cellLevel[nei]);
}
else
{
label bFaceI = faceI - mesh_.nInternalFaces();
start[i] = cellCentres[own];
end[i] = neiCc[bFaceI];
minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
}
}
// Extend segments a bit
{
const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
start -= smallVec;
end += smallVec;
}
// Test for all intersections (with surfaces of higher gap level than
// minLevel) and cache per cell the max surface level and the local normal
// on that surface.
labelList cellMaxLevel(mesh_.nCells(), -1);
vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
pointField cellMaxLocation(mesh_.nCells(), vector::zero);
{
// Per segment the normals of the surfaces
List<vectorList> surfaceLocation;
List<vectorList> surfaceNormal;
// Per segment the list of levels of the surfaces
labelListList surfaceLevel;
surfaces_.findAllHigherIntersections
(
start,
end,
minLevel, // gap level of surface has to be bigger
// than min level of neighbouring cells
surfaces_.gapLevel(),
surfaceLocation,
surfaceNormal,
surfaceLevel
);
// Clear out unnecessary data
start.clear();
end.clear();
minLevel.clear();
//// Extract per cell information on the surface with the highest max
//OBJstream str
//(
// mesh_.time().path()
// / "findAllHigherIntersections_"
// + mesh_.time().timeName()
// + ".obj"
//);
//// All intersections
//OBJstream str2
//(
// mesh_.time().path()
// / "findAllHigherIntersections2_"
// + mesh_.time().timeName()
// + ".obj"
//);
forAll(testFaces, i)
{
label faceI = testFaces[i];
label own = mesh_.faceOwner()[faceI];
const labelList& fLevels = surfaceLevel[i];
const vectorList& fPoints = surfaceLocation[i];
const vectorList& fNormals = surfaceNormal[i];
forAll(fLevels, hitI)
{
checkProximity
(
planarCos,
nAllowRefine,
fLevels[hitI],
fPoints[hitI],
fNormals[hitI],
own,
cellMaxLevel[own],
cellMaxLocation[own],
cellMaxNormal[own],
refineCell,
nRefine
);
}
if (mesh_.isInternalFace(faceI))
{
label nei = mesh_.faceNeighbour()[faceI];
forAll(fLevels, hitI)
{
checkProximity
(
planarCos,
nAllowRefine,
fLevels[hitI],
fPoints[hitI],
fNormals[hitI],
nei,
cellMaxLevel[nei],
cellMaxLocation[nei],
cellMaxNormal[nei],
refineCell,
nRefine
);
}
}
}
}
// 2. Find out a measure of surface curvature and region edges.
// Send over surface region and surface normal to neighbour cell.
labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
pointField neiBndMaxLocation(mesh_.nFaces()-mesh_.nInternalFaces());
vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
{
label bFaceI = faceI-mesh_.nInternalFaces();
label own = mesh_.faceOwner()[faceI];
neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
neiBndMaxLocation[bFaceI] = cellMaxLocation[own];
neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
}
syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLocation);
syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
// Loop over all faces. Could only be checkFaces.. except if they're coupled
// Internal faces
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label own = mesh_.faceOwner()[faceI];
label nei = mesh_.faceNeighbour()[faceI];
if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
{
// Have valid data on both sides. Check planarCos.
if
(
isGap
(
planarCos,
cellMaxLocation[own],
cellMaxNormal[own],
cellMaxLocation[nei],
cellMaxNormal[nei]
)
)
{
// See which side to refine
if (cellLevel[own] < cellMaxLevel[own])
{
if
(
!markForRefine
(
cellMaxLevel[own],
nAllowRefine,
refineCell[own],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching my cell"
<< " limit of " << mesh_.nCells()+7*nRefine
<< endl;
}
break;
}
}
if (cellLevel[nei] < cellMaxLevel[nei])
{
if
(
!markForRefine
(
cellMaxLevel[nei],
nAllowRefine,
refineCell[nei],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching my cell"
<< " limit of " << mesh_.nCells()+7*nRefine
<< endl;
}
break;
}
}
}
}
}
// Boundary faces
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
{
label own = mesh_.faceOwner()[faceI];
label bFaceI = faceI - mesh_.nInternalFaces();
if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
{
// Have valid data on both sides. Check planarCos.
if
(
isGap
(
planarCos,
cellMaxLocation[own],
cellMaxNormal[own],
neiBndMaxLocation[bFaceI],
neiBndMaxNormal[bFaceI]
)
)
{
if
(
!markForRefine
(
cellMaxLevel[own],
nAllowRefine,
refineCell[own],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching my cell"
<< " limit of " << mesh_.nCells()+7*nRefine
<< endl;
}
break;
}
}
}
}
if
(
returnReduce(nRefine, sumOp<label>())
> returnReduce(nAllowRefine, sumOp<label>())
)
{
Info<< "Reached refinement limit." << endl;
}
return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Calculate list of cells to refine. Gets for any edge (start - end)
@ -1221,12 +1692,14 @@ Foam::labelList Foam::meshRefinement::refineCandidates
(
const point& keepPoint,
const scalar curvature,
const scalar planarAngle,
const bool featureRefinement,
const bool featureDistanceRefinement,
const bool internalRefinement,
const bool surfaceRefinement,
const bool curvatureRefinement,
const bool gapRefinement,
const label maxGlobalCells,
const label maxLocalCells
) const
@ -1370,6 +1843,35 @@ Foam::labelList Foam::meshRefinement::refineCandidates
<< ": " << nCurv << " cells." << endl;
}
const scalar planarCos = Foam::cos(degToRad(planarAngle));
if
(
gapRefinement
&& (planarCos >= -1 && planarCos <= 1)
&& (max(surfaces_.gapLevel()) > -1)
)
{
Info<< "Specified gap level : " << max(surfaces_.gapLevel())
<< ", planar angle " << planarAngle << endl;
label nGap = markProximityRefinement
(
planarCos,
nAllowRefine,
neiLevel,
neiCc,
refineCell,
nRefine
);
Info<< "Marked for refinement due to close opposite surfaces "
<< ": " << nGap << " cells." << endl;
}
// Pack cells-to-refine
// ~~~~~~~~~~~~~~~~~~~~

View File

@ -116,10 +116,12 @@ Foam::refinementSurfaces::refinementSurfaces
labelList globalMinLevel(surfI, 0);
labelList globalMaxLevel(surfI, 0);
labelList globalLevelIncr(surfI, 0);
scalarField globalAngle(surfI, -GREAT);
PtrList<dictionary> globalPatchInfo(surfI);
List<Map<label> > regionMinLevel(surfI);
List<Map<label> > regionMaxLevel(surfI);
List<Map<label> > regionLevelIncr(surfI);
List<Map<scalar> > regionAngle(surfI);
List<Map<autoPtr<dictionary> > > regionPatchInfo(surfI);
@ -138,6 +140,31 @@ Foam::refinementSurfaces::refinementSurfaces
const labelPair refLevel(dict.lookup("level"));
globalMinLevel[surfI] = refLevel[0];
globalMaxLevel[surfI] = refLevel[1];
globalLevelIncr[surfI] = dict.lookupOrDefault
(
"gapLevelIncrement",
0
);
if
(
globalMinLevel[surfI] < 0
|| globalMaxLevel[surfI] < globalMinLevel[surfI]
|| globalLevelIncr[surfI] < 0
)
{
FatalErrorIn
(
"refinementSurfaces::refinementSurfaces"
"(const searchableSurfaces&, const dictionary>&"
) << "Illegal level specification for surface "
<< names_[surfI]
<< " : minLevel:" << globalMinLevel[surfI]
<< " maxLevel:" << globalMaxLevel[surfI]
<< " levelIncrement:" << globalLevelIncr[surfI]
<< exit(FatalError);
}
// Global zone names per surface
if (dict.readIfPresent("faceZone", faceZoneNames_[surfI]))
@ -244,6 +271,34 @@ Foam::refinementSurfaces::refinementSurfaces
regionMinLevel[surfI].insert(regionI, refLevel[0]);
regionMaxLevel[surfI].insert(regionI, refLevel[1]);
label levelIncr = regionDict.lookupOrDefault
(
"gapLevelIncrement",
0
);
regionLevelIncr[surfI].insert(regionI, levelIncr);
if
(
refLevel[0] < 0
|| refLevel[1] < refLevel[0]
|| levelIncr < 0
)
{
FatalErrorIn
(
"refinementSurfaces::refinementSurfaces"
"(const searchableSurfaces&, const dictionary>&"
) << "Illegal level specification for surface "
<< names_[surfI] << " region "
<< regionNames[regionI]
<< " : minLevel:" << refLevel[0]
<< " maxLevel:" << refLevel[1]
<< " levelIncrement:" << levelIncr
<< exit(FatalError);
}
if (regionDict.found("perpendicularAngle"))
{
@ -286,6 +341,8 @@ Foam::refinementSurfaces::refinementSurfaces
minLevel_ = 0;
maxLevel_.setSize(nRegions);
maxLevel_ = 0;
gapLevel_.setSize(nRegions);
gapLevel_ = -1;
perpendicularAngle_.setSize(nRegions);
perpendicularAngle_ = -GREAT;
patchInfo_.setSize(nRegions);
@ -301,6 +358,10 @@ Foam::refinementSurfaces::refinementSurfaces
label globalRegionI = regionOffset_[surfI] + i;
minLevel_[globalRegionI] = globalMinLevel[surfI];
maxLevel_[globalRegionI] = globalMaxLevel[surfI];
gapLevel_[globalRegionI] =
maxLevel_[globalRegionI]
+ globalLevelIncr[surfI];
perpendicularAngle_[globalRegionI] = globalAngle[surfI];
if (globalPatchInfo.set(surfI))
{
@ -319,24 +380,9 @@ Foam::refinementSurfaces::refinementSurfaces
minLevel_[globalRegionI] = iter();
maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
// Check validity
if
(
minLevel_[globalRegionI] < 0
|| maxLevel_[globalRegionI] < minLevel_[globalRegionI]
)
{
FatalErrorIn
(
"refinementSurfaces::refinementSurfaces"
"(const searchableSurfaces&, const dictionary>&"
) << "Illegal level or layer specification for surface "
<< names_[surfI]
<< " : minLevel:" << minLevel_[globalRegionI]
<< " maxLevel:" << maxLevel_[globalRegionI]
<< exit(FatalError);
}
gapLevel_[globalRegionI] =
maxLevel_[globalRegionI]
+ regionLevelIncr[surfI][iter.key()];
}
forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
{
@ -714,6 +760,8 @@ void Foam::refinementSurfaces::findAllHigherIntersections
const pointField& end,
const labelList& currentLevel, // current cell refinement level
const labelList& globalRegionLevel,
List<vectorList>& surfaceNormal,
labelListList& surfaceLevel
) const
@ -777,7 +825,7 @@ void Foam::refinementSurfaces::findAllHigherIntersections
label region = globalRegion(surfI, surfRegion[i]);
label pointI = pointMap[i];
if (maxLevel_[region] > currentLevel[pointI])
if (globalRegionLevel[region] > currentLevel[pointI])
{
// Append to pointI info
label sz = surfaceNormal[pointI].size();
@ -785,7 +833,95 @@ void Foam::refinementSurfaces::findAllHigherIntersections
surfaceNormal[pointI][sz] = surfNormal[i];
surfaceLevel[pointI].setSize(sz+1);
surfaceLevel[pointI][sz] = maxLevel_[region];
surfaceLevel[pointI][sz] = globalRegionLevel[region];
}
}
}
}
void Foam::refinementSurfaces::findAllHigherIntersections
(
const pointField& start,
const pointField& end,
const labelList& currentLevel, // current cell refinement level
const labelList& globalRegionLevel,
List<pointList>& surfaceLocation,
List<vectorList>& surfaceNormal,
labelListList& surfaceLevel
) const
{
surfaceLevel.setSize(start.size());
surfaceNormal.setSize(start.size());
surfaceLocation.setSize(start.size());
if (surfaces_.empty())
{
return;
}
// Work arrays
List<List<pointIndexHit> > hitInfo;
labelList pRegions;
vectorField pNormals;
forAll(surfaces_, surfI)
{
allGeometry_[surfaces_[surfI]].findLineAll(start, end, hitInfo);
// Repack hits for surface into flat list
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// To avoid overhead of calling getRegion for every point
label n = 0;
forAll(hitInfo, pointI)
{
n += hitInfo[pointI].size();
}
List<pointIndexHit> surfInfo(n);
labelList pointMap(n);
n = 0;
forAll(hitInfo, pointI)
{
const List<pointIndexHit>& pHits = hitInfo[pointI];
forAll(pHits, i)
{
surfInfo[n] = pHits[i];
pointMap[n] = pointI;
n++;
}
}
labelList surfRegion(n);
vectorField surfNormal(n);
allGeometry_[surfaces_[surfI]].getRegion(surfInfo, surfRegion);
allGeometry_[surfaces_[surfI]].getNormal(surfInfo, surfNormal);
// Extract back into pointwise
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(surfRegion, i)
{
label region = globalRegion(surfI, surfRegion[i]);
label pointI = pointMap[i];
if (globalRegionLevel[region] > currentLevel[pointI])
{
// Append to pointI info
label sz = surfaceNormal[pointI].size();
surfaceLocation[pointI].setSize(sz+1);
surfaceLocation[pointI][sz] = surfInfo[i].hitPoint();
surfaceNormal[pointI].setSize(sz+1);
surfaceNormal[pointI][sz] = surfNormal[i];
surfaceLevel[pointI].setSize(sz+1);
surfaceLevel[pointI][sz] = globalRegionLevel[region];
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -51,6 +51,8 @@ class searchableSurfaces;
class shellSurfaces;
class triSurfaceMesh;
typedef List<point> pointList;
/*---------------------------------------------------------------------------*\
Class refinementSurfaces Declaration
\*---------------------------------------------------------------------------*/
@ -120,6 +122,9 @@ private:
//- From global region number to refinement level
labelList maxLevel_;
//- From global region number to small-gap level
labelList gapLevel_;
//- From global region number to perpendicular angle
scalarField perpendicularAngle_;
@ -226,6 +231,12 @@ public:
return maxLevel_;
}
//- From global region number to small gap refinement level
const labelList& gapLevel() const
{
return gapLevel_;
}
//- From global region number to perpendicular angle
const scalarField& perpendicularAngle() const
{
@ -295,11 +306,25 @@ public:
const pointField& start,
const pointField& end,
const labelList& currentLevel, // current cell refinement level
const labelList& globalRegionLevel, // level per surfregion
List<vectorList>& surfaceNormal,
labelListList& surfaceLevel
) const;
//- Find all intersections of edge. Unsorted order.
void findAllHigherIntersections
(
const pointField& start,
const pointField& end,
const labelList& currentLevel, // current cell refinement level
const labelList& globalRegionLevel, // level per surfregion
List<pointList>& surfaceLocation,
List<vectorList>& surfaceNormal,
labelListList& surfaceLevel
) const;
//- Find intersection nearest to the endpoints. surface1,2 are
// not indices into surfacesToTest but refinement surface indices.
// Returns surface, region on surface (so not global surface)