Skip to content
GitLab
Commit dddd354d authored by Franjo's avatar Franjo
Browse files

Merge branch 'defect-refinementThickness' into developmentPublicRepo

parents 87b24cd2 243c27d5
Hide whitespace changes
Inline Side-by-side
meshLibrary/utilities/helperClasses/parallelHelpers/meshOctreeCubeCoordinatesScalar/meshOctreeCubeCoordinatesScalar.H 0 → 100644
View file @ dddd354d
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | cfMesh: A library for mesh generation
\\ / O peration |
\\ / A nd | Author: Franjo Juretic (franjo.juretic@c-fields.com)
\\/ M anipulation | Copyright (C) Creative Fields, Ltd.
-------------------------------------------------------------------------------
License
This file is part of cfMesh.
cfMesh 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.
cfMesh 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 cfMesh. If not, see <http://www.gnu.org/licenses/>.
Class
meshOctreCubeCoordinatesScalar
Description
A class containing meshOctreeCubeCoordinates and a scalar value.
It is used for exchanging data over processors
SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef meshOctreeCubeCoordinatesScalar_H
#define meshOctreeCubeCoordinatesScalar_H
#include "scalar.H"
#include "meshOctreeCubeCoordinates.H"
#include "contiguous.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class meshOctreeCubeCoordinatesScalar Declaration
\*---------------------------------------------------------------------------*/
class meshOctreeCubeCoordinatesScalar
{
// Private data
//- scalar
scalar sValue_;
//- cube coordinates
meshOctreeCubeCoordinates coordinates_;
public:
// Constructors
//- Null construct
meshOctreeCubeCoordinatesScalar()
:
sValue_(0.0),
coordinates_()
{}
//- Construct from label and cube coordinates
meshOctreeCubeCoordinatesScalar
(
const scalar s,
const meshOctreeCubeCoordinates& cc
)
:
sValue_(s),
coordinates_(cc)
{}
// Destructor
~meshOctreeCubeCoordinatesScalar()
{}
// Member functions
//- return cube label
inline scalar scalarValue() const
{
return sValue_;
}
//- return the value
inline const meshOctreeCubeCoordinates& coordinates() const
{
return coordinates_;
}
// Member operators
inline bool operator<(const meshOctreeCubeCoordinatesScalar& scc) const
{
if( coordinates_ < scc.coordinates_ )
return true;
return false;
}
inline void operator=(const meshOctreeCubeCoordinatesScalar& scc)
{
sValue_ = scc.sValue_;
coordinates_ = scc.coordinates_;
}
inline bool operator==
(
const meshOctreeCubeCoordinatesScalar& scc
) const
{
if( coordinates_ == scc.coordinates_ )
return true;
return false;
}
inline bool operator!=
(
const meshOctreeCubeCoordinatesScalar& scc
) const
{
return !this->operator==(scc);
}
// Friend operators
friend Ostream& operator<<
(
Ostream& os,
const meshOctreeCubeCoordinatesScalar& scc
)
{
os << token::BEGIN_LIST;
os << scc.sValue_ << token::SPACE;
os << scc.coordinates_ << token::END_LIST;
// Check state of Ostream
os.check
(
"operator<<(Ostream&, const meshOctreCubeCoordinatesScalar&"
);
return os;
}
friend Istream& operator>>
(
Istream& is,
meshOctreeCubeCoordinatesScalar& scc
)
{
// Read beginning of meshOctreCubeCoordinatesScalar
is.readBegin("meshOctreCubeCoordinatesScalar");
is >> scc.sValue_;
is >> scc.coordinates_;
// Read end of meshOctreCubeCoordinatesScalar
is.readEnd("meshOctreCubeCoordinatesScalar");
// Check state of Istream
is.check("operator>>(Istream&, meshOctreCubeCoordinatesScalar");
return is;
}
};
//- Specify data associated with meshOctreCubeCoordinatesScalar
//- type is contiguous
template<>
inline bool contiguous<meshOctreeCubeCoordinatesScalar>() {return true;}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
meshLibrary/utilities/octrees/meshOctree/meshOctree.H
View file @ dddd354d
......@@ -138,6 +138,14 @@ public:
//- find a cube containing the vertex
label findLeafContainingVertex(const point&) const;
//- find leaves within the given range from the given point
void findLeavesInSphere
(
const point&,
const scalar,
DynList<label>&
) const;
//- is octree a quadtree or an octree
inline bool isQuadtree() const;
......@@ -334,6 +342,15 @@ public:
LongList<meshOctreeCubeCoordinates>& dataToReceive
) const;
//- exchange requests with other processors generating the octree
void exchangeRequestsWithNeighbourProcessors
(
const LongList<meshOctreeCubeCoordinates>& dataToSend,
const LongList<scalar>& rangesToSend,
LongList<meshOctreeCubeCoordinates>& dataToReceive,
LongList<scalar>& receivedRanges
) const;
//- neighbour processors of the current one
inline const labelList& neiProcs() const;
};
......
meshLibrary/utilities/octrees/meshOctree/meshOctreeCreator/meshOctreeCreatorAdjustOctreeToSurface.C
View file @ dddd354d
......@@ -62,24 +62,22 @@ 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);
List<direction> targetLevel(leaves.size(), direction(0));
scalarList rThickness(leaves.size(), 0.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 +93,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,22 +107,19 @@ 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] =
Foam::max(targetLevel[leafI], surfRefLevel_[triI]);
}
}
if( refine )
{
refineCubes[leafI] = 1;
++nMarked;
changed = true;
}
}
}
......@@ -134,19 +127,20 @@ 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 )
{
nMarked +=
octreeModifier.markAdditionalLayers
(
refineCubes,
nLayers,
targetLevel
);
octreeModifier.refineSelectedBoxesAndAdditionalLayers
(
refineCubes,
rThickness
);
}
else if( changed )
{
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
}
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
# ifdef OCTREETiming
const scalar refTime = omp_get_wtime();
......@@ -155,9 +149,7 @@ void meshOctreeCreator::refineBoundary()
if( Pstream::parRun() )
{
reduce(nMarked, sumOp<label>());
if( nMarked )
if( changed )
{
octreeModifier.distributeLeavesToProcessors();
......@@ -176,7 +168,7 @@ void meshOctreeCreator::refineBoundary()
}
}
} while( nMarked );
} while( changed );
Info << "Finished refining boundary boxes" << endl;
}
......@@ -288,24 +280,23 @@ 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);
List<direction> targetRefLevel(leaves.size(), direction(0));
scalarList rThickness(leaves.size(), 0.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)
{
......@@ -339,16 +330,9 @@ 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]);
useNLayers = true;
}
}
......@@ -357,26 +341,27 @@ 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
);
}
else if( changed )
{
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
}
//- refine boxes
octreeModifier.refineSelectedBoxes(refineCubes, hexRefinement_);
# ifdef OCTREETiming
const scalar refTime = omp_get_wtime();
......@@ -385,8 +370,7 @@ void meshOctreeCreator::refineBoxesContainedInObjects()
if( octree_.neiProcs().size() != 0 )
{
reduce(nMarked, sumOp<label>());
if( nMarked )
if( changed )
{
octreeModifier.distributeLeavesToProcessors();
......@@ -405,7 +389,7 @@ void meshOctreeCreator::refineBoxesContainedInObjects()
}
}
} while( nMarked != 0 );
} while( changed );
//- delete coordinate modifier if it exists
deleteDemandDrivenData(cModPtr);
......@@ -436,7 +420,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]);
......@@ -503,36 +487,40 @@ 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));
scalarField rThickness(leaves.size(), 0.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];
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,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 +536,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,46 +544,38 @@ void meshOctreeCreator::refineBoxesIntersectingSurfaces()
if( oc.intersectsTriangleExact(surf, rootBox, triI) )
{
intersectedLeaves.append(leafI);
# ifdef DEBUGSearch
Info << "Marking leaf " << leafI
<< " with coordinates " << oc
<< " for refinement" << endl;
# endif
if( oc.level() < refLevels[surfI] )
if( oc.level() < surfLevel )
{
# ifdef DEBUGSearch
Info << "Marking leaf " << leafI
<< " with coordinates " << oc
<< " for refinement" << endl;
# endif
if( !refineCubes[leafI] )
++nMarked;
//- mark boxes for refinement
changed = true;
refineCubes[leafI] = 1;
}
}
}
if( refThickness[surfI] > VSMALL )
{
useNLayers = true;
if( sThickness > VSMALL )
{
useNLayers = true;
forAll(intersectedLeaves, i)
{
const label leafI = intersectedLeaves[i];
const meshOctreeCube& oc = *leaves[leafI];
const scalar cs = oc.size(rootBox);
const scalar cs = oc.size(rootBox);
const label numLayers =
ceil(sThickness / cs);
nLayers[leafI] =
Foam::max
(
nLayers[leafI],
max(ceil(refThickness[surfI]/cs), 1)
);
nLayers[leafI] =
Foam::max
(
nLayers[leafI],
Foam::max(numLayers, 1)
);
targetRefLevel[leafI] =
Foam::max
(
targetRefLevel[leafI],
refLevels[surfI]
);
rThickness[leafI] =
max(rThickness[leafI], sThickness);
}
}