Commit e1665232 authored by Franjo's avatar Franjo
Browse files

Faster refinement thickness

parent 862c2cf6
......@@ -66,38 +66,6 @@ meshOctreeCreator::meshOctreeCreator
surfRefThickness_(mo.surface().size(), 0.0)
{}
/*
meshOctreeCreator::meshOctreeCreator
(
meshOctree& mo,
const IOdictionary& dict,
const volScalarField& localCellSize
)
:
octree_(mo),
scalingFactor_(1.0),
meshDict_(dict),
hexRefinement_(false),
globalRefLevel_(0),
surfRefLevel_(mo.surface().size(), direction(0)),
surfRefThickness_(mo.surface().size(), 0.0)
{
FatalErrorIn
(
"meshOctreeCreator::meshOctreeCreator"
"("
"meshOctree& mo,"
"const IOdictionary& dict,"
"const volScalarField& localCellSize"
")"
) << "Not implemented!" << exit(FatalError);
Info << "Constructing octree" << endl;
Info << "Finished constructing octree" << endl;
}
*/
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
meshOctreeCreator::~meshOctreeCreator()
......
......@@ -39,7 +39,6 @@ SourceFiles
#include "DynList.H"
#include "IOdictionary.H"
#include "meshOctreeModifier.H"
//#include "volFields.H"
#include "patchRefinementList.H"
#include <map>
......
......@@ -79,7 +79,7 @@ void meshOctreeCreator::refineBoundary()
//- select boxes which need to be refined
# ifdef USE_OMP
# pragma omp parallel for reduction(+ : nMarked) \
schedule(dynamic, Foam::min(20, leaves.size()/omp_get_num_threads()+1))
schedule(dynamic, 100)
# endif
forAll(leaves, leafI)
{
......@@ -115,7 +115,7 @@ void meshOctreeCreator::refineBoundary()
Foam::max
(
nLayers[leafI],
Foam::max(label(surfRefThickness_[triI]/cs), 1)
Foam::max(ceil(surfRefThickness_[triI]/cs), 1)
);
targetLevel[leafI] =
......@@ -343,7 +343,7 @@ void meshOctreeCreator::refineBoxesContainedInObjects()
Foam::max
(
nLayers[leafI],
Foam::max(label(refThickness[oI]/cs), 1)
Foam::max(ceil(refThickness[oI]/cs), 1)
);
targetRefLevel[leafI] =
......@@ -587,7 +587,7 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
Foam::max
(
nLayers[leafI],
max(label(refThickness[surfI]/cs), 1)
max(ceil(refThickness[surfI]/cs), 1)
);
targetRefLevel[leafI] =
......@@ -822,7 +822,7 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
Foam::max
(
nLayers[leafI],
max(label(refThickness[emI]/cs), 1)
max(ceil(refThickness[emI]/cs), 1)
);
targetRefLevel[leafI] =
......
......@@ -361,7 +361,10 @@ void meshOctreeCreator::setRootCubeSizeAndRefParameters()
if( patchDict.found("additionalRefinementLevels") )
{
nLevel =
readLabel(patchDict.lookup("additionalRefinementLevels"));
readLabel
(
patchDict.lookup("additionalRefinementLevels")
);
}
else if( patchDict.found("cellSize") )
{
......
......@@ -171,193 +171,99 @@ label meshOctreeModifier::markAdditionalLayers
List<direction>& targetRefLevel
) const
{
const FixedList<meshOctreeCubeCoordinates, 26>& rp =
octree_.regularityPositions_;
const LongList<meshOctreeCube*>& leaves = octree_.leaves_;
//- this is needed for parallel runs to reduce the length of messages
labelHashSet transferCoordinates;
FixedList<meshOctreeCube*, 26> neighbours;
//- sort leaves based on the number of addtional layers
label maxLevel = Foam::max(nLayers);
reduce(maxLevel, maxOp<label>());
labelList currentLayer(leaves.size());
forAll(targetRefLevel, leafI)
List<labelLongList> leavesForLayer(maxLevel+1);
forAll(nLayers, leafI)
{
currentLayer[leafI] = targetRefLevel[leafI] ? 1 : 0;
if( nLayers[leafI] < 1 )
continue;
leavesForLayer[nLayers[leafI]].append(leafI);
}
bool changed;
label i(0), nMarked(0);
do
//- set refinement flag to additional boxes marked separately
label nMarked(0);
forAllReverse(leavesForLayer, layerI)
{
++i;
changed = false;
const labelLongList& activeLeaves = leavesForLayer[layerI];
labelList newNLayers = nLayers;
labelList newCurrentLayer = currentLayer;
List<direction> newTargetRefLevel = targetRefLevel;
if( returnReduce(activeLeaves.size(), sumOp<label>()) == 0 )
continue;
LongList<meshOctreeCubeCoordinates> processorChecks;
labelLongList processorNLayers;
labelLongList processorTargetLevels;
//- find the max required refinement level
direction maxLevel(0);
forAll(leaves, leafI)
# ifdef USE_OMP
# pragma omp parallel
# endif
{
if( currentLayer[leafI] != i )
continue;
const meshOctreeCube* oc = leaves[leafI];
direction localMax(0);
forAll(rp, posI)
{
const meshOctreeCubeCoordinates cc
(
oc->coordinates() + rp[posI]
);
# ifdef USE_OMP
# pragma omp for schedule(dynamic, 50)
# endif
forAll(activeLeaves, i)
localMax = Foam::max(localMax, targetRefLevel[activeLeaves[i]]);
const label neiLabel = octree_.findLeafLabelForPosition(cc);
# ifdef USE_OMP
# pragma omp critical
# endif
maxLevel = Foam::max(localMax, maxLevel);
}
if( neiLabel > -1 )
{
neighbours[posI] = leaves[neiLabel];
}
else if( neiLabel == -1 )
{
neighbours[posI] = NULL;
}
else if( neiLabel == meshOctreeCubeBasic::OTHERPROC )
{
neighbours[posI] = NULL;
reduce(maxLevel, maxOp<direction>());
if( !transferCoordinates.found(leafI) )
{
processorChecks.append(oc->coordinates());
processorNLayers.append(nLayers[leafI]);
processorTargetLevels.append(targetRefLevel[leafI]);
transferCoordinates.insert(leafI);
}
}
}
//- mark additional boxes for refinement
for(direction levelI=maxLevel;levelI>0;--levelI)
{
List<direction> markedBoxes(leaves.size(), direction(0));
forAll(neighbours, neiI)
label counter(0);
# ifdef USE_OMP
# pragma omp parallel for reduction(+:counter)
# endif
forAll(activeLeaves, lI)
{
const meshOctreeCube* nei = neighbours[neiI];
if( !nei ) continue;
if( !nei->isLeaf() ) continue;
const label leafI = activeLeaves[lI];
if( i <= nLayers[leafI] )
if( refineBox[leafI] && (targetRefLevel[leafI] == levelI) )
{
newNLayers[nei->cubeLabel()] =
Foam::max(nLayers[nei->cubeLabel()], nLayers[leafI]);
newTargetRefLevel[nei->cubeLabel()] =
Foam::max
(
targetRefLevel[nei->cubeLabel()],
targetRefLevel[leafI]
);
changed = true;
newCurrentLayer[nei->cubeLabel()] = i+1;
markedBoxes[leafI] = 1;
++counter;
}
}
}
if( octree_.neiProcs().size() )
{
std::map<label, LongList<meshOctreeCubeCoordinates> > eCoords;
std::map<label, labelLongList > eNLayers;
std::map<label, labelLongList > eTargetLevels;
forAll(octree_.neiProcs(), neiI)
{
const label neiProc = octree_.neiProcs()[neiI];
eCoords[neiProc] = processorChecks;
eNLayers[neiProc] = processorNLayers;
eTargetLevels[neiProc] = processorTargetLevels;
}
LongList<meshOctreeCubeCoordinates> receivedCoords;
help::exchangeMap(eCoords, receivedCoords);
labelLongList receivedNLayers;
help::exchangeMap(eNLayers, receivedNLayers);
if( returnReduce(counter, sumOp<label>()) == 0 )
continue;
labelLongList receivedTargetLevels;
help::exchangeMap(eTargetLevels, receivedTargetLevels);
//- mark additional cells at this refinement level
markAdditionalLayers(markedBoxes, direction(layerI));
if( receivedCoords.size() != receivedNLayers.size() )
FatalErrorIn
(
"label meshOctreeModifier::markAdditionalLayers("
"List<direction>&, labelList&, List<direction>&) const"
) << "Invalid list size" << abort(FatalError);
if( receivedCoords.size() != receivedTargetLevels.size() )
FatalErrorIn
(
"label meshOctreeModifier::markAdditionalLayers("
"List<direction>&, labelList&, List<direction>&) const"
) << "Invalid list size" << abort(FatalError);
//- check consistency with received cube coordinates
forAll(receivedCoords, ccI)
//- update the main list
# ifdef USE_OMP
# pragma omp parallel for schedule(dynamic, 100) \
reduction(+:nMarked)
# endif
forAll(markedBoxes, leafI)
{
forAll(rp, posI)
{
const meshOctreeCubeCoordinates cc
(
receivedCoords[ccI] + rp[posI]
);
const meshOctreeCube* nei =
octree_.findCubeForPosition(cc);
if( !nei ) continue;
if( !nei->isLeaf() ) continue;
if( !markedBoxes[leafI] )
continue;
if( leaves[leafI]->level() >= levelI )
continue;
if( i <= receivedNLayers[ccI] )
{
newNLayers[nei->cubeLabel()] =
Foam::max
(
nLayers[nei->cubeLabel()],
receivedNLayers[ccI]
);
newTargetRefLevel[nei->cubeLabel()] =
Foam::max
(
targetRefLevel[nei->cubeLabel()],
direction(receivedTargetLevels[ccI])
);
changed = true;
newCurrentLayer[nei->cubeLabel()] = i+1;
}
if( !refineBox[leafI] )
{
refineBox[leafI] |= 1;
++nMarked;
}
}
}
nLayers.transfer(newNLayers);
currentLayer.transfer(newCurrentLayer);
targetRefLevel.transfer(newTargetRefLevel);
//- exchange information with all processors
reduce(changed, maxOp<bool>());
} while( changed );
//- update refine data
forAll(targetRefLevel, leafI)
{
if
(
!refineBox[leafI] &&
(targetRefLevel[leafI] > leaves[leafI]->level())
)
{
refineBox[leafI] = 1;
++nMarked;
}
}
return nMarked;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment