Commit 264432b6 authored by Franjo's avatar Franjo
Browse files

Merge branch 'development' into task-qualityControls

parents a6960c7f 87b7c279
#if defined(__GNUC__)
# if defined(darwinIntel64)
OMP_FLAGS =
# else
ifeq (Gcc,$(findstring Gcc,$(WM_COMPILER)))
OMP_FLAGS = -DUSE_OMP -fopenmp
# endif
#else
else
OMP_FLAGS =
#endif
endif
EXE_INC = \
$(OMP_FLAGS) \
......
......@@ -188,7 +188,7 @@ void cartesian2DMeshGenerator::refBoundaryLayers()
void cartesian2DMeshGenerator::replaceBoundaries()
{
renameBoundaryPatches rbp(mesh_, meshDict_);
renameBoundaryPatches rbp(mesh_, meshDict_, true);
}
void cartesian2DMeshGenerator::renumberMesh()
......@@ -198,8 +198,6 @@ void cartesian2DMeshGenerator::renumberMesh()
void cartesian2DMeshGenerator::generateMesh()
{
try
{
if( controller_.runCurrentStep("templateGeneration") )
{
createCartesianMesh();
......@@ -247,18 +245,6 @@ void cartesian2DMeshGenerator::generateMesh()
replaceBoundaries();
controller_.workflowCompleted();
}
catch(const std::string& message)
{
Info << message << endl;
}
catch(...)
{
WarningIn
(
"void cartesian2DMeshGenerator::generateMesh()"
) << "Meshing process terminated!" << endl;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -283,6 +269,8 @@ cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time& time)
mesh_(time),
controller_(mesh_)
{
try
{
if( true )
{
checkMeshDict cmd(meshDict_);
......@@ -317,7 +305,10 @@ cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time& time)
if( surfacePtr_->featureEdges().size() != 0 )
{
//- get rid of duplicate triangles as they cause strange problems
triSurfaceCleanupDuplicateTriangles(const_cast<triSurf&>(*surfacePtr_));
triSurfaceCleanupDuplicateTriangles
(
const_cast<triSurf&>(*surfacePtr_)
);
//- create surface patches based on the feature edges
//- and update the meshDict based on the given data
......@@ -347,6 +338,18 @@ cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time& time)
meshOctreeCreator(*octreePtr_, meshDict_).createOctreeBoxes();
generateMesh();
}
catch(const std::string& message)
{
Info << message << endl;
}
catch(...)
{
WarningIn
(
"cartesian2DMeshGenerator::cartesian2DMeshGenerator(const Time&)"
) << "Meshing process terminated!" << endl;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
......
......@@ -247,8 +247,6 @@ void cartesianMeshGenerator::renumberMesh()
void cartesianMeshGenerator::generateMesh()
{
try
{
if( controller_.runCurrentStep("templateGeneration") )
{
createCartesianMesh();
......@@ -298,18 +296,6 @@ void cartesianMeshGenerator::generateMesh()
replaceBoundaries();
controller_.workflowCompleted();
}
catch(const std::string& message)
{
Info << message << endl;
}
catch(...)
{
WarningIn
(
"void cartesianMeshGenerator::generateMesh()"
) << "Meshing process terminated!" << endl;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -334,6 +320,8 @@ cartesianMeshGenerator::cartesianMeshGenerator(const Time& time)
mesh_(time),
controller_(mesh_)
{
try
{
if( true )
{
checkMeshDict cmd(meshDict_);
......@@ -385,6 +373,19 @@ cartesianMeshGenerator::cartesianMeshGenerator(const Time& time)
meshOctreeCreator(*octreePtr_, meshDict_).createOctreeBoxes();
generateMesh();
}
catch(const std::string& message)
{
Info << "Here" << endl;
Info << message << endl;
}
catch(...)
{
WarningIn
(
"cartesianMeshGenerator::cartesianMeshGenerator(const Time&)"
) << "Meshing process terminated!" << endl;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
......
......@@ -235,8 +235,6 @@ void tetMeshGenerator::renumberMesh()
void tetMeshGenerator::generateMesh()
{
try
{
if( controller_.runCurrentStep("templateGeneration") )
{
createTetMesh();
......@@ -286,18 +284,6 @@ void tetMeshGenerator::generateMesh()
replaceBoundaries();
controller_.workflowCompleted();
}
catch(const std::string& message)
{
Info << message << endl;
}
catch(...)
{
WarningIn
(
"void tetMeshGenerator::generateMesh()"
) << "Meshing process terminated!" << endl;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -323,6 +309,8 @@ tetMeshGenerator::tetMeshGenerator(const Time& time)
mesh_(time),
controller_(mesh_)
{
try
{
if( true )
{
checkMeshDict cmd(meshDict_);
......@@ -372,6 +360,18 @@ tetMeshGenerator::tetMeshGenerator(const Time& time)
meshOctreeCreator(*octreePtr_, meshDict_).createOctreeBoxes();
generateMesh();
}
catch(const std::string& message)
{
Info << message << endl;
}
catch(...)
{
WarningIn
(
"tetMeshGenerator::tetMeshGenerator(const Time&)"
) << "Meshing process terminated!" << endl;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
......
......@@ -523,16 +523,23 @@ inline Istream& operator>>
meshOctreeCubeCoordinates& cc
)
{
// Read beginning of meshOctreeCubeCoordinates
//- Read beginning of meshOctreeCubeCoordinates
is.readBegin("meshOctreeCubeCoordinates");
//- read the level in the octree structure
label l;
is >> l;
cc.level_ = l;
//- read x, y and z coordinates
is.readBegin("meshOctreeCubeCoordinates");
is >> cc.posX_;
is >> cc.posY_;
is >> cc.posZ_;
is.readEnd("meshOctreeCubeCoordinates");
// Read end of meshOctreeCubeCoordinates
is.readEnd("meshOctreeCubeCoordinates");
......@@ -551,10 +558,15 @@ inline Ostream& operator<<
os << token::BEGIN_LIST;
os << label(cc.level_) << token::SPACE;
os << token::BEGIN_LIST;
os << cc.posX_ << token::SPACE;
os << cc.posY_ << token::SPACE;
os << cc.posZ_ << token::END_LIST;
os << token::END_LIST;
// Check state of Ostream
os.check("operator<<(Ostream&, const meshOctreeCubeCoordinates");
......
......@@ -26,7 +26,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "meshOctreeModifier.H"
#include "helperFunctionsPar.H"
#include "helperFunctions.H"
#include "triSurf.H"
#include <map>
......@@ -203,8 +203,6 @@ void meshOctreeModifier::addLayerFromNeighbouringProcessors()
return;
const LongList<meshOctreeCube*>& leaves = octree_.leaves_;
const labelList& neiProcs = octree_.neiProcs_;
const List<Pair<meshOctreeCubeCoordinates> >& neiRange = octree_.neiRange_;
forAll(leaves, leafI)
if( leaves[leafI]->procNo() != Pstream::myProcNo() )
......@@ -212,6 +210,9 @@ void meshOctreeModifier::addLayerFromNeighbouringProcessors()
Info << "Adding an additional layer of cells" << endl;
const labelList& neiProcs = octree_.neiProcs_;
const List<Pair<meshOctreeCubeCoordinates> >& neiRange = octree_.neiRange_;
meshOctreeCubeCoordinates minCoord, maxCoord;
std::map<label, LongList<meshOctreeCubeBasic> > toProcs;
forAll(neiProcs, i)
......@@ -227,19 +228,25 @@ void meshOctreeModifier::addLayerFromNeighbouringProcessors()
forAll(neiProcs, procI)
{
if(
if
(
(maxCoord >= neiRange[procI].first()) &&
(minCoord <= neiRange[procI].second())
)
{
toProcs[neiProcs[procI]].append(*leaves[leafI]);
}
}
}
# ifdef OCTREE_DEBUG
for(label i=0;i<Pstream::nProcs();++i)
{
if( i == Pstream::myProcNo() )
{
Pout << "Neighbour processors " << neiProcs << endl;
Pout << "Neighbour range " << neiRange << endl;
std::map<label, LongList<meshOctreeCubeBasic> >::iterator it;
for(it=toProcs.begin();it!=toProcs.end();++it)
{
......
......@@ -489,8 +489,12 @@ void meshOctreeModifier::refineSelectedBoxes
}
}
//- update the list of leaves
createListOfLeaves();
//- update the communication pattern
updateCommunicationPattern();
# ifdef OCTREETiming
Info << "Time for actual refinement " << (omp_get_wtime()-regTime) << endl;
# endif
......
......@@ -30,7 +30,6 @@ Description
#include "polyMeshGenAddressing.H"
#include "helperFunctions.H"
#include "polyMeshGenChecks.H"
#include "meshOptimizer.H"
// #define DEBUGSearch
......@@ -96,13 +95,13 @@ void symmetryPlaneOptimisation::detectSymmetryPlanes()
const point c = it->second.first / it->second.second;
const std::pair<vector, label>& ns = normalSum[it->first];
const point n = ns.first / ns.second;
const vector n = ns.first / ns.second;
symmetryPlanes_.insert(std::make_pair(it->first, plane(c, n)));
}
}
void symmetryPlaneOptimisation::pointInPlanes(VRWGraph& pointInPlanes) const
bool symmetryPlaneOptimisation::pointInPlanes(VRWGraph& pointInPlanes) const
{
const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
const pointFieldPMG& points = mesh_.points();
......@@ -111,6 +110,8 @@ void symmetryPlaneOptimisation::pointInPlanes(VRWGraph& pointInPlanes) const
pointInPlanes.clear();
pointInPlanes.setSize(points.size());
bool foundProblematic(false);
forAll(boundaries, patchI)
{
if( boundaries[patchI].patchType() == "symmetryPlane" )
......@@ -121,11 +122,28 @@ void symmetryPlaneOptimisation::pointInPlanes(VRWGraph& pointInPlanes) const
{
const face& f = faces[faceI];
const point c = f.centre(points);
scalar maxDist(0.0);
forAll(f, pI)
maxDist = max(maxDist, mag(points[f[pI]] - c));
forAll(f, pI)
{
std::map<label, plane>::const_iterator it =
symmetryPlanes_.find(patchI);
if( it != symmetryPlanes_.end() )
{
const scalar d = it->second.distance(points[f[pI]]);
if( d > 0.5 * maxDist )
foundProblematic = true;
}
pointInPlanes.appendIfNotIn(f[pI], patchI);
}
}
}
}
if( Pstream::parRun() )
{
......@@ -174,6 +192,10 @@ void symmetryPlaneOptimisation::pointInPlanes(VRWGraph& pointInPlanes) const
pointInPlanes.appendIfNotIn(pointI, receivedData[counter++]);
}
}
reduce(foundProblematic, maxOp<bool>());
return foundProblematic;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -198,7 +220,17 @@ void symmetryPlaneOptimisation::optimizeSymmetryPlanes()
pointFieldPMG& points = mesh_.points();
VRWGraph pointInPlane;
pointInPlanes(pointInPlane);
if( pointInPlanes(pointInPlane) )
{
Warning
<< "Detected large deviation from some symmetry planes."
<< "Please check your settings and ensure that all your patches"
<< " with type symmetryPlane are planar. If there exists a patch"
<< " that is not planar please use the symmetry type, instead."
<< " Skipping symmetryPlane optimisation..." << endl;
return;
}
forAll(pointInPlane, pointI)
{
......@@ -211,21 +243,49 @@ void symmetryPlaneOptimisation::optimizeSymmetryPlanes()
"void symmetryPlaneOptimisation::optimizeSymmetryPlanes()"
) << "Point " << pointI << " is in more than three symmetry"
<< " planes. Cannot move it" << endl;
continue;
}
point& p = points[pointI];
vector disp(vector::zero);
for(label plI=0;plI<nPlanes;++plI)
if( nPlanes == 1 )
{
//- point is in a plane
std::map<label, plane>::const_iterator it =
symmetryPlanes_.find(pointInPlane(pointI, 0));
const point newP = it->second.nearestPoint(points[pointI]);
disp += newP - p;
p = it->second.nearestPoint(p);
}
else if( nPlanes == 2 )
{
//- point is at the edge between two planes
const plane& pl0 =
symmetryPlanes_.find(pointInPlane(pointI, 0))->second;
const plane& pl1 =
symmetryPlanes_.find(pointInPlane(pointI, 1))->second;
p += disp;
const plane::ray il = pl0.planeIntersect(pl1);
vector n = il.dir();
n /= (mag(n) + VSMALL);
const scalar lProj = (p - il.refPoint()) & n;
p = il.refPoint() + lProj * n;
}
else if( nPlanes == 3 )
{
//- points is a corner between three planes
const plane& pl0 =
symmetryPlanes_.find(pointInPlane(pointI, 0))->second;
const plane& pl1 =
symmetryPlanes_.find(pointInPlane(pointI, 1))->second;
const plane& pl2 =
symmetryPlanes_.find(pointInPlane(pointI, 2))->second;
p = pl0.planePlaneIntersect(pl1, pl2);
}
}
labelHashSet badFaces;
......
......@@ -68,7 +68,7 @@ class symmetryPlaneOptimisation
void detectSymmetryPlanes();
//- point-planes addressing
void pointInPlanes(VRWGraph&) const;
bool pointInPlanes(VRWGraph&) const;
public:
......
......@@ -26,6 +26,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "renameBoundaryPatches.H"
#include "symmetryPlaneOptimisation.H"
#include "demandDrivenData.H"
#include "IOdictionary.H"
......@@ -248,12 +249,34 @@ void renameBoundaryPatches::calculateNewBoundary()
Info << "Finished renaming boundary patches" << endl;
}
void renameBoundaryPatches::checkEmptyPatches()
{
polyMeshGenModifier meshModifier(mesh_);
forAll(mesh_.boundaries(), patchI)
{
boundaryPatch& patch = meshModifier.boundariesAccess()[patchI];
if( patch.patchType() == "empty" )
{
patch.patchType() = "wall";
}
}
}
void renameBoundaryPatches::checkSymmetryPlanes()
{
symmetryPlaneOptimisation symmSmother(mesh_);
symmSmother.optimizeSymmetryPlanes();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
renameBoundaryPatches::renameBoundaryPatches
(
polyMeshGen& mesh,
const IOdictionary& meshDict
const IOdictionary& meshDict,
const bool allowEmptyPatches
)
:
mesh_(mesh),
......@@ -261,6 +284,13 @@ renameBoundaryPatches::renameBoundaryPatches
{
if( meshDict.found("renameBoundary") )
calculateNewBoundary();
if( !allowEmptyPatches )
{
checkEmptyPatches();
}
checkSymmetryPlanes();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
......
......@@ -62,6 +62,13 @@ class renameBoundaryPatches
//- calculate new boundary and replace the existing one
void calculateNewBoundary();
//- find empty patches in 3D meshes and change them to wall
void checkEmptyPatches();
//- check if there exist any symmetry planes
//- and make sure they get perfectly flat
void checkSymmetryPlanes();
//- Disallow default bitwise copy construct
renameBoundaryPatches(const renameBoundaryPatches&);
......@@ -76,7 +83,8 @@ public:
renameBoundaryPatches
(
polyMeshGen& mesh,
const IOdictionary& meshDict
const IOdictionary& meshDict,
const bool allowEmptyPatches = false
);
// Destructor
......
......@@ -31,6 +31,7 @@ Description
#include "checkMeshDict.H"
#include <map>
#include <stdexcept>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -206,9 +207,28 @@ const triSurf* triSurfacePatchManipulator::surfaceWithPatches
}
}
bool foundProblematic(false);
forAll(patchToNewPatches, patchI)
{
const word& pName = origPatches[patchI].name();
if
(
(patchTypes[pName] == "symmetryPlane") &&
(patchToNewPatches[patchI].size() > 1)
)
{
SeriousError << "Symmetry plane patch with a name " << pName
<< " is decomposed into "
<< patchToNewPatches[patchI].size()
<< " separate parts. Please split the patch into parts "
<< "bounded by feature edges"
<< " or change the type to symmetry."
<< endl;
foundProblematic = true;
}