Skip to content

GitLab

Commit 1a8e96b1 authored by Franjo's avatar Franjo
Browse files

Updated ensureCorrectRegularity + improvements of refinement thickness

parent 6332402a
Hide whitespace changes
Inline Side-by-side
meshLibrary/utilities/octrees/meshOctree/meshOctreeCreator/meshOctreeCreatorAdjustOctreeToSurface.C
View file @ 1a8e96b1
......@@ -62,24 +62,23 @@ void meshOctreeCreator::refineBoundary()
//- refine DATA boxes to the given level
Info << "Refining boundary boxes to the given size" << endl;
label nMarked;
bool changed;
do
{
nMarked = 0;
changed = false;
# ifdef OCTREETiming
const scalar startIter = omp_get_wtime();
# endif
labelList refineCubes(leaves.size(), 0);
labelList nLayers(leaves.size(), 0);
scalarList rThickness(leaves.size(), 0.0);
List<direction> targetLevel(leaves.size(), direction(0));
bool useNLayers(false);
//- select boxes which need to be refined
# ifdef USE_OMP
# pragma omp parallel for reduction(+ : nMarked) \
schedule(dynamic, 100)
# pragma omp parallel for schedule(dynamic, 100)
# endif
forAll(leaves, leafI)
{
......@@ -95,8 +94,6 @@ void meshOctreeCreator::refineBoundary()
const VRWGraph& containedTriangles =
oc.slotPtr()->containedTriangles_;
const scalar cs = oc.size(octree_.rootBox());
bool refine(false);
forAllRow(containedTriangles, elRowI, tI)
{
......@@ -111,11 +108,11 @@ void meshOctreeCreator::refineBoundary()
{
useNLayers = true;
nLayers[leafI] =
rThickness[leafI] =
Foam::max
(
nLayers[leafI],
Foam::max(ceil(surfRefThickness_[triI]/cs), 1)
rThickness[leafI],
surfRefThickness_[triI]
);
targetLevel[leafI] =
......@@ -126,7 +123,7 @@ void meshOctreeCreator::refineBoundary()
if( refine )
{
refineCubes[leafI] = 1;
++nMarked;
changed = true;
}
}
}
......@@ -134,20 +131,21 @@ void meshOctreeCreator::refineBoundary()
//- mark additional boxes for refinement to achieve
//- correct refinement distance
reduce(useNLayers, maxOp<label>());
if( useNLayers )
reduce(changed, maxOp<bool>());
if( useNLayers && changed )
{
Info << "1. Marking additional layers of boxes " << endl;
nMarked +=
octreeModifier.markAdditionalLayers
(
refineCubes,
nLayers,
targetLevel
);
octreeModifier.refineSelectedBoxesAndAdditionalLayers
(
refineCubes,
rThickness,
targetLevel
);
}
else if( changed )
{
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
}
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
# ifdef OCTREETiming
const scalar refTime = omp_get_wtime();
......@@ -156,9 +154,7 @@ void meshOctreeCreator::refineBoundary()
if( Pstream::parRun() )
{
reduce(nMarked, sumOp<label>());
if( nMarked )
if( changed )
{
octreeModifier.distributeLeavesToProcessors();
......@@ -177,7 +173,7 @@ void meshOctreeCreator::refineBoundary()
}
}
} while( nMarked );
} while( changed );
Info << "Finished refining boundary boxes" << endl;
}
......@@ -289,24 +285,24 @@ void meshOctreeCreator::refineBoxesContainedInObjects()
meshOctreeModifier octreeModifier(octree_);
const LongList<meshOctreeCube*>& leaves = octreeModifier.leavesAccess();
label nMarked;
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);
scalarList rThickness(leaves.size(), 0.0);
List<direction> targetRefLevel(leaves.size(), direction(0));
bool useNLayers(false);
//- select boxes which need to be refined
# ifdef USE_OMP
# pragma omp parallel for if( leaves.size() > 1000 ) \
reduction( + : nMarked) schedule(dynamic, 20)
schedule(dynamic, 20)
# endif
forAll(leaves, leafI)
{
......@@ -340,13 +336,8 @@ void meshOctreeCreator::refineBoxesContainedInObjects()
if( refThickness[oI] > VSMALL )
{
const scalar cs = bb.max().x() - bb.min().x();
nLayers[leafI] =
Foam::max
(
nLayers[leafI],
Foam::max(ceil(refThickness[oI]/cs), 1)
);
rThickness[leafI] =
Foam::max(rThickness[leafI], refThickness[oI]);
targetRefLevel[leafI] =
Foam::max(targetRefLevel[leafI], refLevels[oI]);
......@@ -358,26 +349,28 @@ void meshOctreeCreator::refineBoxesContainedInObjects()
if( refine )
{
refineCubes[leafI] = 1;
++nMarked;
changed = true;
}
}
//- mark additional boxes for refinement to achieve
//- correct refinement distance
reduce(useNLayers, maxOp<label>());
if( useNLayers )
reduce(changed, maxOp<bool>());
if( useNLayers && changed )
{
nMarked +=
octreeModifier.markAdditionalLayers
(
refineCubes,
nLayers,
targetRefLevel
);
octreeModifier.refineSelectedBoxesAndAdditionalLayers
(
refineCubes,
rThickness,
targetRefLevel
);
}
else if( changed )
{
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
}
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
# ifdef OCTREETiming
const scalar refTime = omp_get_wtime();
......@@ -386,8 +379,7 @@ void meshOctreeCreator::refineBoxesContainedInObjects()
if( octree_.neiProcs().size() != 0 )
{
reduce(nMarked, sumOp<label>());
if( nMarked )
if( changed )
{
octreeModifier.distributeLeavesToProcessors();
......@@ -406,7 +398,7 @@ void meshOctreeCreator::refineBoxesContainedInObjects()
}
}
} while( nMarked != 0 );
} while( changed );
//- delete coordinate modifier if it exists
deleteDemandDrivenData(cModPtr);
......@@ -437,7 +429,7 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
scalarList refThickness(surfaces.size(), 0.0);
//- load surface meshes into memory
forAll(surfaceMeshesPtr, surfI)
forAll(surfaces, surfI)
{
const dictionary& dict = surfDict.subDict(surfaces[surfI]);
......@@ -518,6 +510,7 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
labelList refineCubes(leaves.size(), 0);
labelList nLayers(leaves.size(), 0);
List<direction> targetRefLevel(leaves.size(), direction(0));
scalarField rThickness(leaves.size(), 0.0);
bool useNLayers(false);
//- select boxes that need to be refined
......@@ -526,6 +519,9 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
const triSurf& surf = surfaceMeshesPtr[surfI];
const pointField& points = surf.points();
const direction surfLevel = refLevels[surfI];
const scalar sThickness = refThickness[surfI];
# ifdef USE_OMP
# pragma omp parallel for \
reduction( + : nMarked) schedule(dynamic, 10) private(leavesInBox)
......@@ -558,39 +554,42 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
if( oc.intersectsTriangleExact(surf, rootBox, triI) )
{
if( oc.level() < refLevels[surfI] )
{
# ifdef DEBUGSearch
Info << "Marking leaf " << leafI
<< " with coordinates " << oc
<< " for refinement" << endl;
# endif
# ifdef DEBUGSearch
Info << "Marking leaf " << leafI
<< " with coordinates " << oc
<< " for refinement" << endl;
# endif
if( oc.level() < surfLevel )
{
//- mark boxes for refinement
changed = true;
refineCubes[leafI] = 1;
}
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]
);
}
if( sThickness > VSMALL )
{
useNLayers = true;
const scalar cs = oc.size(rootBox);
const label numLayers =
ceil(sThickness / cs);
nLayers[leafI] =
Foam::max
(
nLayers[leafI],
Foam::max(numLayers, 1)
);
rThickness[leafI] = max(rThickness[leafI], sThickness);
targetRefLevel[leafI] =
Foam::max
(
targetRefLevel[leafI],
surfLevel
);
}
}
}
......@@ -600,19 +599,25 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
//- mark additional boxes for refinement to achieve
//- correct refinement distance
reduce(useNLayers, maxOp<label>());
reduce(changed, maxOp<bool>());
if( useNLayers )
if( useNLayers && changed )
{
octreeModifier.markAdditionalLayers
octreeModifier.refineSelectedBoxesAndAdditionalLayers
(
refineCubes,
nLayers,
rThickness,
targetRefLevel
);
}
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
if( nMarked != 0 )
changed = true;
}
else if( changed )
{
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
}
# ifdef OCTREETiming
const scalar refTime = omp_get_wtime();
......@@ -654,8 +659,6 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
Info << "Refining boxes intersecting edge meshes" << endl;
label nMarked;
//- read edge meshes and calculate the refinement level for each
//- edge mesh
const dictionary& edgeDict = meshDictPtr_->subDict("edgeMeshRefinement");
......@@ -684,6 +687,7 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
edgeMeshesPtr.set(emI, new edgeMesh(fName));
}
label nMarked;
direction addLevel(0);
if( dict.found("cellSize") )
{
......@@ -735,16 +739,17 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
const LongList<meshOctreeCube*>& leaves = octreeModifier.leavesAccess();
DynList<label> leavesInBox, intersectedLeaves;
bool changed;
do
{
changed = false;
# ifdef OCTREETiming
const scalar startIter = omp_get_wtime();
# endif
nMarked = 0;
labelList refineCubes(leaves.size(), 0);
labelList nLayers(leaves.size(), 0);
scalarList rThickness(leaves.size(), 0.0);
List<direction> targetRefLevel(leaves.size(), direction(0));
bool useNLayers(false);
......@@ -756,8 +761,7 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
const edgeList& edges = edgeMesh.edges();
# ifdef USE_OMP
# pragma omp parallel for \
reduction( + : nMarked) schedule(dynamic, 10) \
# pragma omp parallel for schedule(dynamic, 10) \
private(leavesInBox,intersectedLeaves)
# endif
forAll(edges, edgeI)
......@@ -797,9 +801,8 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
<< " for refinement" << endl;
# endif
if( !refineCubes[leafI] )
++nMarked;
refineCubes[leafI] = 1;
changed = true;
}
}
}
......@@ -811,15 +814,9 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
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[emI]/cs), 1)
);
rThickness[leafI] =
Foam::max(rThickness[leafI], refThickness[emI]);
targetRefLevel[leafI] =
Foam::max
......@@ -835,19 +832,21 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
//- mark additional boxes for refinement to achieve
//- correct refinement distance
reduce(useNLayers, maxOp<label>());
if( useNLayers )
reduce(changed, maxOp<bool>());
if( useNLayers && changed )
{
nMarked +=
octreeModifier.markAdditionalLayers
(
refineCubes,
nLayers,
targetRefLevel
);
octreeModifier.refineSelectedBoxesAndAdditionalLayers
(
refineCubes,
rThickness,
targetRefLevel
);
}
else if( changed )
{
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
}
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
# ifdef OCTREETiming
const scalar refTime = omp_get_wtime();
......@@ -856,8 +855,7 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
if( octree_.neiProcs().size() != 0 )
{
reduce(nMarked, sumOp<label>());
if( nMarked )
if( changed )
{
octreeModifier.distributeLeavesToProcessors();
......@@ -875,7 +873,7 @@ void meshOctreeCreator::refineBoxesIntersectingEdgeMeshes()
# endif
}
}
} while( nMarked != 0 );
} while( changed );
Info << "Finished refinement of boxes intersecting edge meshes" << endl;
}
......
meshLibrary/utilities/octrees/meshOctree/meshOctreeModifier/meshOctreeModifier.H
View file @ 1a8e96b1
......@@ -130,6 +130,13 @@ public:
const label nLayers = 1
) const;
//- mark additional layers around the leaves selected for refinement
void markAdditionalLayersOfFaceNeighbours
(
labelList& refineBox,
const label nLayers = 1
) const;
//- mark additional layers around the leaves selected for refinement
//- given on a box-by-box basis
//- returns the number of boxes selected for refinement
......@@ -149,6 +156,14 @@ public:
const bool hexRefinement = false
);
//- refine selected boxes and the boxes within the given range
void refineSelectedBoxesAndAdditionalLayers
(
labelList& refineBox,
const scalarList& refThickness,
const List<direction>& targetRefLevel
);
// functions for parallel runs
//- distribute leaves of the initial octree to processors
//- each processor creates a list of neighbouring processors
......
meshLibrary/utilities/octrees/meshOctree/meshOctreeModifier/meshOctreeModifierEnsureCorrectRegularity.C
View file @ 1a8e96b1
......@@ -43,8 +43,6 @@ namespace Foam
void meshOctreeModifier::ensureCorrectRegularity(labelList& refineBox)
{
const FixedList<meshOctreeCubeCoordinates, 26>& rp =
octree_.regularityPositions_;
const LongList<meshOctreeCube*>& leaves = octree_.leaves_;
//- this is needed for parallel runs to reduce the bandwidth
......@@ -57,59 +55,38 @@ void meshOctreeModifier::ensureCorrectRegularity(labelList& refineBox)
front.append(leafI);
}
FixedList<meshOctreeCube*, 26> neighbours;
DynList<label> neighbours;
label nMarked;
do
{
nMarked = 0;
transferCoordinates.clear();
LongList<meshOctreeCubeCoordinates> processorChecks;
# ifdef USE_OMP
# pragma omp parallel if( front.size() > 1000 ) \
private(neighbours) reduction(+ : nMarked)
# pragma omp parallel private(neighbours)
# endif
{
labelLongList tFront;
# ifdef USE_OMP
# pragma omp for
# pragma omp for schedule(dynamic, 50)
# endif
forAll(front, i)
tFront.append(front[i]);
# ifdef USE_OMP
# pragma omp barrier
# endif
front.clear();
while( tFront.size() != 0 )
{
const label leafI = tFront.removeLastElement();
const label leafI = front[i];
const meshOctreeCube* oc = leaves[leafI];
forAll(rp, posI)
{
const meshOctreeCubeCoordinates cc
(
oc->coordinates() + rp[posI]
);
neighbours.clear();
octree_.findAllLeafNeighbours(*oc, neighbours);
const label neiLabel = octree_.findLeafLabelForPosition(cc);
forAll(neighbours, neiI)
{
const label nei = neighbours[neiI];
if( neiLabel > -1 )
if( nei == 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
......@@ -120,25 +97,50 @@ void meshOctreeModifier::ensureCorrectRegularity(labelList& refineBox)