Commit 575843ce authored by Franjo's avatar Franjo
Browse files

Revised additional refinement layers

parent 87b24cd2
......@@ -136,6 +136,7 @@ void meshOctreeCreator::refineBoundary()
reduce(useNLayers, maxOp<label>());
if( useNLayers )
{
Info << "1. Marking additional layers of boxes " << endl;
nMarked +=
octreeModifier.markAdditionalLayers
(
......@@ -503,22 +504,23 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
const vector tol = SMALL * rootBox.span();
meshOctreeModifier octreeModifier(octree_);
const LongList<meshOctreeCube*>& leaves = octreeModifier.leavesAccess();
DynList<label> leavesInBox, intersectedLeaves;
DynList<label> leavesInBox;
bool changed;
do
{
# ifdef OCTREETiming
const scalar startIter = omp_get_wtime();
# endif
nMarked = 0;
changed = false;
labelList refineCubes(leaves.size(), 0);
labelList nLayers(leaves.size(), 0);
List<direction> targetRefLevel(leaves.size(), direction(0));
bool useNLayers(false);
//- select boxes which need to be refined
//- select boxes that need to be refined
forAll(surfaceMeshesPtr, surfI)
{
const triSurf& surf = surfaceMeshesPtr[surfI];
......@@ -526,13 +528,13 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
# ifdef USE_OMP
# pragma omp parallel for \
reduction( + : nMarked) schedule(dynamic, 10) \
private(leavesInBox,intersectedLeaves)
reduction( + : nMarked) schedule(dynamic, 10) private(leavesInBox)
# endif
forAll(surf, triI)
{
//- find the bounding box of the current triangle
const labelledTri& tri = surf[triI];
boundBox triBB(points[tri[0]], points[tri[0]]);
for(label pI=1;pI<3;++pI)
{
......@@ -548,7 +550,6 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
octree_.findLeavesContainedInBox(triBB, leavesInBox);
//- check which of the leaves are intersected by the triangle
intersectedLeaves.clear();
forAll(leavesInBox, i)
{
const label leafI = leavesInBox[i];
......@@ -557,8 +558,6 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
if( oc.intersectsTriangleExact(surf, rootBox, triI) )
{
intersectedLeaves.append(leafI);
if( oc.level() < refLevels[surfI] )
{
# ifdef DEBUGSearch
......@@ -567,36 +566,32 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
<< " for refinement" << endl;
# endif
if( !refineCubes[leafI] )
++nMarked;
changed = true;
refineCubes[leafI] = 1;
}
}
}
if( refThickness[surfI] > VSMALL )
{
useNLayers = true;
forAll(intersectedLeaves, i)
{
const label leafI = intersectedLeaves[i];
const meshOctreeCube& oc = *leaves[leafI];
const scalar cs = oc.size(rootBox);
nLayers[leafI] =
Foam::max
(
nLayers[leafI],
max(ceil(refThickness[surfI]/cs), 1)
);
targetRefLevel[leafI] =
Foam::max
(
targetRefLevel[leafI],
refLevels[surfI]
);
if( refThickness[surfI] > VSMALL )
{
useNLayers = true;
const scalar cs = oc.size(rootBox);
const label numLayers =
ceil(refThickness[surfI] / cs);
nLayers[leafI] =
Foam::max
(
nLayers[leafI],
Foam::max(numLayers, 1)
);
targetRefLevel[leafI] =
Foam::max
(
targetRefLevel[leafI],
refLevels[surfI]
);
}
}
}
}
}
......@@ -605,15 +600,15 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
//- mark additional boxes for refinement to achieve
//- correct refinement distance
reduce(useNLayers, maxOp<label>());
if( useNLayers )
{
nMarked +=
octreeModifier.markAdditionalLayers
(
refineCubes,
nLayers,
targetRefLevel
);
octreeModifier.markAdditionalLayers
(
refineCubes,
nLayers,
targetRefLevel
);
}
//- refine boxes
......@@ -626,8 +621,8 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
if( octree_.neiProcs().size() != 0 )
{
reduce(nMarked, sumOp<label>());
if( nMarked )
reduce(changed, maxOp<bool>());
if( changed )
{
octreeModifier.distributeLeavesToProcessors();
......@@ -645,7 +640,7 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
# endif
}
}
} while( nMarked != 0 );
} while( changed );
Info << "Finished refinement of boxes intersecting surface meshes" << endl;
}
......
......@@ -500,13 +500,6 @@ void meshOctreeCreator::createOctreeBoxes()
{
loadDistribution(true);
}
//- delete octree data which is not needed any more
// if( Pstream::parRun() )
// {
// meshOctreeModifier om(octree_);
// om.reduceMemoryConsumption();
// }
}
void meshOctreeCreator::createOctreeWithRefinedBoundary
......
......@@ -52,75 +52,60 @@ void meshOctreeModifier::markAdditionalLayers
const label nLayers
) const
{
const FixedList<meshOctreeCubeCoordinates, 26>& rp =
octree_.regularityPositions_;
const LongList<meshOctreeCube*>& leaves = octree_.leaves_;
//- this is needed for parallel runs to reduce the bandwidth
labelHashSet transferCoordinates;
FixedList<meshOctreeCube*, 26> neighbours;
DynList<label> neiLeaves;
for(label i=1;i<=nLayers;++i)
{
LongList<meshOctreeCubeCoordinates> processorChecks;
transferCoordinates.clear();
labelLongList activeLeaves;
forAll(leaves, leafI)
if( refineBox[leafI] == i )
activeLeaves.append(leafI);
# ifdef USE_OMP
# pragma omp parallel for if( leaves.size() > 1000 ) \
private(neighbours) schedule(dynamic, 20)
# pragma omp parallel for private(neiLeaves) schedule(dynamic, 20)
# endif
forAll(leaves, leafI)
forAll(activeLeaves, lI)
{
if( refineBox[leafI] != i )
continue;
const label leafI = activeLeaves[lI];
const meshOctreeCube* oc = leaves[leafI];
const meshOctreeCubeCoordinates oc = leaves[leafI]->coordinates();
forAll(rp, posI)
{
const meshOctreeCubeCoordinates cc
(
oc->coordinates() + rp[posI]
);
octree_.findAllLeafNeighbours(oc, neiLeaves);
const label neiLabel = octree_.findLeafLabelForPosition(cc);
forAll(neiLeaves, posI)
{
const label neiLabel = neiLeaves[posI];
if( neiLabel > -1 )
if( neiLabel == meshOctreeCubeBasic::OTHERPROC )
{
neighbours[posI] = leaves[neiLabel];
}
else if( neiLabel == -1 )
{
neighbours[posI] = NULL;
}
else if( neiLabel == meshOctreeCubeBasic::OTHERPROC )
{
neighbours[posI] = NULL;
# ifdef USE_OMP
# pragma omp critical
# endif
{
if( !transferCoordinates.found(leafI) )
{
processorChecks.append(oc->coordinates());
processorChecks.append(oc);
transferCoordinates.insert(leafI);
}
}
continue;
}
}
forAll(neighbours, neiI)
{
const meshOctreeCube* nei = neighbours[neiI];
if( !nei ) continue;
if( !nei->isLeaf() ) continue;
if( nei->level() > oc->level() ) continue;
if( neiLabel < 0 )
continue;
if( !refineBox[nei->cubeLabel()] )
{
refineBox[nei->cubeLabel()] = i+1;
}
if( !refineBox[neiLabel] )
refineBox[neiLabel] = i+1;
}
}
......@@ -136,28 +121,19 @@ void meshOctreeModifier::markAdditionalLayers
//- check consistency with received cube coordinates
# ifdef USE_OMP
# pragma omp parallel for if( receivedCoords.size() > 1000 ) \
schedule(dynamic, 20)
schedule(dynamic, 20) private(neiLeaves)
# endif
forAll(receivedCoords, ccI)
{
forAll(rp, posI)
{
const meshOctreeCubeCoordinates cc
(
receivedCoords[ccI] + rp[posI]
);
const meshOctreeCube* nei =
octree_.findCubeForPosition(cc);
octree_.findAllLeafNeighbours(receivedCoords[ccI], neiLeaves);
if( !nei ) continue;
if( !nei->isLeaf() ) continue;
if( nei->level() > cc.level() ) continue;
forAll(neiLeaves, posI)
{
if( neiLeaves[posI] < 0 )
continue;
if( !refineBox[nei->cubeLabel()] )
{
refineBox[nei->cubeLabel()] = i+1;
}
if( !refineBox[neiLeaves[posI]] )
refineBox[neiLeaves[posI]] = i+1;
}
}
}
......@@ -180,7 +156,7 @@ label meshOctreeModifier::markAdditionalLayers
List<labelLongList> leavesForLayer(maxLevel+1);
forAll(nLayers, leafI)
{
if( nLayers[leafI] < 1 )
if( !refineBox[leafI] || (nLayers[leafI] < 1) )
continue;
leavesForLayer[nLayers[leafI]].append(leafI);
......@@ -245,7 +221,7 @@ label meshOctreeModifier::markAdditionalLayers
continue;
//- mark additional cells at this refinement level
markAdditionalLayers(markedBoxes, direction(layerI));
markAdditionalLayers(markedBoxes, layerI);
//- update the main list
# ifdef USE_OMP
......@@ -268,6 +244,8 @@ label meshOctreeModifier::markAdditionalLayers
}
}
reduce(nMarked, sumOp<label>());
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