diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index 67a9c3a9444a3803f8c43dd270b2b7330219b32a..7b4d8577c4e7ace8293bf50a703c234bcc3193f2 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -1075,6 +1075,26 @@ int main(int argc, char *argv[])
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+ // Optionally read limit shells
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ const dictionary limitDict(refineDict.subOrEmptyDict("limitRegions"));
+
+ if (!limitDict.empty())
+ {
+ Info<< "Reading limit shells." << endl;
+ }
+
+ shellSurfaces limitShells(allGeometry, limitDict);
+
+ if (!limitDict.empty())
+ {
+ Info<< "Read refinement shells in = "
+ << mesh.time().cpuTimeIncrement() << " s" << nl << endl;
+ }
+
+
+
// Read feature meshes
// ~~~~~~~~~~~~~~~~~~~
@@ -1105,7 +1125,8 @@ int main(int argc, char *argv[])
overwrite, // overwrite mesh files?
surfaces, // for surface intersection refinement
features, // for feature edges/point based refinement
- shells // for volume (inside/outside) refinement
+ shells, // for volume (inside/outside) refinement
+ limitShells // limit of volume refinement
);
Info<< "Calculated surface intersections in = "
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 9b759b5ca4f6628dd9e11117d62697a0b12e0aef..ea073cd9f070647735f623262a850d4be8f4d220 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -95,10 +95,10 @@ castellatedMeshControls
// actually be a lot less.
maxGlobalCells 2000000;
- // The surface refinement loop might spend lots of iterations refining just a
- // few cells. This setting will cause refinement to stop if <= minimumRefine
- // are selected for refinement. Note: it will at least do one iteration
- // (unless the number of cells to refine is 0)
+ // The surface refinement loop might spend lots of iterations refining just
+ // a few cells. This setting will cause refinement to stop if
+ // <= minimumRefine cells are selected for refinement. Note: it will
+ // at least do one iteration (unless the number of cells to refine is 0)
minRefinementCells 0;
// Allow a certain level of imbalance during refining
@@ -249,6 +249,13 @@ castellatedMeshControls
//}
}
+
+ // Limit refinement in geometric region
+ limitRegions
+ {
+ }
+
+
// Mesh selection
// ~~~~~~~~~~~~~~
diff --git a/src/OpenFOAM/algorithms/indexedOctree/volumeType.H b/src/OpenFOAM/algorithms/indexedOctree/volumeType.H
index 7bdfc043e6024e5b3a2417dbd4d8960b4927515a..3d845c414bef3fda5df79c57cb64d471f2dd3218 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/volumeType.H
+++ b/src/OpenFOAM/algorithms/indexedOctree/volumeType.H
@@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -95,6 +95,12 @@ public:
t_(t)
{}
+ //- Construct from integer
+ explicit volumeType(const int t)
+ :
+ t_(static_cast<volumeType::type>(t))
+ {}
+
// Member Functions
diff --git a/src/mesh/autoMesh/Make/files b/src/mesh/autoMesh/Make/files
index 081d70afbd52185b77011cfa13f68f8f07fea44e..29001bafb78ee672699f01948483be8bf7772d89 100644
--- a/src/mesh/autoMesh/Make/files
+++ b/src/mesh/autoMesh/Make/files
@@ -15,6 +15,7 @@ $(autoHexMesh)/meshRefinement/meshRefinement.C
$(autoHexMesh)/meshRefinement/meshRefinementMerge.C
$(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
$(autoHexMesh)/meshRefinement/meshRefinementRefine.C
+$(autoHexMesh)/meshRefinement/meshRefinementGapRefine.C
$(autoHexMesh)/meshRefinement/patchFaceOrientation.C
$(autoHexMesh)/refinementFeatures/refinementFeatures.C
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index be908b76bbc9588fb48a586579613ac1539a62bd..e7e2217f27d3e9bb737c6bee289fef31e6f848b8 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -105,7 +105,9 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
false, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement
+ false, // smallFeatureRefinement
false, // gapRefinement
+ false, // bigGapRefinement
refineParams.maxGlobalCells(),
refineParams.maxLocalCells()
)
@@ -179,6 +181,127 @@ Foam::label Foam::autoRefineDriver::featureEdgeRefine
}
+Foam::label Foam::autoRefineDriver::smallFeatureRefine
+(
+ const refinementParameters& refineParams,
+ const label maxIter
+)
+{
+ const fvMesh& mesh = meshRefiner_.mesh();
+
+
+ label iter = 0;
+
+ // See if any surface has an extendedGapLevel
+ labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
+ labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
+
+ if (max(surfaceMaxLevel) == 0 && max(shellMaxLevel) == 0)
+ {
+ return iter;
+ }
+
+ for (; iter < maxIter; iter++)
+ {
+ Info<< nl
+ << "Small surface feature refinement iteration " << iter << nl
+ << "--------------------------------------------" << nl
+ << endl;
+
+
+ // Determine cells to refine
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ labelList candidateCells
+ (
+ meshRefiner_.refineCandidates
+ (
+ refineParams.locationsInMesh(),
+ refineParams.curvature(),
+ refineParams.planarAngle(),
+
+ false, // featureRefinement
+ false, // featureDistanceRefinement
+ false, // internalRefinement
+ false, // surfaceRefinement
+ false, // curvatureRefinement
+ true, // smallFeatureRefinement
+ false, // gapRefinement
+ false, // bigGapRefinement
+ refineParams.maxGlobalCells(),
+ refineParams.maxLocalCells()
+ )
+ );
+
+ 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)
+ {
+ 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
+ (
+ "small feature refinement iteration " + name(iter),
+ decomposer_,
+ distributor_,
+ cellsToRefine,
+ refineParams.maxLoadUnbalance()
+ );
+ }
+ else
+ {
+ meshRefiner_.refineAndBalance
+ (
+ "small feature refinement iteration " + name(iter),
+ decomposer_,
+ distributor_,
+ cellsToRefine,
+ refineParams.maxLoadUnbalance()
+ );
+ }
+ }
+ return iter;
+}
+
+
Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
(
const refinementParameters& refineParams,
@@ -218,7 +341,9 @@ Foam::label Foam::autoRefineDriver::surfaceOnlyRefine
false, // internalRefinement
true, // surfaceRefinement
true, // curvatureRefinement
+ false, // smallFeatureRefinement
false, // gapRefinement
+ false, // bigGapRefinement
refineParams.maxGlobalCells(),
refineParams.maxLocalCells()
)
@@ -352,7 +477,9 @@ Foam::label Foam::autoRefineDriver::gapOnlyRefine
false, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement
+ false, // smallFeatureRefinement
true, // gapRefinement
+ false, // bigGapRefinement
refineParams.maxGlobalCells(),
refineParams.maxLocalCells()
)
@@ -527,6 +654,146 @@ Foam::label Foam::autoRefineDriver::gapOnlyRefine
}
+Foam::label Foam::autoRefineDriver::bigGapOnlyRefine
+(
+ const refinementParameters& refineParams,
+ const label maxIter
+)
+{
+ const fvMesh& mesh = meshRefiner_.mesh();
+
+ label iter = 0;
+
+ // See if any surface has an extendedGapLevel
+ labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
+ labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
+
+ label overallMaxLevel(max(max(surfaceMaxLevel), max(shellMaxLevel)));
+
+ if (overallMaxLevel == 0)
+ {
+ return iter;
+ }
+
+
+ for (; iter < maxIter; iter++)
+ {
+ Info<< nl
+ << "Big gap refinement iteration " << iter << nl
+ << "------------------------------" << nl
+ << endl;
+
+
+ // Determine cells to refine
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ labelList candidateCells
+ (
+ meshRefiner_.refineCandidates
+ (
+ refineParams.locationsInMesh(),
+ refineParams.curvature(),
+ refineParams.planarAngle(),
+
+ false, // featureRefinement
+ false, // featureDistanceRefinement
+ false, // internalRefinement
+ false, // surfaceRefinement
+ false, // curvatureRefinement
+ false, // smallFeatureRefinement
+ false, // gapRefinement
+ true, // bigGapRefinement
+ refineParams.maxGlobalCells(),
+ refineParams.maxLocalCells()
+ )
+ );
+
+
+ if (debug&meshRefinement::MESH)
+ {
+ Pout<< "Dumping " << candidateCells.size()
+ << " cells to cellSet candidateCellsFromBigGap." << endl;
+ cellSet c(mesh, "candidateCellsFromBigGap", 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 >= overallMaxLevel
+ && 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
+ (
+ "big gap refinement iteration " + name(iter),
+ decomposer_,
+ distributor_,
+ cellsToRefine,
+ refineParams.maxLoadUnbalance()
+ );
+ }
+ else
+ {
+ meshRefiner_.refineAndBalance
+ (
+ "big gap refinement iteration " + name(iter),
+ decomposer_,
+ distributor_,
+ cellsToRefine,
+ refineParams.maxLoadUnbalance()
+ );
+ }
+ }
+ return iter;
+}
+
+
Foam::label Foam::autoRefineDriver::danglingCellRefine
(
const refinementParameters& refineParams,
@@ -1108,7 +1375,9 @@ Foam::label Foam::autoRefineDriver::shellRefine
true, // internalRefinement
false, // surfaceRefinement
false, // curvatureRefinement
+ false, // smallFeatureRefinement
false, // gapRefinement
+ false, // bigGapRefinement
refineParams.maxGlobalCells(),
refineParams.maxLocalCells()
)
@@ -1630,6 +1899,13 @@ void Foam::autoRefineDriver::doRefine
0 // min cells to refine
);
+ // Refine cells that contain a gap
+ smallFeatureRefine
+ (
+ refineParams,
+ 100 // maxIter
+ );
+
// Refine based on surface
surfaceOnlyRefine
(
@@ -1650,6 +1926,13 @@ void Foam::autoRefineDriver::doRefine
1 // nBufferLayers
);
+ // Refine consistently across narrow gaps (a form of shell refinement)
+ bigGapOnlyRefine
+ (
+ refineParams,
+ 100 // maxIter
+ );
+
// Internal mesh refinement
shellRefine
(
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
index 11ed1d9454d0b5c7d6db927e0fa25c0d9fdb1ef4..e952af9fd60e00aa5aaaa26b08ce694b6ce2c4af 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H
@@ -84,6 +84,13 @@ class autoRefineDriver
const label minRefine
);
+ //- Refine all cells containing small surface features
+ label smallFeatureRefine
+ (
+ const refinementParameters& refineParams,
+ const label maxIter
+ );
+
//- Refine all cells interacting with the surface
label surfaceOnlyRefine
(
@@ -98,6 +105,13 @@ class autoRefineDriver
const label maxIter
);
+ //- Refine all cells in large gaps
+ label bigGapOnlyRefine
+ (
+ const refinementParameters& refineParams,
+ const label maxIter
+ );
+
//- Refine cells with almost all sides refined
label danglingCellRefine
(
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
index e4bfa8840c994ca689be4dcfb46ece5ae3e445fd..93818eb8d66013d54885341c5c18712431194694 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
@@ -215,6 +215,55 @@ void Foam::meshRefinement::calcNeighbourData
}
+void Foam::meshRefinement::calcCellCellRays
+(
+ const pointField& neiCc,
+ const labelList& neiLevel,
+ const labelList& testFaces,
+ pointField& start,
+ pointField& end,
+ labelList& minLevel
+) const
+{
+ const labelList& cellLevel = meshCutter_.cellLevel();
+ const pointField& cellCentres = mesh_.cellCentres();
+
+ start.setSize(testFaces.size());
+ end.setSize(testFaces.size());
+ minLevel.setSize(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(ROOTSMALL*(end-start));
+ start -= smallVec;
+ end += smallVec;
+ }
+}
+
+
// Find intersections of edges (given by their two endpoints) with surfaces.
// Returns first intersection if there are more than one.
void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
@@ -1187,7 +1236,8 @@ Foam::meshRefinement::meshRefinement
const bool overwrite,
const refinementSurfaces& surfaces,
const refinementFeatures& features,
- const shellSurfaces& shells
+ const shellSurfaces& shells,
+ const shellSurfaces& limitShells
)
:
mesh_(mesh),
@@ -1197,6 +1247,7 @@ Foam::meshRefinement::meshRefinement
surfaces_(surfaces),
features_(features),
shells_(shells),
+ limitShells_(limitShells),
meshCutter_
(
mesh,
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 31ef9f2e1a286d2bbc2cc9fa7c35a15bef20ebef..b2cbbd983f22d4d8de36ffdd8edc7ea6314bd278 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -34,6 +34,7 @@ Description
SourceFiles
meshRefinement.C
meshRefinementBaffles.C
+ meshRefinementGapRefine.C
meshRefinementMerge.C
meshRefinementProblemCells.C
meshRefinementRefine.C
@@ -53,6 +54,7 @@ SourceFiles
#include "pointIndexHit.H"
#include "wordPairHashTable.H"
#include "surfaceZonesInfo.H"
+#include "volumeType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -176,6 +178,9 @@ private:
//- All shell-refinement interaction
const shellSurfaces& shells_;
+ //- All limit-refinement interaction
+ const shellSurfaces& limitShells_;
+
//- Refinement engine
hexRef8 meshCutter_;
@@ -226,10 +231,18 @@ private:
);
//- Calculate coupled boundary end vector and refinement level
- void calcNeighbourData
+ void calcNeighbourData(labelList& neiLevel, pointField& neiCc) const;
+
+ //- Calculate rays from cell-centre to cell-centre and corresponding
+ // min cell refinement level
+ void calcCellCellRays
(
- labelList& neiLevel,
- pointField& neiCc
+ const pointField& neiCc,
+ const labelList& neiLevel,
+ const labelList& testFaces,
+ pointField& start,
+ pointField& end,
+ labelList& minLevel
) const;
//- Remove cells. Put exposedFaces into exposedPatchIDs.
@@ -241,12 +254,12 @@ private:
removeCells& cellRemover
);
- // Get cells which are inside any closed surface. Note that
- // all closed surfaces
- // will have already been oriented to have keepPoint outside.
+ //- Get cells which are inside any closed surface. Note that
+ // all closed surfaces
+ // will have already been oriented to have keepPoint outside.
labelList getInsideCells(const word&) const;
- // Do all to remove inside cells
+ //- Do all to remove inside cells
autoPtr<mapPolyMesh> removeInsideCells
(
const string& msg,
@@ -302,6 +315,14 @@ private:
label& nRefine
) const;
+ //- Unmark cells for refinement based on limit-shells. Return number
+ // of limited cells.
+ label unmarkInternalRefinement
+ (
+ labelList& refineCell,
+ label& nRefine
+ ) const;
+
//- Collect faces that are intersected and whose neighbours aren't
// yet marked for refinement.
labelList getRefineCandidateFaces
@@ -319,6 +340,83 @@ private:
label& nRefine
) const;
+ //- Mark cells intersected by the surface if they are inside
+ // close gaps
+ label markSurfaceGapRefinement
+ (
+ const scalar planarCos,
+ const label nAllowRefine,
+ const labelList& neiLevel,
+ const pointField& neiCc,
+
+ labelList& refineCell,
+ label& nRefine
+ ) const;
+
+ //- Generate single ray for cell to detect gap
+ bool generateRay
+ (
+ const bool useSurfaceNormal,
+
+ const point& nearPoint,
+ const vector& nearNormal,
+ const FixedList<label, 3>& gapInfo,
+ const volumeType& mode,
+
+ const point& cc,
+ const label cLevel,
+
+ point& start,
+ point& end,
+ point& start2,
+ point& end2
+ ) const;
+
+ //- Select candidate cells (cells inside a shell with gapLevel
+ // specified)
+ void selectGapCandidates
+ (
+ const labelList& refineCell,
+ const label nRefine,
+
+ labelList& cellMap,
+ List<FixedList<label, 3> >& shellGapInfo,
+ List<volumeType>& shellGapMode
+ ) const;
+
+ //- Merge gap information coming from shell and from surface
+ // (surface wins)
+ void mergeGapInfo
+ (
+ const FixedList<label, 3>& shellGapInfo,
+ const volumeType shellGapMode,
+ const FixedList<label, 3>& surfGapInfo,
+ const volumeType surfGapMode,
+ FixedList<label, 3>& gapInfo,
+ volumeType& gapMode
+ ) const;
+
+ //- Mark cells for non-surface intersection based gap refinement
+ label markInternalGapRefinement
+ (
+ const scalar planarCos,
+ const label nAllowRefine,
+ labelList& refineCell,
+ label& nRefine
+ ) const;
+
+ //- Refine cells containing small gaps
+ label markSmallFeatureRefinement
+ (
+ const scalar planarCos,
+ const label nAllowRefine,
+ const labelList& neiLevel,
+ const pointField& neiCc,
+
+ labelList& refineCell,
+ label& nRefine
+ ) const;
+
//- Helper: count number of normals1 that are in normals2
label countMatches
(
@@ -387,6 +485,7 @@ private:
void getIntersections
(
const labelList& surfacesToTest,
+ const labelList& neiLevel,
const pointField& neiCc,
const labelList& testFaces,
@@ -567,7 +666,7 @@ private:
const labelList& cellToZone,
const labelList& neiCellZone,
const labelList& faceToZone,
- const boolList& meshFlipMap,
+ const PackedBoolList& meshFlipMap,
polyTopoChange& meshMod
) const;
@@ -631,7 +730,7 @@ private:
const labelList& nMasterFaces,
const labelList& faceToZone,
const Map<label>& zoneToOrientation,
- boolList& meshFlipMap
+ PackedBoolList& meshFlipMap
) const;
@@ -657,6 +756,7 @@ public:
const bool overwrite,
const refinementSurfaces&,
const refinementFeatures&,
+ const shellSurfaces&,
const shellSurfaces&
);
@@ -855,7 +955,9 @@ public:
const bool internalRefinement,
const bool surfaceRefinement,
const bool curvatureRefinement,
+ const bool smallFeatureRefinement,
const bool gapRefinement,
+ const bool bigGapRefinement,
const label maxGlobalCells,
const label maxLocalCells
) const;
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 6f500c117a10304be1cf1330ea45baec56674473..302882a3fb1d0e745125f61eb7404c74462b9bbd 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -175,6 +175,7 @@ Foam::label Foam::meshRefinement::createBaffle
void Foam::meshRefinement::getIntersections
(
const labelList& surfacesToTest,
+ const labelList& neiLevel,
const pointField& neiCc,
const labelList& testFaces,
@@ -198,8 +199,6 @@ void Foam::meshRefinement::getIntersections
<< str().name() << nl << endl;
}
- const pointField& cellCentres = mesh_.cellCentres();
-
globalRegion1.setSize(mesh_.nFaces());
globalRegion1 = -1;
@@ -213,29 +212,17 @@ void Foam::meshRefinement::getIntersections
pointField start(testFaces.size());
pointField end(testFaces.size());
- forAll(testFaces, i)
{
- label faceI = testFaces[i];
-
- label own = mesh_.faceOwner()[faceI];
-
- if (mesh_.isInternalFace(faceI))
- {
- start[i] = cellCentres[own];
- end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
- }
- else
- {
- start[i] = cellCentres[own];
- end[i] = neiCc[faceI-mesh_.nInternalFaces()];
- }
- }
-
- // Extend segments a bit
- {
- const vectorField smallVec(ROOTSMALL*(end-start));
- start -= smallVec;
- end += smallVec;
+ labelList minLevel;
+ calcCellCellRays
+ (
+ neiCc,
+ neiLevel,
+ testFaces,
+ start,
+ end,
+ minLevel
+ );
}
@@ -319,6 +306,7 @@ void Foam::meshRefinement::getBafflePatches
getIntersections
(
surfaceZonesInfo::getUnnamedSurfaces(surfaces_.surfZones()),
+ neiLevel,
neiCc,
testFaces,
globalRegion1,
@@ -2608,7 +2596,7 @@ void Foam::meshRefinement::consistentOrientation
const labelList& nMasterFacesPerEdge,
const labelList& faceToZone,
const Map<label>& zoneToOrientation,
- boolList& meshFlipMap
+ PackedBoolList& meshFlipMap
) const
{
const polyBoundaryMesh& bm = mesh_.boundaryMesh();
@@ -2858,7 +2846,7 @@ void Foam::meshRefinement::zonify
const labelList& cellToZone,
const labelList& neiCellZone,
const labelList& faceToZone,
- const boolList& meshFlipMap,
+ const PackedBoolList& meshFlipMap,
polyTopoChange& meshMod
) const
{
@@ -3902,6 +3890,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
getIntersections
(
identity(surfaces_.surfaces().size()), // surfacesToTest,
+ neiLevel,
neiCc,
intersectedFaces(), // testFaces
globalRegion1,
@@ -4297,7 +4286,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
}
- // 2. Combine faceZoneNames allocated on different processors
+ // 2.Combine faceZoneNames allocated on different processors
Pstream::mapCombineGather(zonesToFaceZone, eqOp<word>());
Pstream::mapCombineScatter(zonesToFaceZone);
@@ -4438,7 +4427,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// - do a consistent orientation
// - check number of faces with consistent orientation
// - if <0 flip the whole patch
- boolList meshFlipMap(mesh_.nFaces(), false);
+ PackedBoolList meshFlipMap(mesh_.nFaces(), false);
{
// Collect all data on zone faces without cellZones on either side.
const indirectPrimitivePatch patch
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C
new file mode 100644
index 0000000000000000000000000000000000000000..62a285928c128659b52964cd70621f37b96d3e2c
--- /dev/null
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementGapRefine.C
@@ -0,0 +1,1546 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "meshRefinement.H"
+#include "Time.H"
+#include "refinementSurfaces.H"
+#include "refinementFeatures.H"
+#include "shellSurfaces.H"
+#include "OBJstream.H"
+#include "triSurfaceMesh.H"
+#include "treeDataCell.H"
+#include "searchableSurfaces.H"
+#include "DynamicField.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::label Foam::meshRefinement::markSurfaceGapRefinement
+(
+ 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();
+
+ // Get the gap level for the shells
+ const labelList maxLevel(shells_.maxGapLevel());
+
+ label oldNRefine = nRefine;
+
+ if (max(maxLevel) > 0)
+ {
+ // Use cached surfaceIndex_ to detect if any intersection. If so
+ // re-intersect to determine level wanted.
+
+ // Collect candidate faces
+ // ~~~~~~~~~~~~~~~~~~~~~~~
+
+ labelList testFaces(getRefineCandidateFaces(refineCell));
+
+ // Collect segments
+ // ~~~~~~~~~~~~~~~~
+
+ pointField start(testFaces.size());
+ pointField end(testFaces.size());
+ labelList minLevel(testFaces.size());
+
+ calcCellCellRays
+ (
+ neiCc,
+ neiLevel,
+ testFaces,
+ start,
+ end,
+ minLevel
+ );
+
+
+ // Collect cells to test for inside/outside in shell
+ labelList cellToCompact(mesh_.nCells(), -1);
+ labelList bFaceToCompact(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
+ List<FixedList<label, 3> > shellGapInfo;
+ List<volumeType> shellGapMode;
+ {
+ DynamicField<point> compactToCc(mesh_.nCells()/10);
+ DynamicList<label> compactToLevel(compactToCc.capacity());
+ forAll(testFaces, i)
+ {
+ label faceI = testFaces[i];
+ label own = mesh_.faceOwner()[faceI];
+ if (cellToCompact[own] == -1)
+ {
+ cellToCompact[own] = compactToCc.size();
+ compactToCc.append(cellCentres[own]);
+ compactToLevel.append(cellLevel[own]);
+ }
+ if (mesh_.isInternalFace(faceI))
+ {
+ label nei = mesh_.faceNeighbour()[faceI];
+ if (cellToCompact[nei] == -1)
+ {
+ cellToCompact[nei] = compactToCc.size();
+ compactToCc.append(cellCentres[nei]);
+ compactToLevel.append(cellLevel[nei]);
+ }
+ }
+ else
+ {
+ label bFaceI = faceI - mesh_.nInternalFaces();
+ if (bFaceToCompact[bFaceI] == -1)
+ {
+ bFaceToCompact[bFaceI] = compactToCc.size();
+ compactToCc.append(neiCc[bFaceI]);
+ compactToLevel.append(neiLevel[bFaceI]);
+ }
+ }
+ }
+
+ shells_.findHigherGapLevel
+ (
+ compactToCc,
+ compactToLevel,
+ shellGapInfo,
+ shellGapMode
+ );
+ }
+
+
+ //OBJstream str(mesh_.time().timePath()/"markSurfaceRefinement.obj");
+ //Info<< "Dumping rays to " << str.name( ) << endl;
+
+ const List<FixedList<label, 3> >& extendedGapLevel =
+ surfaces_.extendedGapLevel();
+ const List<volumeType>& extendedGapMode =
+ surfaces_.extendedGapMode();
+
+ labelList surface1;
+ List<pointIndexHit> hit1;
+ labelList region1;
+ vectorField normal1;
+
+ {
+ labelList surface2;
+ List<pointIndexHit> hit2;
+ labelList region2;
+ vectorField normal2;
+
+ surfaces_.findNearestIntersection
+ (
+ identity(surfaces_.surfaces().size()),
+ start,
+ end,
+
+ surface1,
+ hit1,
+ region1,
+ normal1,
+
+ surface2,
+ hit2,
+ region2,
+ normal2
+ );
+ }
+
+
+ start.setSize(2*surface1.size());
+ end.setSize(2*surface1.size());
+ labelList map(2*surface1.size());
+ pointField start2(2*surface1.size());
+ pointField end2(2*surface1.size());
+ labelList cellMap(2*surface1.size());
+
+ label compactI = 0;
+
+ forAll(surface1, i)
+ {
+ label surfI = surface1[i];
+
+ if (surfI != -1)
+ {
+ label globalRegionI = surfaces_.globalRegion(surfI, region1[i]);
+
+ label faceI = testFaces[i];
+ const point& surfPt = hit1[i].hitPoint();
+
+ label own = mesh_.faceOwner()[faceI];
+ if (cellToCompact[own] != -1)
+ {
+ // Combine info from shell and surface
+ label compactI = cellToCompact[own];
+ FixedList<label, 3> gapInfo;
+ volumeType gapMode;
+ mergeGapInfo
+ (
+ shellGapInfo[compactI],
+ shellGapMode[compactI],
+ extendedGapLevel[globalRegionI],
+ extendedGapMode[globalRegionI],
+ gapInfo,
+ gapMode
+ );
+
+ const point& cc = cellCentres[own];
+ bool okRay = generateRay
+ (
+ false,
+ surfPt,
+ normal1[i],
+ gapInfo,
+ gapMode,
+ surfPt+((cc-surfPt)&normal1[i])*normal1[i],
+ cellLevel[own],
+
+ start[compactI],
+ end[compactI],
+ start2[compactI],
+ end2[compactI]
+ );
+ cellMap[compactI] = own;
+
+ if (okRay)
+ {
+ map[compactI++] = i;
+ }
+ }
+ if (mesh_.isInternalFace(faceI))
+ {
+ label nei = mesh_.faceNeighbour()[faceI];
+ if (cellToCompact[nei] != -1)
+ {
+ // Combine info from shell and surface
+ label compactI = cellToCompact[nei];
+ FixedList<label, 3> gapInfo;
+ volumeType gapMode;
+ mergeGapInfo
+ (
+ shellGapInfo[compactI],
+ shellGapMode[compactI],
+ extendedGapLevel[globalRegionI],
+ extendedGapMode[globalRegionI],
+ gapInfo,
+ gapMode
+ );
+
+ const point& cc = cellCentres[nei];
+ bool okRay = generateRay
+ (
+ false,
+ surfPt,
+ normal1[i],
+ gapInfo,
+ gapMode,
+ surfPt+((cc-surfPt)&normal1[i])*normal1[i],
+ cellLevel[nei],
+
+ start[compactI],
+ end[compactI],
+ start2[compactI],
+ end2[compactI]
+ );
+ cellMap[compactI] = nei;
+
+ if (okRay)
+ {
+ map[compactI++] = i;
+ }
+ }
+ }
+ else
+ {
+ label bFaceI = faceI - mesh_.nInternalFaces();
+
+ if (bFaceToCompact[bFaceI] != -1)
+ {
+ // Combine info from shell and surface
+ label compactI = bFaceToCompact[bFaceI];
+ FixedList<label, 3> gapInfo;
+ volumeType gapMode;
+ mergeGapInfo
+ (
+ shellGapInfo[compactI],
+ shellGapMode[compactI],
+ extendedGapLevel[globalRegionI],
+ extendedGapMode[globalRegionI],
+ gapInfo,
+ gapMode
+ );
+
+ const point& cc = neiCc[bFaceI];
+ bool okRay = generateRay
+ (
+ false,
+ surfPt,
+ normal1[i],
+ gapInfo,
+ gapMode,
+ surfPt+((cc-surfPt)&normal1[i])*normal1[i],
+ neiLevel[bFaceI],
+
+ start[compactI],
+ end[compactI],
+ start2[compactI],
+ end2[compactI]
+ );
+ cellMap[compactI] = -1;
+
+ if (okRay)
+ {
+ map[compactI++] = i;
+ }
+ }
+ }
+ }
+ }
+
+ //Info<< "Retesting " << returnReduce(compactI, sumOp<label>())
+ // << " out of " << returnReduce(start.size(), sumOp<label>())
+ // << endl;
+
+ start.setSize(compactI);
+ end.setSize(compactI);
+ start2.setSize(compactI);
+ end2.setSize(compactI);
+ map.setSize(compactI);
+ cellMap.setSize(compactI);
+
+ testFaces = UIndirectList<label>(testFaces, map)();
+ minLevel = UIndirectList<label>(minLevel, map)();
+ surface1.clear();
+ hit1.clear();
+ region1.clear();
+ normal1 = UIndirectList<vector>(normal1, map)();
+
+ // Do intersections in first direction
+ labelList surfaceHit;
+ vectorField surfaceNormal;
+ surfaces_.findNearestIntersection
+ (
+ start,
+ end,
+ surfaceHit,
+ surfaceNormal
+ );
+ {
+ // Do intersections in second direction and merge
+ labelList surfaceHit2;
+ vectorField surfaceNormal2;
+ surfaces_.findNearestIntersection
+ (
+ start2,
+ end2,
+ surfaceHit2,
+ surfaceNormal2
+ );
+ forAll(surfaceHit, i)
+ {
+ if (surfaceHit[i] == -1 && surfaceHit2[i] != -1)
+ {
+ surfaceHit[i] = surfaceHit2[i];
+ surfaceNormal[i] = surfaceNormal2[i];
+ }
+ }
+ }
+
+
+ forAll(surfaceHit, i)
+ {
+ label surfI = surfaceHit[i];
+
+ if (surfI != -1)
+ {
+ // Found intersection with surface. Check opposite normal.
+
+ label cellI = cellMap[i];
+
+ if (cellI != -1 && mag(normal1[i]&surfaceNormal[i]) > planarCos)
+ {
+ if
+ (
+ !markForRefine
+ (
+ surfI,
+ nAllowRefine,
+ refineCell[cellI],
+ nRefine
+ )
+ )
+ {
+ break;
+ }
+ }
+ }
+ }
+
+ if
+ (
+ returnReduce(nRefine, sumOp<label>())
+ > returnReduce(nAllowRefine, sumOp<label>())
+ )
+ {
+ Info<< "Reached refinement limit." << endl;
+ }
+ }
+
+ return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
+//Foam::labelList Foam::meshRefinement::extractHits
+//(
+// const List<pointIndexHit>& info
+//)
+//{
+// labelList compactMap(info.size());
+// label compactI = 0;
+// forAll(info, i)
+// {
+// if (info[i].hit())
+// {
+// compactMap[compactI++] = i;
+// }
+// }
+// compactMap.setSize(compactI);
+// return compactMap;
+//}
+//
+//
+//void Foam::meshRefinement::mergeHits
+//(
+// const pointField& samples,
+// const List<pointIndexHit>& extraInfo,
+// List<pointIndexHit>& info
+//)
+//{
+// if (extraInfo.size() != info.size())
+// {
+// FatalErrorIn
+// (
+// "meshRefinement::mergeHits(const List<pointIndexHit>&"
+// ", List<pointIndexHit>& info)"
+// ) << "Sizes differ " << extraInfo.size() << ' ' << info.size()
+// << exit(FatalError);
+// }
+//
+// // Merge oppositeInfo2 into oppositeInfo
+// forAll(extraInfo, i)
+// {
+// if (!info[i].hit())
+// {
+// if (extraInfo[i].hit())
+// {
+// info[i] = extraInfo[i];
+// }
+// }
+// else if (extraInfo[i].hit())
+// {
+// const point& nearPt = samples[i]; //nearInfo[i].hitPoint();
+//
+// scalar d = magSqr(info[i].hitPoint()-nearPt);
+// scalar d2 = magSqr(extraInfo[i].hitPoint()-nearPt);
+//
+// if (d2 < d)
+// {
+// info[i] = extraInfo[i];
+// }
+// }
+// }
+//}
+//
+//
+//Foam::tmp<Foam::pointField> Foam::meshRefinement::extractPoints
+//(
+// const List<pointIndexHit>& info
+//)
+//{
+// tmp<pointField> tfld(new pointField(info.size()));
+// pointField& fld = tfld();
+//
+// forAll(info, i)
+// {
+// fld[i] = info[i].rawPoint();
+// }
+// return tfld;
+//}
+//
+//
+//Foam::meshRefinement::findNearestOppositeOp::findNearestOppositeOp
+//(
+// const indexedOctree<treeDataTriSurface>& tree,
+// const point& oppositePoint,
+// const vector& oppositeNormal,
+// const scalar minCos
+//)
+//:
+// tree_(tree),
+// oppositePoint_(oppositePoint),
+// oppositeNormal_(oppositeNormal),
+// minCos_(minCos)
+//{}
+//
+//
+//void Foam::meshRefinement::findNearestOppositeOp::operator()
+//(
+// const labelUList& indices,
+// const point& sample,
+// scalar& nearestDistSqr,
+// label& minIndex,
+// point& nearestPoint
+//) const
+//{
+// const treeDataTriSurface& shape = tree_.shapes();
+// const triSurface& patch = shape.patch();
+// const pointField& points = patch.points();
+//
+// forAll(indices, i)
+// {
+// const label index = indices[i];
+// const labelledTri& f = patch[index];
+//
+// pointHit nearHit = f.nearestPoint(sample, points);
+// scalar distSqr = sqr(nearHit.distance());
+//
+// if (distSqr < nearestDistSqr)
+// {
+// // Nearer. Check if
+// // - a bit way from other hit
+// // - in correct search cone
+// vector d(nearHit.rawPoint()-oppositePoint_);
+// scalar normalDist(d&oppositeNormal_);
+//
+// if (normalDist > Foam::sqr(SMALL) && normalDist/mag(d) > minCos_)
+// {
+// nearestDistSqr = distSqr;
+// minIndex = index;
+// nearestPoint = nearHit.rawPoint();
+// }
+// }
+// }
+//}
+//
+//
+//void Foam::meshRefinement::searchCone
+//(
+// const label surfI,
+// labelList& nearMap, // cells
+// scalarField& nearGap, // gap size
+// List<pointIndexHit>& nearInfo, // nearest point on surface
+// List<pointIndexHit>& oppositeInfo // detected point on gap (or miss)
+//) const
+//{
+// const labelList& cellLevel = meshCutter_.cellLevel();
+// const pointField& cellCentres = mesh_.cellCentres();
+// const scalar edge0Len = meshCutter_.level0EdgeLength();
+//
+// const labelList& surfaceIndices = surfaces_.surfaces();
+// const List<FixedList<label, 3> >& extendedGapLevel =
+// surfaces_.extendedGapLevel();
+// const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+//
+//
+// label geomI = surfaceIndices[surfI];
+// const searchableSurface& geom = surfaces_.geometry()[geomI];
+//
+// const triSurfaceMesh& s = refCast<const triSurfaceMesh>(geom);
+// const indexedOctree<treeDataTriSurface>& tree = s.tree();
+//
+//
+// const scalar searchCos(Foam::cos(degToRad(30)));
+//
+// // Normals for ray shooting and inside/outside detection
+// vectorField nearNormal;
+// geom.getNormal(nearInfo, nearNormal);
+// // Regions
+// labelList nearRegion;
+// geom.getRegion(nearInfo, nearRegion);
+//
+//
+// // Now loop over all near points and search in the half cone
+// labelList map(nearInfo.size());
+// label compactI = 0;
+//
+// oppositeInfo.setSize(nearInfo.size());
+//
+// forAll(nearInfo, i)
+// {
+// label globalRegionI =
+// surfaces_.globalRegion(surfI, nearRegion[i]);
+//
+// // Get updated gap information now we have the region
+// label nGapCells = extendedGapLevel[globalRegionI][0];
+// label minLevel = extendedGapLevel[globalRegionI][1];
+// label maxLevel = extendedGapLevel[globalRegionI][2];
+// volumeType mode = extendedGapMode[globalRegionI];
+//
+// label cellI = nearMap[i];
+// label cLevel = cellLevel[cellI];
+//
+// if (cLevel >= minLevel && cLevel < maxLevel)
+// {
+// scalar cellSize = edge0Len/pow(2, cLevel);
+//
+// // Update gap size
+// nearGap[i] = nGapCells*cellSize;
+//
+// const point& nearPt = nearInfo[i].hitPoint();
+// vector v(cellCentres[cellI]-nearPt);
+// scalar magV = mag(v);
+//
+// // Like with ray shooting we want to
+// // - find triangles up to nearGap away on the wanted side of the
+// // surface
+// // - find triangles up to 0.5*cellSize away on the unwanted side
+// // of the surface. This is for cells straddling the surface
+// // where
+// // the cell centre might be on the wrong side of the surface
+//
+// // Tbd: check that cell centre is inbetween the gap hits
+// // (only if the cell is far enough away)
+//
+// scalar posNormalSize = 0.0;
+// scalar negNormalSize = 0.0;
+//
+// if (mode == volumeType::OUTSIDE)
+// {
+// posNormalSize = nearGap[i];
+// if (magV < 0.5*cellSize)
+// {
+// negNormalSize = 0.5*cellSize;
+// }
+// }
+// else if (mode == volumeType::INSIDE)
+// {
+// if (magV < 0.5*cellSize)
+// {
+// posNormalSize = 0.5*cellSize;
+// }
+// negNormalSize = nearGap[i];
+// }
+// else
+// {
+// posNormalSize = nearGap[i];
+// negNormalSize = nearGap[i];
+// }
+//
+// // Test with positive normal
+// oppositeInfo[compactI] = tree.findNearest
+// (
+// nearPt,
+// sqr(posNormalSize),
+// findNearestOppositeOp
+// (
+// tree,
+// nearPt,
+// nearNormal[i],
+// searchCos
+// )
+// );
+//
+// if (oppositeInfo[compactI].hit())
+// {
+// map[compactI++] = i;
+// }
+// else
+// {
+// // Test with negative normal
+// oppositeInfo[compactI] = tree.findNearest
+// (
+// nearPt,
+// sqr(negNormalSize),
+// findNearestOppositeOp
+// (
+// tree,
+// nearPt,
+// -nearNormal[i],
+// searchCos
+// )
+// );
+//
+// if (oppositeInfo[compactI].hit())
+// {
+// map[compactI++] = i;
+// }
+// }
+// }
+// }
+//
+// Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+// << " hits on the correct side out of "
+// << returnReduce(map.size(), sumOp<label>()) << endl;
+// map.setSize(compactI);
+// oppositeInfo.setSize(compactI);
+//
+// nearMap = UIndirectList<label>(nearMap, map)();
+// nearGap = UIndirectList<scalar>(nearGap, map)();
+// nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
+// nearNormal = UIndirectList<vector>(nearNormal, map)();
+//
+// // Exclude hits which aren't opposite enough. E.g. you might find
+// // a point on a perpendicular wall - but this does not consistute a gap.
+// vectorField oppositeNormal;
+// geom.getNormal(oppositeInfo, oppositeNormal);
+//
+// compactI = 0;
+// forAll(oppositeInfo, i)
+// {
+// if ((nearNormal[i] & oppositeNormal[i]) < -0.707)
+// {
+// map[compactI++] = i;
+// }
+// }
+//
+// Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+// << " hits opposite the nearest out of "
+// << returnReduce(map.size(), sumOp<label>()) << endl;
+// map.setSize(compactI);
+//
+// nearMap = UIndirectList<label>(nearMap, map)();
+// nearGap = UIndirectList<scalar>(nearGap, map)();
+// nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
+// oppositeInfo = UIndirectList<pointIndexHit>(oppositeInfo, map)();
+//}
+
+
+bool Foam::meshRefinement::generateRay
+(
+ const bool useSurfaceNormal,
+
+ const point& nearPoint,
+ const vector& nearNormal,
+ const FixedList<label, 3>& gapInfo,
+ const volumeType& mode,
+
+ const point& cc,
+ const label cLevel,
+
+ point& start,
+ point& end,
+ point& start2,
+ point& end2
+) const
+{
+ // We want to handle the following cases:
+ // - surface: small gap (marked with 'surface'). gap might be
+ // on inside or outside of surface.
+ // - A: cell well inside the gap.
+ // - B: cell well outside the gap.
+ // - C: cell straddling the gap. cell centre might be inside
+ // or outside
+ //
+ // +---+
+ // | B |
+ // +---+
+ //
+ // +------+
+ // | |
+ // | C |
+ // --------|------|----surface
+ // +------+
+ //
+ // +---+
+ // | A |
+ // +---+
+ //
+ //
+ // --------------------surface
+ //
+ // So:
+ // - find nearest point on surface
+ // - in situation A,B decide if on wanted side of surface
+ // - detect if locally a gap (and the cell inside the gap) by
+ // shooting a ray from the point on the surface in the direction
+ // of
+ // - A,B: the cell centre
+ // - C: the surface normal and/or negative surface normal
+ // and see we hit anything
+ //
+ // Variations of this scheme:
+ // - always shoot in the direction of the surface normal. This needs
+ // then an additional check to make sure the cell centre is
+ // somewhere inside the gap
+ // - instead of ray shooting use a 'constrained' nearest search
+ // by e.g. looking inside a search cone (implemented in searchCone).
+ // The problem with this constrained nearest is that it still uses
+ // the absolute nearest point on each triangle and only afterwards
+ // checks if it is inside the search cone.
+
+
+ // Decide which near points are good:
+ // - with updated minLevel and maxLevel and nearGap make sure
+ // the cell is still a candidate
+ // NOTE: inside the gap the nearest point on the surface will
+ // be HALF the gap size - otherwise we would have found
+ // a point on the opposite side
+ // - if the mode is both sides
+ // - or if the hit is inside the current cell (situation 'C',
+ // magV < 0.5cellSize)
+ // - or otherwise if on the correct side
+
+ bool okRay = false;
+
+ if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
+ {
+ scalar cellSize = meshCutter_.level0EdgeLength()/pow(2, cLevel);
+
+ // Calculate gap size
+ scalar nearGap = gapInfo[0]*cellSize;
+
+ // Distance to nearest
+ vector v(cc-nearPoint);
+ scalar magV = mag(v);
+
+ if (useSurfaceNormal || magV < 0.5*cellSize)
+ {
+ const vector& n = nearNormal;
+
+ // Situation 'C' above: cell too close. Use surface
+ // normal to shoot rays
+
+ if (mode == volumeType::OUTSIDE)
+ {
+ start = nearPoint+1e-6*n;
+ end = nearPoint+nearGap*n;
+ // Small second vector just to make sure to
+ // refine cell
+ start2 = nearPoint-1e-6*n;
+ end2 = nearPoint-0.5*cellSize*n;
+ }
+ else if (mode == volumeType::INSIDE)
+ {
+ start = nearPoint-1e-6*n;
+ end = nearPoint-nearGap*n;
+ // Small second vector just to make sure to
+ // refine cell
+ start2 = nearPoint+1e-6*n;
+ end2 = nearPoint+0.5*cellSize*n;
+ }
+ else if (mode == volumeType::MIXED)
+ {
+ start = nearPoint+1e-6*n;
+ end = nearPoint+nearGap*n;
+
+ start2 = nearPoint-1e-6*n;
+ end2 = nearPoint-nearGap*n;
+ }
+
+ okRay = true;
+ }
+ else
+ {
+ // Situation 'A' or 'B' above: cell well away. Test if
+ // cell on correct side of surface and shoot ray through
+ // cell centre. Note: no need to shoot ray in other
+ // direction since we're trying to detect cell inside
+ // the gap.
+
+ scalar s = (v&nearNormal);
+
+ if
+ (
+ (mode == volumeType::MIXED)
+ || (mode == volumeType::OUTSIDE && s > SMALL)
+ || (mode == volumeType::INSIDE && s < -SMALL)
+ )
+ {
+ // Use vector through cell centre
+ vector n(v/(magV+ROOTVSMALL));
+
+ start = nearPoint+1e-6*n;
+ end = nearPoint+nearGap*n;
+ // Dummy second vector
+ start2 = start;
+ end2 = start2;
+
+ okRay = true;
+ }
+ }
+ }
+
+ return okRay;
+}
+
+
+//void Foam::meshRefinement::shootRays
+//(
+// const label surfI,
+// labelList& nearMap,
+// scalarField& nearGap,
+// List<pointIndexHit>& nearInfo,
+// List<pointIndexHit>& oppositeInfo
+//) const
+//{
+// const labelList& cellLevel = meshCutter_.cellLevel();
+// const scalar edge0Len = meshCutter_.level0EdgeLength();
+//
+// const labelList& surfaceIndices = surfaces_.surfaces();
+// const List<FixedList<label, 3> >& extendedGapLevel =
+// surfaces_.extendedGapLevel();
+// const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+//
+//
+// label geomI = surfaceIndices[surfI];
+// const searchableSurface& geom = surfaces_.geometry()[geomI];
+//
+//
+// // Normals for ray shooting and inside/outside detection
+// vectorField nearNormal;
+// geom.getNormal(nearInfo, nearNormal);
+// // Regions
+// labelList nearRegion;
+// geom.getRegion(nearInfo, nearRegion);
+//
+// labelList map(nearInfo.size());
+//
+// pointField start(nearInfo.size());
+// pointField end(nearInfo.size());
+// pointField start2(nearInfo.size());
+// pointField end2(nearInfo.size());
+//
+// label compactI = 0;
+// forAll(nearInfo, i)
+// {
+// label globalRegionI =
+// surfaces_.globalRegion(surfI, nearRegion[i]);
+//
+// label cellI = nearMap[i];
+// label cLevel = cellLevel[cellI];
+// scalar cellSize = edge0Len/pow(2, cLevel);
+//
+// // Update gap size
+// nearGap[i] = extendedGapLevel[globalRegionI][0]*cellSize;
+//
+// // Construct one or two rays to test for oppositeness
+// bool okRay = generateRay
+// (
+// false,
+// nearInfo[i].hitPoint(),
+// nearNormal[i],
+// extendedGapLevel[globalRegionI],
+// extendedGapMode[globalRegionI],
+//
+// cellCentres[cellI],
+// cellLevel[cellI],
+//
+// start[compactI],
+// end[compactI],
+// start2[compactI],
+// end2[compactI]
+// );
+// if (okRay)
+// {
+// map[compactI++] = i;
+// }
+// }
+//
+// Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+// << " hits on the correct side out of "
+// << returnReduce(map.size(), sumOp<label>()) << endl;
+// map.setSize(compactI);
+// start.setSize(compactI);
+// end.setSize(compactI);
+// start2.setSize(compactI);
+// end2.setSize(compactI);
+//
+// nearMap = UIndirectList<label>(nearMap, map)();
+// nearGap = UIndirectList<scalar>(nearGap, map)();
+// nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
+//
+// geom.findLineAny(start, end, oppositeInfo);
+//
+// List<pointIndexHit> oppositeInfo2;
+// geom.findLineAny(start2, end2, oppositeInfo2);
+// mergeHits(extractPoints(nearInfo), oppositeInfo2, oppositeInfo);
+//}
+
+
+void Foam::meshRefinement::selectGapCandidates
+(
+ const labelList& refineCell,
+ const label nRefine,
+
+ labelList& cellMap,
+ List<FixedList<label, 3> >& shellGapInfo,
+ List<volumeType>& shellGapMode
+) const
+{
+ const labelList& cellLevel = meshCutter_.cellLevel();
+ const pointField& cellCentres = mesh_.cellCentres();
+
+ // Collect cells to test
+ cellMap.setSize(cellLevel.size()-nRefine);
+ label compactI = 0;
+
+ forAll(cellLevel, cellI)
+ {
+ if (refineCell[cellI] == -1)
+ {
+ cellMap[compactI++] = cellI;
+ }
+ }
+ Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+ << " unmarked cells out of "
+ << mesh_.globalData().nTotalCells() << endl;
+ cellMap.setSize(compactI);
+
+ // Do test to see whether cells are inside/outside shell with
+ // applicable specification (minLevel <= celllevel < maxLevel)
+ shells_.findHigherGapLevel
+ (
+ pointField(cellCentres, cellMap),
+ UIndirectList<label>(cellLevel, cellMap)(),
+ shellGapInfo,
+ shellGapMode
+ );
+
+ // Compact out hits
+
+ labelList map(shellGapInfo.size());
+ compactI = 0;
+ forAll(shellGapInfo, i)
+ {
+ if (shellGapInfo[i][2] > 0)
+ {
+ map[compactI++] = i;
+ }
+ }
+
+ Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+ << " cells inside gap shells out of "
+ << mesh_.globalData().nTotalCells() << endl;
+
+ map.setSize(compactI);
+ cellMap = UIndirectList<label>(cellMap, map)();
+ shellGapInfo = UIndirectList<FixedList<label, 3> >(shellGapInfo, map)();
+ shellGapMode = UIndirectList<volumeType>(shellGapMode, map)();
+}
+
+
+void Foam::meshRefinement::mergeGapInfo
+(
+ const FixedList<label, 3>& shellGapInfo,
+ const volumeType shellGapMode,
+ const FixedList<label, 3>& surfGapInfo,
+ const volumeType surfGapMode,
+ FixedList<label, 3>& gapInfo,
+ volumeType& gapMode
+) const
+{
+ if (surfGapInfo[0] == 0)
+ {
+ gapInfo = shellGapInfo;
+ gapMode = shellGapMode;
+ }
+ else if (shellGapInfo[0] == 0)
+ {
+ gapInfo = surfGapInfo;
+ gapMode = surfGapMode;
+ }
+ else
+ {
+ // Both specify a level. Does surface level win? Or does information
+ // need to be merged?
+
+ //gapInfo[0] = max(surfGapInfo[0], shellGapInfo[0]);
+ //gapInfo[1] = min(surfGapInfo[1], shellGapInfo[1]);
+ //gapInfo[2] = max(surfGapInfo[2], shellGapInfo[2]);
+ gapInfo = surfGapInfo;
+ gapMode = surfGapMode;
+ }
+}
+
+
+Foam::label Foam::meshRefinement::markInternalGapRefinement
+(
+ const scalar planarCos,
+ const label nAllowRefine,
+
+ labelList& refineCell,
+ label& nRefine
+) const
+{
+ const labelList& cellLevel = meshCutter_.cellLevel();
+ const pointField& cellCentres = mesh_.cellCentres();
+ const scalar edge0Len = meshCutter_.level0EdgeLength();
+
+ const List<FixedList<label, 3> >& extendedGapLevel =
+ surfaces_.extendedGapLevel();
+ const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+
+ // Get the gap level for the shells
+ const labelList maxLevel(shells_.maxGapLevel());
+
+ label oldNRefine = nRefine;
+
+ if (max(maxLevel) > 0)
+ {
+ // Collect cells to test
+ labelList cellMap;
+ List<FixedList<label, 3> > shellGapInfo;
+ List<volumeType> shellGapMode;
+ selectGapCandidates
+ (
+ refineCell,
+ nRefine,
+
+ cellMap,
+ shellGapInfo,
+ shellGapMode
+ );
+
+
+ // Find nearest point and normal on the surfaces
+ List<pointIndexHit> nearInfo;
+ vectorField nearNormal;
+ labelList nearSurface;
+ labelList nearRegion;
+ {
+ // Now we have both the cell-level and the gap size information. Use
+ // this to calculate the gap size
+ scalarField gapSize(cellMap.size());
+ forAll(cellMap, i)
+ {
+ label cellI = cellMap[i];
+ scalar cellSize = edge0Len/pow(2, cellLevel[cellI]);
+ gapSize[i] = shellGapInfo[i][0]*cellSize;
+ }
+
+ surfaces_.findNearestRegion
+ (
+ identity(surfaces_.surfaces().size()),
+ pointField(cellCentres, cellMap),
+ sqr(gapSize),
+ nearSurface,
+ nearInfo,
+ nearRegion,
+ nearNormal
+ );
+ }
+
+
+
+ labelList map(nearInfo.size());
+ pointField start(nearInfo.size());
+ pointField end(nearInfo.size());
+ pointField start2(nearInfo.size());
+ pointField end2(nearInfo.size());
+
+ label compactI = 0;
+
+ forAll(nearInfo, i)
+ {
+ if (nearInfo[i].hit())
+ {
+ label globalRegionI = surfaces_.globalRegion
+ (
+ nearSurface[i],
+ nearRegion[i]
+ );
+
+ // Combine info from shell and surface
+ FixedList<label, 3> gapInfo;
+ volumeType gapMode;
+ mergeGapInfo
+ (
+ shellGapInfo[i],
+ shellGapMode[i],
+ extendedGapLevel[globalRegionI],
+ extendedGapMode[globalRegionI],
+ gapInfo,
+ gapMode
+ );
+
+
+ // Construct one or two rays to test for oppositeness
+ bool okRay = generateRay
+ (
+ false,
+ nearInfo[i].hitPoint(),
+ nearNormal[i],
+ gapInfo,
+ gapMode,
+
+ cellCentres[cellMap[i]],
+ cellLevel[cellMap[i]],
+
+ start[compactI],
+ end[compactI],
+ start2[compactI],
+ end2[compactI]
+ );
+ if (okRay)
+ {
+ map[compactI++] = i;
+ }
+ }
+ }
+
+ Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+ << " cells for testing out of "
+ << mesh_.globalData().nTotalCells() << endl;
+ map.setSize(compactI);
+ start.setSize(compactI);
+ end.setSize(compactI);
+ start2.setSize(compactI);
+ end2.setSize(compactI);
+
+ cellMap = UIndirectList<label>(cellMap, map)();
+ nearNormal = UIndirectList<vector>(nearNormal, map)();
+ shellGapInfo.clear();
+ shellGapMode.clear();
+ nearInfo.clear();
+ nearSurface.clear();
+ nearRegion.clear();
+
+
+ // Do intersections in first direction
+ labelList surfaceHit;
+ vectorField surfaceNormal;
+ surfaces_.findNearestIntersection
+ (
+ start,
+ end,
+ surfaceHit,
+ surfaceNormal
+ );
+ {
+ // Do intersections in second direction and merge
+ labelList surfaceHit2;
+ vectorField surfaceNormal2;
+ surfaces_.findNearestIntersection
+ (
+ start2,
+ end2,
+ surfaceHit2,
+ surfaceNormal2
+ );
+ forAll(surfaceHit, i)
+ {
+ if (surfaceHit[i] == -1 && surfaceHit2[i] != -1)
+ {
+ surfaceHit[i] = surfaceHit2[i];
+ surfaceNormal[i] = surfaceNormal2[i];
+ }
+ }
+ }
+
+
+ forAll(surfaceHit, i)
+ {
+ label surfI = surfaceHit[i];
+
+ if (surfI != -1 && mag(nearNormal[i]&surfaceNormal[i]) > planarCos)
+ {
+ // Found intersection with surface. Check opposite normal.
+
+ label cellI = cellMap[i];
+
+ if
+ (
+ !markForRefine
+ (
+ surfI,
+ nAllowRefine,
+ refineCell[cellI],
+ nRefine
+ )
+ )
+ {
+ break;
+ }
+ }
+ }
+
+ if
+ (
+ returnReduce(nRefine, sumOp<label>())
+ > returnReduce(nAllowRefine, sumOp<label>())
+ )
+ {
+ Info<< "Reached refinement limit." << endl;
+ }
+ }
+
+ return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
+Foam::label Foam::meshRefinement::markSmallFeatureRefinement
+(
+ 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();
+
+ const labelList& surfaceIndices = surfaces_.surfaces();
+ const List<FixedList<label, 3> >& extendedGapLevel =
+ surfaces_.extendedGapLevel();
+ const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
+
+ label oldNRefine = nRefine;
+
+ // Check that we're using any gap refinement
+ labelList shellMaxLevel(shells_.maxGapLevel());
+
+ if (max(shellMaxLevel) == 0)
+ {
+ return 0;
+ }
+
+ //- Force calculation of tetBasePt
+ (void)mesh_.tetBasePtIs();
+ (void)mesh_.cellTree();
+
+
+ forAll(surfaceIndices, surfI)
+ {
+ label geomI = surfaceIndices[surfI];
+ const searchableSurface& geom = surfaces_.geometry()[geomI];
+
+
+ // Get the element index in a roundabout way. Problem is e.g.
+ // distributed surface where local indices differ from global
+ // ones (needed for getRegion call)
+
+ pointField ctrs;
+ labelList region;
+ vectorField normal;
+ {
+ // Representative local coordinates and bounding sphere
+ scalarField radiusSqr;
+ geom.boundingSpheres(ctrs, radiusSqr);
+
+ List<pointIndexHit> info;
+ geom.findNearest(ctrs, radiusSqr, info);
+
+ forAll(info, i)
+ {
+ if (!info[i].hit())
+ {
+ FatalErrorIn("meshRefinement::markSmallFeatureRefinement")
+ << "fc:" << ctrs[i]
+ << " radius:" << radiusSqr[i]
+ << exit(FatalError);
+ }
+ }
+
+ geom.getRegion(info, region);
+ geom.getNormal(info, normal);
+ }
+
+ // Do test to see whether triangles are inside/outside shell with
+ // applicable specification (minLevel <= celllevel < maxLevel)
+ List<FixedList<label, 3> > shellGapInfo;
+ List<volumeType> shellGapMode;
+ shells_.findHigherGapLevel
+ (
+ ctrs,
+ labelList(ctrs.size(), 0),
+ shellGapInfo,
+ shellGapMode
+ );
+
+
+ labelList map(ctrs.size());
+ labelList cellMap(ctrs.size());
+ label compactI = 0;
+
+ pointField start(ctrs.size());
+ pointField end(ctrs.size());
+ pointField start2(ctrs.size());
+ pointField end2(ctrs.size());
+
+ forAll(ctrs, i)
+ {
+ if (shellGapInfo[i][2] > 0)
+ {
+ label globalRegionI = surfaces_.globalRegion(surfI, region[i]);
+
+ // Combine info from shell and surface
+ FixedList<label, 3> gapInfo;
+ volumeType gapMode;
+ mergeGapInfo
+ (
+ shellGapInfo[i],
+ shellGapMode[i],
+ extendedGapLevel[globalRegionI],
+ extendedGapMode[globalRegionI],
+ gapInfo,
+ gapMode
+ );
+
+ //- Option 1: use octree nearest searching inside polyMesh
+ //label cellI = mesh_.findCell(pt);
+
+ //- Option 2: use octree 'inside' searching inside polyMesh. Is
+ // much faster.
+ label cellI = -1;
+ const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
+ if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
+ {
+ cellI = tree.findInside(ctrs[i]);
+ }
+
+ if (cellI != -1 && refineCell[cellI] == -1)
+ {
+ // Construct one or two rays to test for oppositeness
+ // Note that we always want to use the surface normal
+ // and not the vector from cell centre to surface point
+
+ bool okRay = generateRay
+ (
+ true, // always use surface normal
+ ctrs[i],
+ normal[i],
+ gapInfo,
+ gapMode,
+
+ cellCentres[cellI],
+ cellLevel[cellI],
+
+ start[compactI],
+ end[compactI],
+ start2[compactI],
+ end2[compactI]
+ );
+ if (okRay)
+ {
+ cellMap[compactI] = cellI;
+ map[compactI++] = i;
+ }
+ }
+ }
+ }
+
+ Info<< "Selected " << returnReduce(compactI, sumOp<label>())
+ << " cells containing triangle centres out of "
+ << mesh_.globalData().nTotalCells() << endl;
+ map.setSize(compactI);
+ cellMap.setSize(compactI);
+ start.setSize(compactI);
+ end.setSize(compactI);
+ start2.setSize(compactI);
+ end2.setSize(compactI);
+
+ ctrs.clear();
+ region.clear();
+ shellGapInfo.clear();
+ shellGapMode.clear();
+ normal = UIndirectList<vector>(normal, map)();
+
+
+ // Do intersections in first direction.
+ // Note passing in of -1 for cell level such that it does not
+ // discard surfaces with level 0. (since can be overridden with
+ // gap level)
+ labelList surfaceHit;
+ vectorField surfaceNormal;
+ surfaces_.findNearestIntersection
+ (
+ start,
+ end,
+ surfaceHit,
+ surfaceNormal
+ );
+ {
+ // Do intersections in second direction and merge
+ labelList surfaceHit2;
+ vectorField surfaceNormal2;
+ surfaces_.findNearestIntersection
+ (
+ start2,
+ end2,
+ surfaceHit2,
+ surfaceNormal2
+ );
+ forAll(surfaceHit, i)
+ {
+ if (surfaceHit[i] == -1 && surfaceHit2[i] != -1)
+ {
+ surfaceHit[i] = surfaceHit2[i];
+ surfaceNormal[i] = surfaceNormal2[i];
+ }
+ }
+ }
+
+
+ label nGaps = 0;
+
+ forAll(surfaceHit, i)
+ {
+ label surfI = surfaceHit[i];
+
+ if (surfI != -1 && mag(normal[i]&surfaceNormal[i]) > planarCos)
+ {
+ nGaps++;
+
+ if
+ (
+ !markForRefine
+ (
+ surfI,
+ nAllowRefine,
+ refineCell[cellMap[i]],
+ nRefine
+ )
+ )
+ {
+ break;
+ }
+ }
+ }
+
+ Info<< "For surface " << geom.name() << " found "
+ << returnReduce(nGaps, sumOp<label>())
+ << " cells in small gaps" << endl;
+
+ if
+ (
+ returnReduce(nRefine, sumOp<label>())
+ > returnReduce(nAllowRefine, sumOp<label>())
+ )
+ {
+ Info<< "Reached refinement limit." << endl;
+ }
+ }
+
+ return returnReduce(nRefine-oldNRefine, sumOp<label>());
+}
+
+
+// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index b27c4fae3f80a1b06c78b56993663662bb6acee8..8429a3c1175334cae79260499d718e9277757dba 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -742,7 +742,7 @@ Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
}
}
- // Do test to see whether cells is inside/outside shell with higher level
+ // Do test to see whether cells are inside/outside shell with higher level
labelList maxLevel;
features_.findHigherLevel(testCc, testLevels, maxLevel);
@@ -820,7 +820,7 @@ Foam::label Foam::meshRefinement::markInternalRefinement
}
}
- // Do test to see whether cells is inside/outside shell with higher level
+ // Do test to see whether cells are inside/outside shell with higher level
labelList maxLevel;
shells_.findHigherLevel(testCc, testLevels, maxLevel);
@@ -869,6 +869,56 @@ Foam::label Foam::meshRefinement::markInternalRefinement
}
+Foam::label Foam::meshRefinement::unmarkInternalRefinement
+(
+ labelList& refineCell,
+ label& nRefine
+) const
+{
+ const labelList& cellLevel = meshCutter_.cellLevel();
+ const pointField& cellCentres = mesh_.cellCentres();
+
+ label oldNRefine = nRefine;
+
+ // Collect cells to test
+ pointField testCc(nRefine);
+ labelList testLevels(nRefine);
+ label testI = 0;
+
+ forAll(cellLevel, cellI)
+ {
+ if (refineCell[cellI] >= 0)
+ {
+ testCc[testI] = cellCentres[cellI];
+ testLevels[testI] = cellLevel[cellI];
+ testI++;
+ }
+ }
+
+ // Do test to see whether cells are inside/outside shell with higher level
+ labelList levelShell;
+ limitShells_.findLevel(testCc, testLevels, levelShell);
+
+ // Mark for refinement. Note that we didn't store the original cellID so
+ // now just reloop in same order.
+ testI = 0;
+ forAll(cellLevel, cellI)
+ {
+ if (refineCell[cellI] >= 0)
+ {
+ if (levelShell[testI] != -1)
+ {
+ refineCell[cellI] = -1;
+ nRefine--;
+ }
+ testI++;
+ }
+ }
+
+ return returnReduce(oldNRefine-nRefine, sumOp<label>());
+}
+
+
// Collect faces that are intersected and whose neighbours aren't yet marked
// for refinement.
Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
@@ -922,7 +972,6 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
- const pointField& cellCentres = mesh_.cellCentres();
label oldNRefine = nRefine;
@@ -941,36 +990,15 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
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(ROOTSMALL*(end-start));
- start -= smallVec;
- end += smallVec;
- }
+ calcCellCellRays
+ (
+ neiCc,
+ neiLevel,
+ testFaces,
+ start,
+ end,
+ minLevel
+ );
// Do test for higher intersections
@@ -1130,6 +1158,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
pointField end(testFaces.size());
labelList minLevel(testFaces.size());
+ // Note: uses isMasterFace otherwise could be call to calcCellCellRays
forAll(testFaces, i)
{
label faceI = testFaces[i];
@@ -1719,7 +1748,6 @@ Foam::label Foam::meshRefinement::markProximityRefinement
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
- const pointField& cellCentres = mesh_.cellCentres();
label oldNRefine = nRefine;
@@ -1739,36 +1767,15 @@ Foam::label Foam::meshRefinement::markProximityRefinement
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(ROOTSMALL*(end-start));
- start -= smallVec;
- end += smallVec;
- }
+ calcCellCellRays
+ (
+ neiCc,
+ neiLevel,
+ testFaces,
+ start,
+ end,
+ minLevel
+ );
// Test for all intersections (with surfaces of higher gap level than
@@ -2045,7 +2052,9 @@ Foam::labelList Foam::meshRefinement::refineCandidates
const bool internalRefinement,
const bool surfaceRefinement,
const bool curvatureRefinement,
+ const bool smallFeatureRefinement,
const bool gapRefinement,
+ const bool bigGapRefinement,
const label maxGlobalCells,
const label maxLocalCells
) const
@@ -2096,6 +2105,8 @@ Foam::labelList Foam::meshRefinement::refineCandidates
calcNeighbourData(neiLevel, neiCc);
+ const scalar planarCos = Foam::cos(degToRad(planarAngle));
+
// Cells pierced by feature lines
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -2163,6 +2174,30 @@ Foam::labelList Foam::meshRefinement::refineCandidates
);
Info<< "Marked for refinement due to surface intersection "
<< ": " << nSurf << " cells." << endl;
+
+ // Refine intersected-cells only inside gaps. See
+ // markInternalGapRefinement to refine all cells inside gaps.
+ if
+ (
+ planarCos >= -1
+ && planarCos <= 1
+ && max(shells_.maxGapLevel()) > 0
+ )
+ {
+ label nGapSurf = markSurfaceGapRefinement
+ (
+ planarCos,
+ nAllowRefine,
+ neiLevel,
+ neiCc,
+
+ refineCell,
+ nRefine
+ );
+ Info<< "Marked for refinement due to surface intersection"
+ << " (at gaps)"
+ << ": " << nGapSurf << " cells." << endl;
+ }
}
// Refinement based on curvature of surface
@@ -2190,7 +2225,33 @@ Foam::labelList Foam::meshRefinement::refineCandidates
}
- const scalar planarCos = Foam::cos(degToRad(planarAngle));
+ // Refinement based on features smaller than cell size
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ if
+ (
+ smallFeatureRefinement
+ && (planarCos >= -1 && planarCos <= 1)
+ && max(shells_.maxGapLevel()) > 0
+ )
+ {
+ label nGap = markSmallFeatureRefinement
+ (
+ planarCos,
+ nAllowRefine,
+ neiLevel,
+ neiCc,
+
+ refineCell,
+ nRefine
+ );
+ Info<< "Marked for refinement due to close opposite surfaces "
+ << ": " << nGap << " cells." << endl;
+ }
+
+
+ // Refinement based on gap (only neighbouring cells)
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if
(
@@ -2217,6 +2278,44 @@ Foam::labelList Foam::meshRefinement::refineCandidates
}
+ // Refinement based on gaps larger than cell size
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ if
+ (
+ bigGapRefinement
+ && (planarCos >= -1 && planarCos <= 1)
+ && max(shells_.maxGapLevel()) > 0
+ )
+ {
+ // Refine based on gap information provided by shell and nearest
+ // surface
+ label nGap = markInternalGapRefinement
+ (
+ planarCos,
+ nAllowRefine,
+ refineCell,
+ nRefine
+ );
+ Info<< "Marked for refinement due to opposite surfaces"
+ << " "
+ << ": " << nGap << " cells." << endl;
+ }
+
+
+ // Limit refinement
+ // ~~~~~~~~~~~~~~~~
+
+ {
+ label nUnmarked = unmarkInternalRefinement(refineCell, nRefine);
+ if (nUnmarked > 0)
+ {
+ Info<< "Unmarked for refinement due to limit shells"
+ << " : " << nUnmarked << " cells." << endl;
+ }
+ }
+
+
// Pack cells-to-refine
// ~~~~~~~~~~~~~~~~~~~~
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index 7d47a612c9dc35a655594fdfe6a449b7fb490264..6f36ce59041f6d1dcecee87e8cd2b4bc3cbd3751 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -157,11 +157,22 @@ Foam::refinementSurfaces::refinementSurfaces
labelList globalMinLevel(surfI, 0);
labelList globalMaxLevel(surfI, 0);
labelList globalLevelIncr(surfI, 0);
+
+ FixedList<label, 3> nullGapLevel;
+ nullGapLevel[0] = 0;
+ nullGapLevel[1] = 0;
+ nullGapLevel[2] = 0;
+
+ List<FixedList<label, 3> > globalGapLevel(surfI);
+ List<volumeType> globalGapMode(surfI);
+
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<FixedList<label, 3> > > regionGapLevel(surfI);
+ List<Map<volumeType> > regionGapMode(surfI);
List<Map<scalar> > regionAngle(surfI);
List<Map<autoPtr<dictionary> > > regionPatchInfo(surfI);
@@ -212,6 +223,44 @@ Foam::refinementSurfaces::refinementSurfaces
<< exit(FatalIOError);
}
+
+ // Optional gapLevel specification
+
+ globalGapLevel[surfI] = dict.lookupOrDefault
+ (
+ "gapLevel",
+ nullGapLevel
+ );
+ globalGapMode[surfI] = volumeType::names
+ [
+ dict.lookupOrDefault<word>
+ (
+ "gapMode",
+ volumeType::names[volumeType::MIXED]
+ )
+ ];
+ if
+ (
+ globalGapMode[surfI] == volumeType::UNKNOWN
+ || globalGapLevel[surfI][0] < 0
+ || globalGapLevel[surfI][1] < 0
+ || globalGapLevel[surfI][2] < 0
+ || globalGapLevel[surfI][1] > globalGapLevel[surfI][2]
+ )
+ {
+ FatalIOErrorIn
+ (
+ "refinementSurfaces::refinementSurfaces"
+ "(const searchableSurfaces&, const dictionary>&",
+ dict
+ ) << "Illegal gapLevel specification for surface "
+ << names_[surfI]
+ << " : gapLevel:" << globalGapLevel[surfI]
+ << " gapMode:" << volumeType::names[globalGapMode[surfI]]
+ << exit(FatalIOError);
+ }
+
+
const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
// Surface zones
@@ -275,6 +324,54 @@ Foam::refinementSurfaces::refinementSurfaces
<< exit(FatalIOError);
}
+
+
+ // Optional gapLevel specification
+
+ FixedList<label, 3> gapSpec
+ (
+ regionDict.lookupOrDefault
+ (
+ "gapLevel",
+ nullGapLevel
+ )
+ );
+ regionGapLevel[surfI].insert(regionI, gapSpec);
+ volumeType gapModeSpec
+ (
+ volumeType::names
+ [
+ regionDict.lookupOrDefault<word>
+ (
+ "gapMode",
+ volumeType::names[volumeType::MIXED]
+ )
+ ]
+ );
+ regionGapMode[surfI].insert(regionI, gapModeSpec);
+ if
+ (
+ gapModeSpec == volumeType::UNKNOWN
+ || gapSpec[0] < 0
+ || gapSpec[1] < 0
+ || gapSpec[2] < 0
+ || gapSpec[1] > gapSpec[2]
+ )
+ {
+ FatalIOErrorIn
+ (
+ "refinementSurfaces::refinementSurfaces"
+ "(const searchableSurfaces&,"
+ " const dictionary>&",
+ dict
+ ) << "Illegal gapLevel specification for surface "
+ << names_[surfI]
+ << " : gapLevel:" << gapSpec
+ << " gapMode:" << volumeType::names[gapModeSpec]
+ << exit(FatalIOError);
+ }
+
+
if (regionDict.found("perpendicularAngle"))
{
regionAngle[surfI].insert
@@ -331,6 +428,10 @@ Foam::refinementSurfaces::refinementSurfaces
maxLevel_ = 0;
gapLevel_.setSize(nRegions);
gapLevel_ = -1;
+ extendedGapLevel_.setSize(nRegions);
+ extendedGapLevel_ = nullGapLevel;
+ extendedGapMode_.setSize(nRegions);
+ extendedGapMode_ = volumeType::UNKNOWN;
perpendicularAngle_.setSize(nRegions);
perpendicularAngle_ = -GREAT;
patchInfo_.setSize(nRegions);
@@ -349,7 +450,8 @@ Foam::refinementSurfaces::refinementSurfaces
gapLevel_[globalRegionI] =
maxLevel_[globalRegionI]
+ globalLevelIncr[surfI];
-
+ extendedGapLevel_[globalRegionI] = globalGapLevel[surfI];
+ extendedGapMode_[globalRegionI] = globalGapMode[surfI];
perpendicularAngle_[globalRegionI] = globalAngle[surfI];
if (globalPatchInfo.set(surfI))
{
@@ -371,6 +473,10 @@ Foam::refinementSurfaces::refinementSurfaces
gapLevel_[globalRegionI] =
maxLevel_[globalRegionI]
+ regionLevelIncr[surfI][iter.key()];
+ extendedGapLevel_[globalRegionI] =
+ regionGapLevel[surfI][iter.key()];
+ extendedGapMode_[globalRegionI] =
+ regionGapMode[surfI][iter.key()];
}
forAllConstIter(Map<scalar>, regionAngle[surfI], iter)
{
@@ -482,12 +588,28 @@ Foam::refinementSurfaces::refinementSurfaces
// }
+Foam::labelList Foam::refinementSurfaces::maxGapLevel() const
+{
+ labelList surfaceMax(surfaces_.size(), 0);
+
+ forAll(surfaces_, surfI)
+ {
+ const wordList& regionNames = allGeometry_[surfaces_[surfI]].regions();
+
+ forAll(regionNames, regionI)
+ {
+ label globalI = globalRegion(surfI, regionI);
+ const FixedList<label, 3>& gapInfo = extendedGapLevel_[globalI];
+ surfaceMax[surfI] = max(surfaceMax[surfI], gapInfo[2]);
+ }
+ }
+ return surfaceMax;
+}
+
+
// Precalculate the refinement level for every element of the searchable
// surface.
-void Foam::refinementSurfaces::setMinLevelFields
-(
- const shellSurfaces& shells
-)
+void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
{
forAll(surfaces_, surfI)
{
@@ -1208,6 +1330,49 @@ void Foam::refinementSurfaces::findNearestIntersection
}
+void Foam::refinementSurfaces::findNearestIntersection
+(
+ const pointField& start,
+ const pointField& end,
+
+ labelList& surface1,
+ vectorField& normal1
+) const
+{
+ // Initialize arguments
+ surface1.setSize(start.size());
+ surface1 = -1;
+ normal1.setSize(start.size());
+ normal1 = vector::zero;
+
+ // Current end of segment to test.
+ pointField nearest(end);
+ // Work array
+ List<pointIndexHit> nearestInfo(start.size());
+ labelList region;
+ vectorField normal;
+
+ forAll(surfaces_, surfI)
+ {
+ const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
+
+ // See if any intersection between start and current nearest
+ geom.findLine(start, nearest, nearestInfo);
+ geom.getNormal(nearestInfo, normal);
+
+ forAll(nearestInfo, pointI)
+ {
+ if (nearestInfo[pointI].hit())
+ {
+ surface1[pointI] = surfI;
+ normal1[pointI] = normal[pointI];
+ nearest[pointI] = nearestInfo[pointI].hitPoint();
+ }
+ }
+ }
+}
+
+
void Foam::refinementSurfaces::findAnyIntersection
(
const pointField& start,
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index fadb9aaa2694e6b1d1b3dfa8c02e138a1c642307..72f20c17773bdd7d1637c0d9975acf14dd10b524 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -42,6 +42,7 @@ SourceFiles
#include "vectorList.H"
#include "pointIndexHit.H"
#include "surfaceZonesInfo.H"
+#include "volumeType.H"
#include "pointList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -85,6 +86,12 @@ class refinementSurfaces
//- From global region number to small-gap level
labelList gapLevel_;
+ //- From global region number to small-gap level specification
+ List<FixedList<label, 3> > extendedGapLevel_;
+
+ //- From global region number to side of surface to detect
+ List<volumeType> extendedGapMode_;
+
//- From global region number to perpendicular angle
scalarField perpendicularAngle_;
@@ -188,6 +195,23 @@ public:
return gapLevel_;
}
+ //- From global region number to specification of gap and its
+ // refinement: 3 labels specifying
+ // - minimum wanted number of cells in the gap
+ // - minimum cell level when to start trying to detect gaps
+ // - maximum cell level to refine to (so do not detect gaps if
+ // cell >= maximum level)
+ const List<FixedList<label, 3> >& extendedGapLevel() const
+ {
+ return extendedGapLevel_;
+ }
+
+ //- From global region number to side of surface to detect
+ const List<volumeType>& extendedGapMode() const
+ {
+ return extendedGapMode_;
+ }
+
//- From global region number to perpendicular angle
const scalarField& perpendicularAngle() const
{
@@ -226,6 +250,9 @@ public:
return minLevel_.size();
}
+ //- Per surface the maximum extendedGapLevel over all its regions
+ labelList maxGapLevel() const;
+
//- Calculate minLevelFields
void setMinLevelFields
(
@@ -314,6 +341,15 @@ public:
vectorField& normal2
) const;
+ //- Find nearest (to start only) intersection of edge
+ void findNearestIntersection
+ (
+ const pointField& start,
+ const pointField& end,
+ labelList& surfaces,
+ vectorField& normal
+ ) const;
+
//- Used for debugging only: find intersection of edge.
void findAnyIntersection
(
diff --git a/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
index cd398025534b2394e07b2790afd40f725aa65b83..128d2892ec25a9bc85e2f53519f4d35878d68fd1 100644
--- a/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -124,7 +124,7 @@ void Foam::shellSurfaces::setAndCheckLevels
}
else
{
- if (!allGeometry_[shells_[shellI]].hasVolumeType())
+ if (!shell.hasVolumeType())
{
FatalErrorIn
(
@@ -152,6 +152,46 @@ void Foam::shellSurfaces::setAndCheckLevels
}
+void Foam::shellSurfaces::checkGapLevels
+(
+ const dictionary& shellDict,
+ const label shellI,
+ const List<FixedList<label, 3> >& levels
+)
+{
+ const searchableSurface& shell = allGeometry_[shells_[shellI]];
+
+ forAll(levels, regionI)
+ {
+ const FixedList<label, 3>& info = levels[regionI];
+
+ if (info[2] > 0)
+ {
+ if (modes_[shellI] == DISTANCE)
+ {
+ FatalIOErrorIn
+ (
+ "shellSurfaces::shellSurfaces(..)",
+ shellDict
+ ) << "'gapLevel' specification cannot be used with mode "
+ << refineModeNames_[DISTANCE]
+ << " for shell " << shell.name()
+ << exit(FatalIOError);
+ }
+ }
+ }
+
+ // Hardcode for region 0
+ if (levels[0][0] > 0)
+ {
+ Info<< "Refinement level up to " << levels[0][2]
+ << " for all cells in gaps for shell "
+ << shell.name() << endl;
+ }
+}
+
+
+
// Specifically orient triSurfaces using a calculated point outside.
// Done since quite often triSurfaces not of consistent orientation which
// is (currently) necessary for sideness calculation
@@ -288,18 +328,18 @@ void Foam::shellSurfaces::findHigherLevel
);
// Update maxLevel
- forAll(nearInfo, candidateI)
+ forAll(nearInfo, i)
{
- if (nearInfo[candidateI].hit())
+ if (nearInfo[i].hit())
{
// Check which level it actually is in.
label minDistI = findLower
(
distances,
- mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
+ mag(nearInfo[i].hitPoint()-candidates[i])
);
- label pointI = candidateMap[candidateI];
+ label pointI = candidateMap[i];
// pt is inbetween shell[minDistI] and shell[minDistI+1]
maxLevel[pointI] = levels[minDistI+1];
@@ -356,6 +396,197 @@ void Foam::shellSurfaces::findHigherLevel
}
+void Foam::shellSurfaces::findHigherGapLevel
+(
+ const pointField& pt,
+ const labelList& ptLevel,
+ const label shellI,
+ labelList& gapShell,
+ List<FixedList<label, 3> >& gapInfo,
+ List<volumeType>& gapMode
+) const
+{
+ //TBD: hardcoded for region 0 information
+ const FixedList<label, 3>& info = extendedGapLevel_[shellI][0];
+ volumeType mode = extendedGapMode_[shellI][0];
+
+ if (info[2] == 0)
+ {
+ return;
+ }
+
+
+ // Collect all those points that have a current maxLevel less than the
+ // shell.
+
+ labelList candidateMap(pt.size());
+ label candidateI = 0;
+
+ forAll(ptLevel, pointI)
+ {
+ if (ptLevel[pointI] >= info[1] && ptLevel[pointI] < info[2])
+ {
+ candidateMap[candidateI++] = pointI;
+ }
+ }
+ candidateMap.setSize(candidateI);
+
+ // Do the expensive nearest test only for the candidate points.
+ List<volumeType> volType;
+ allGeometry_[shells_[shellI]].getVolumeType
+ (
+ pointField(pt, candidateMap),
+ volType
+ );
+
+ forAll(volType, i)
+ {
+ label pointI = candidateMap[i];
+
+ bool isInside = (volType[i] == volumeType::INSIDE);
+
+ if
+ (
+ (
+ (modes_[shellI] == INSIDE && isInside)
+ || (modes_[shellI] == OUTSIDE && !isInside)
+ )
+ && info[2] > gapInfo[pointI][2]
+ )
+ {
+ gapShell[pointI] = shellI;
+ gapInfo[pointI] = info;
+ gapMode[pointI] = mode;
+ }
+ }
+}
+
+
+void Foam::shellSurfaces::findLevel
+(
+ const pointField& pt,
+ const label shellI,
+ labelList& minLevel,
+ labelList& shell
+) const
+{
+ const labelList& levels = levels_[shellI];
+
+ if (modes_[shellI] == DISTANCE)
+ {
+ // Distance mode.
+
+ const scalarField& distances = distances_[shellI];
+
+ // Collect all those points that have a current level equal/greater
+ // (any of) the shell. Also collect the furthest distance allowable
+ // to any shell with a higher level.
+
+ pointField candidates(pt.size());
+ labelList candidateMap(pt.size());
+ scalarField candidateDistSqr(pt.size());
+ label candidateI = 0;
+
+ forAll(shell, pointI)
+ {
+ if (shell[pointI] == -1)
+ {
+ forAllReverse(levels, levelI)
+ {
+ if (levels[levelI] <= minLevel[pointI])
+ {
+ candidates[candidateI] = pt[pointI];
+ candidateMap[candidateI] = pointI;
+ candidateDistSqr[candidateI] = sqr(distances[levelI]);
+ candidateI++;
+ break;
+ }
+ }
+ }
+ }
+ candidates.setSize(candidateI);
+ candidateMap.setSize(candidateI);
+ candidateDistSqr.setSize(candidateI);
+
+ // Do the expensive nearest test only for the candidate points.
+ List<pointIndexHit> nearInfo;
+ allGeometry_[shells_[shellI]].findNearest
+ (
+ candidates,
+ candidateDistSqr,
+ nearInfo
+ );
+
+ // Update maxLevel
+ forAll(nearInfo, i)
+ {
+ if (nearInfo[i].hit())
+ {
+ // Check which level it actually is in.
+ label minDistI = findLower
+ (
+ distances,
+ mag(nearInfo[i].hitPoint()-candidates[i])
+ );
+
+ label pointI = candidateMap[i];
+
+ // pt is inbetween shell[minDistI] and shell[minDistI+1]
+ shell[pointI] = shellI;
+ minLevel[pointI] = levels[minDistI+1];
+ }
+ }
+ }
+ else
+ {
+ // Inside/outside mode
+
+ // Collect all those points that have a current maxLevel less than the
+ // shell.
+
+ pointField candidates(pt.size());
+ labelList candidateMap(pt.size());
+ label candidateI = 0;
+
+ forAll(shell, pointI)
+ {
+ if (shell[pointI] == -1 && levels[0] <= minLevel[pointI])
+ {
+ candidates[candidateI] = pt[pointI];
+ candidateMap[candidateI] = pointI;
+ candidateI++;
+ }
+ }
+ candidates.setSize(candidateI);
+ candidateMap.setSize(candidateI);
+
+ // Do the expensive nearest test only for the candidate points.
+ List<volumeType> volType;
+ allGeometry_[shells_[shellI]].getVolumeType(candidates, volType);
+
+ forAll(volType, i)
+ {
+ if
+ (
+ (
+ modes_[shellI] == INSIDE
+ && volType[i] == volumeType::INSIDE
+ )
+ || (
+ modes_[shellI] == OUTSIDE
+ && volType[i] == volumeType::OUTSIDE
+ )
+ )
+ {
+ label pointI = candidateMap[i];
+ shell[pointI] = shellI;
+ minLevel[pointI] = levels[0];
+ }
+ }
+ }
+}
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::shellSurfaces::shellSurfaces
@@ -387,6 +618,15 @@ Foam::shellSurfaces::shellSurfaces
distances_.setSize(shellI);
levels_.setSize(shellI);
+ extendedGapLevel_.setSize(shellI);
+ extendedGapMode_.setSize(shellI);
+
+ FixedList<label, 3> nullGapLevel;
+ nullGapLevel[0] = 0;
+ nullGapLevel[1] = 0;
+ nullGapLevel[2] = 0;
+
+
HashSet<word> unmatchedKeys(shellsDict.toc());
shellI = 0;
@@ -407,6 +647,85 @@ Foam::shellSurfaces::shellSurfaces
// Read pairs of distance+level
setAndCheckLevels(shellI, dict.lookup("levels"));
+
+
+ // Gap specification
+ // ~~~~~~~~~~~~~~~~~
+
+
+ // Shell-wide gap level specification
+ const searchableSurface& surface = allGeometry_[geomI];
+ const wordList& regionNames = surface.regions();
+
+ FixedList<label, 3> gapSpec
+ (
+ dict.lookupOrDefault
+ (
+ "gapLevel",
+ nullGapLevel
+ )
+ );
+ extendedGapLevel_[shellI].setSize(regionNames.size());
+ extendedGapLevel_[shellI] = gapSpec;
+
+ volumeType gapModeSpec
+ (
+ volumeType::names
+ [
+ dict.lookupOrDefault<word>
+ (
+ "gapMode",
+ volumeType::names[volumeType::MIXED]
+ )
+ ]
+ );
+ extendedGapMode_[shellI].setSize(regionNames.size());
+ extendedGapMode_[shellI] = gapModeSpec;
+
+
+ // Override on a per-region basis?
+
+ if (dict.found("regions"))
+ {
+ const dictionary& regionsDict = dict.subDict("regions");
+ forAll(regionNames, regionI)
+ {
+ if (regionsDict.found(regionNames[regionI]))
+ {
+ // Get the dictionary for region
+ const dictionary& regionDict = regionsDict.subDict
+ (
+ regionNames[regionI]
+ );
+ FixedList<label, 3> gapSpec
+ (
+ regionDict.lookupOrDefault
+ (
+ "gapLevel",
+ nullGapLevel
+ )
+ );
+ extendedGapLevel_[shellI][regionI] = gapSpec;
+
+ volumeType gapModeSpec
+ (
+ volumeType::names
+ [
+ regionDict.lookupOrDefault<word>
+ (
+ "gapMode",
+ volumeType::names[volumeType::MIXED]
+ )
+ ]
+ );
+ extendedGapMode_[shellI][regionI] = gapModeSpec;
+ }
+ }
+ }
+
+
+ checkGapLevels(dict, shellI, extendedGapLevel_[shellI]);
+
shellI++;
}
}
@@ -444,6 +763,22 @@ Foam::label Foam::shellSurfaces::maxLevel() const
}
+Foam::labelList Foam::shellSurfaces::maxGapLevel() const
+{
+ labelList surfaceMax(extendedGapLevel_.size(), 0);
+
+ forAll(extendedGapLevel_, shellI)
+ {
+ const List<FixedList<label, 3> >& levels = extendedGapLevel_[shellI];
+ forAll(levels, i)
+ {
+ surfaceMax[shellI] = max(surfaceMax[shellI], levels[i][2]);
+ }
+ }
+ return surfaceMax;
+}
+
+
void Foam::shellSurfaces::findHigherLevel
(
const pointField& pt,
@@ -461,4 +796,66 @@ void Foam::shellSurfaces::findHigherLevel
}
+void Foam::shellSurfaces::findHigherGapLevel
+(
+ const pointField& pt,
+ const labelList& ptLevel,
+ labelList& gapShell,
+ List<FixedList<label, 3> >& gapInfo,
+ List<volumeType>& gapMode
+) const
+{
+ gapShell.setSize(pt.size());
+ gapShell = -1;
+
+ FixedList<label, 3> nullGapLevel;
+ nullGapLevel[0] = 0;
+ nullGapLevel[1] = 0;
+ nullGapLevel[2] = 0;
+
+ gapInfo.setSize(pt.size());
+ gapInfo = nullGapLevel;
+
+ gapMode.setSize(pt.size());
+ gapMode = volumeType::MIXED;
+
+ forAll(shells_, shellI)
+ {
+ findHigherGapLevel(pt, ptLevel, shellI, gapShell, gapInfo, gapMode);
+ }
+}
+
+
+void Foam::shellSurfaces::findHigherGapLevel
+(
+ const pointField& pt,
+ const labelList& ptLevel,
+ List<FixedList<label, 3> >& gapInfo,
+ List<volumeType>& gapMode
+) const
+{
+ labelList gapShell;
+ findHigherGapLevel(pt, ptLevel, gapShell, gapInfo, gapMode);
+}
+
+
+void Foam::shellSurfaces::findLevel
+(
+ const pointField& pt,
+ const labelList& ptLevel,
+ labelList& shell
+) const
+{
+ shell.setSize(pt.size());
+ shell = -1;
+
+ labelList minLevel(ptLevel);
+
+ forAll(shells_, shellI)
+ {
+ findLevel(pt, shellI, minLevel, shell);
+ }
+}
+
+
// ************************************************************************* //
diff --git a/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H b/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H
index 9b0e4c9e050be2dc3c43e9f84da3489ebb2f283e..13817de16f0b10d124f0ee90d54d904d0d5c33ba 100644
--- a/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H
+++ b/src/mesh/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H
@@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
- \\/ M anipulation |
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -85,6 +85,15 @@ private:
labelListList levels_;
+ // Gap level refinement
+
+ //- Per shell, per region the small-gap level specification
+ List<List<FixedList<label, 3> > > extendedGapLevel_;
+
+ //- Per shell, per region the small-gap level specification
+ List<List<volumeType> > extendedGapMode_;
+
+
// Private data
//- refineMode names
@@ -93,15 +102,24 @@ private:
// Private Member Functions
- //- Helper function for initialisation.
+ //- Helper function for initialisation of levels
void setAndCheckLevels
(
const label shellI,
const List<Tuple2<scalar, label> >&
);
+ //- Helper function for checking of gap information
+ void checkGapLevels
+ (
+ const dictionary&,
+ const label shellI,
+ const List<FixedList<label, 3> >& levels
+ );
+
void orient();
+ //- Find first shell with a level higher than maxLevel
void findHigherLevel
(
const pointField& pt,
@@ -109,6 +127,27 @@ private:
labelList& maxLevel
) const;
+ //- Update highest min gap level
+ void findHigherGapLevel
+ (
+ const pointField& pt,
+ const labelList& ptLevel,
+ const label shellI,
+ labelList& gapShell,
+ List<FixedList<label, 3> >& gapInfo,
+ List<volumeType>& gapMode
+ ) const;
+
+ //- Find first shell with a level lower or equal to minLevel. Update
+ // minLevel and shell.
+ void findLevel
+ (
+ const pointField& pt,
+ const label shellI,
+ labelList& minLevel,
+ labelList& shell
+ ) const;
+
public:
// Constructors
@@ -125,16 +164,11 @@ public:
// Access
- //const List<scalarField>& distances() const
- //{
- // return distances_;
- //}
- //
- ////- Per shell per distance the refinement level
- //const labelListList& levels() const
- //{
- // return levels_;
- //}
+ //- Indices of surfaces that are shells
+ const labelList& shells() const
+ {
+ return shells_;
+ }
// Query
@@ -142,6 +176,9 @@ public:
//- Highest shell level
label maxLevel() const;
+ //- Highest shell gap level
+ labelList maxGapLevel() const;
+
//- Find shell level higher than ptLevel
void findHigherLevel
(
@@ -150,6 +187,33 @@ public:
labelList& maxLevel
) const;
+ //- Find a shell whose minimum gap level is >= ptLevel
+ void findHigherGapLevel
+ (
+ const pointField& pt,
+ const labelList& ptLevel,
+ labelList& gapShell,
+ List<FixedList<label, 3> >& gapInfo,
+ List<volumeType>& gapMode
+ ) const;
+
+ //- Find a shell whose minimum gap level is >= ptLevel. gapInfo
+ // is (0 0 0) if no shell found
+ void findHigherGapLevel
+ (
+ const pointField& pt,
+ const labelList& ptLevel,
+ List<FixedList<label, 3> >& gapInfo,
+ List<volumeType>& gapMode
+ ) const;
+
+ //- Find first shell (or -1) with level equal or lower than ptLevel.
+ void findLevel
+ (
+ const pointField& pt,
+ const labelList& ptLevel,
+ labelList& shell
+ ) const;
};
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 0c1b794f7c6d45d76d37d2448d8cc5b8e6725278..bd0abc24a9be4c307f0476f391e3ef3183699781 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -63,6 +63,7 @@ indexedOctree/treeDataTriSurface.C
searchableSurface = searchableSurface
$(searchableSurface)/searchableBox.C
+$(searchableSurface)/searchableRotatedBox.C
$(searchableSurface)/searchableCylinder.C
$(searchableSurface)/searchableDisk.C
$(searchableSurface)/searchablePlane.C
@@ -75,6 +76,7 @@ $(searchableSurface)/searchableSurfacesQueries.C
$(searchableSurface)/searchableSurfaceWithGaps.C
$(searchableSurface)/triSurfaceMesh.C
$(searchableSurface)/closedTriSurfaceMesh.C
+$(searchableSurface)/subTriSurfaceMesh.C
topoSets = sets/topoSets
$(topoSets)/cellSet.C
diff --git a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
index 9928f39f4738c01dec1f2042c995c108fbf20c19..0d19394afdc2bcd98ea6b4242671f634a3cdfcbb 100644
--- a/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
+++ b/src/parallel/decompose/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.C
@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
- \\/ M anipulation |
+ \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -408,7 +408,14 @@ void Foam::hierarchGeomDecomp::sortComponent
{
// No need for binary searching of bin size
localSize = label(current.size()/n_[compI]);
- rightCoord = sortedCoord[leftIndex+localSize];
+ if (leftIndex+localSize < sortedCoord.size())
+ {
+ rightCoord = sortedCoord[leftIndex+localSize];
+ }
+ else
+ {
+ rightCoord = maxCoord;
+ }
}
else
{