Commit 965e43d5 authored by franjo_j@hotmail.com's avatar franjo_j@hotmail.com
Browse files

A mostly-working version of boundary layer refinement. Interrupted

layers require additional improvement.


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@36 fdcce57e-7e00-11e2-b579-49867b4cea03
parent bbbd2032
......@@ -42,6 +42,7 @@ Description
#include "meshSurfaceOptimizer.H"
#include "topologicalCleaner.H"
#include "boundaryLayers.H"
#include "refineBoundaryLayers.H"
#include "renameBoundaryPatches.H"
#include "checkMeshDict.H"
#include "checkCellConnectionsOverFaces.H"
......@@ -241,12 +242,13 @@ void cartesianMeshGenerator::optimiseMeshSurface()
# endif
}
void cartesianMeshGenerator::generateBoudaryLayers()
void cartesianMeshGenerator::generateBoundaryLayers()
{
boundaryLayers bl(mesh_);
if( meshDict_.found("boundaryLayers") )
/* if( meshDict_.found("boundaryLayers") )
{
if( )
wordList createLayers;
if( meshDict_.isDict("boundaryLayers") )
......@@ -268,6 +270,9 @@ void cartesianMeshGenerator::generateBoudaryLayers()
//bl.createOTopologyLayers();
bl.addLayerForAllPatches();
}
*/
bl.addLayerForAllPatches();
# ifdef DEBUG
# ifdef DEBUGfpma
......@@ -280,6 +285,15 @@ void cartesianMeshGenerator::generateBoudaryLayers()
# endif
}
void cartesianMeshGenerator::refBoundaryLayers()
{
refineBoundaryLayers refLayers(mesh_);
refineBoundaryLayers::readSettings(meshDict_, refLayers);
refLayers.refineLayers();
}
void cartesianMeshGenerator::optimiseFinalMesh()
{
//- final optimisation
......@@ -338,10 +352,12 @@ void cartesianMeshGenerator::generateMesh()
optimiseMeshSurface();
generateBoudaryLayers();
generateBoundaryLayers();
optimiseFinalMesh();
refBoundaryLayers();
renumberMesh();
replaceBoundaries();
......
......@@ -89,7 +89,10 @@ class cartesianMeshGenerator
void optimiseMeshSurface();
//- add boundary layers
void generateBoudaryLayers();
void generateBoundaryLayers();
//- refine boundary layers
void refBoundaryLayers();
//- mesh optimisation
void optimiseFinalMesh();
......
......@@ -76,6 +76,12 @@ refineBoundaryLayers::~refineBoundaryLayers()
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void refineBoundaryLayers::avoidRefinement()
{
globalNumLayers_ = 1;
numLayersForPatch_.clear();
}
void refineBoundaryLayers::setGlobalNumberOfLayers(const label nLayers)
{
if( nLayers < 2 )
......@@ -197,6 +203,19 @@ void refineBoundaryLayers::setInteruptForPatch(const word& patchName)
void refineBoundaryLayers::refineLayers()
{
bool refinePatch(false);
for
(
std::map<word, label>::const_iterator it=numLayersForPatch_.begin();
it!=numLayersForPatch_.end();
++it
)
if( it->second > 1 )
refinePatch = true;
if( !(globalNumLayers_ > 1 || refinePatch) )
return;
Info << "Starting refining boundary layers" << endl;
if( done_ )
......@@ -232,6 +251,99 @@ void refineBoundaryLayers::refineLayers()
Info << "Finished refining boundary layers" << endl;
}
void refineBoundaryLayers::readSettings
(
const dictionary& meshDict,
refineBoundaryLayers& refLayers
)
{
if( meshDict.isDict("boundaryLayers") )
{
const dictionary& bndLayers = meshDict.subDict("boundaryLayers");
//- read global properties
if( bndLayers.found("nLayers") )
{
const label nLayers = readLabel(bndLayers.lookup("nLayers"));
refLayers.setGlobalNumberOfLayers(nLayers);
}
if( bndLayers.found("thicknessRatio") )
{
const scalar ratio = readScalar(bndLayers.lookup("thicknessRatio"));
refLayers.setGlobalThicknessRatio(ratio);
}
if( bndLayers.found("maxFirstLayerThickness") )
{
const scalar maxFirstThickness =
readScalar(bndLayers.lookup("maxFirstLayerThickness"));
refLayers.setGlobalMaxThicknessOfFirstLayer(maxFirstThickness);
}
//- patch-based properties
if( bndLayers.isDict("patchBoundaryLayers") )
{
const dictionary& patchBndLayers =
bndLayers.subDict("patchBoundaryLayers");
const wordList patchNames = patchBndLayers.toc();
forAll(patchNames, patchI)
{
const word pName = patchNames[patchI];
if( patchBndLayers.isDict(pName) )
{
const dictionary& patchDict =
patchBndLayers.subDict(pName);
if( patchDict.found("nLayers") )
{
const label nLayers =
readLabel(patchDict.lookup("nLayers"));
refLayers.setNumberOfLayersForPatch(pName, nLayers);
}
if( patchDict.found("thicknessRatio") )
{
const scalar ratio =
readScalar(patchDict.lookup("thicknessRatio"));
refLayers.setThicknessRatioForPatch(pName, ratio);
}
if( patchDict.found("maxFirstLayerThickness") )
{
const scalar maxFirstThickness =
readScalar
(
patchDict.lookup("maxFirstLayerThickness")
);
refLayers.setMaxThicknessOfFirstLayerForPatch
(
pName,
maxFirstThickness
);
}
if( patchDict.found("allowDiscontinuity") )
{
const bool allowDiscontinuity =
readBool(patchDict.lookup("allowDiscontinuity"));
if( allowDiscontinuity )
refLayers.setInteruptForPatch(pName);
}
}
else
{
Warning << "Cannot refine layer for patch "
<< patchNames[patchI] << endl;
}
}
}
}
else
{
//- the layer will not be refined
refLayers.avoidRefinement();
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
......
......@@ -149,7 +149,7 @@ class refineBoundaryLayers
(
const face& f,
const FixedList<label, 2>& nLayersInDirection,
DynList<DynList<label, 4> >& newFaces
DynList<DynList<label, 4>, 128>& newFaces
);
//- map split edges onto a cell
......@@ -196,6 +196,9 @@ public:
~refineBoundaryLayers();
// Public member functions
//- set no refinement flag
void avoidRefinement();
//- set the global number of boundary layers
void setGlobalNumberOfLayers(const label nLayers);
......@@ -237,8 +240,10 @@ public:
//- performs refinement based on the given settings
void refineLayers();
};
//- read the settings from dictionary
static void readSettings(const dictionary&, refineBoundaryLayers&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -32,7 +32,7 @@ Description
#include "FixedList.H"
#include "helperFunctions.H"
#define DEBUGLayer
//#define DEBUGLayer
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -45,7 +45,7 @@ void refineBoundaryLayers::refineFace
(
const face& f,
const FixedList<label, 2>& nLayersInDirection,
DynList<DynList<label, 4> >& newFaces
DynList<DynList<label, 4>, 128>& newFaces
)
{
//- this face must be a quad
......@@ -415,8 +415,7 @@ void refineBoundaryLayers::generateNewFaces()
}
//- refine the face
DynList<DynList<label, 4> > newFacesForFace;
DynList<DynList<label, 4>, 128> newFacesForFace;
refineFace(f, nRefinementInDirection, newFacesForFace);
//- store decomposed faces
......@@ -426,11 +425,13 @@ void refineBoundaryLayers::generateNewFaces()
newFaces_.appendList(newFacesForFace[fI]);
}
# ifdef DEBUGLayer
Info << "Internal face " << faceI << " with points " << f
<< " is refined " << endl;
forAllRow(facesFromFace_, faceI, i)
Info << "New face " << i << " is "
<< newFaces_[facesFromFace_(faceI, i)] << endl;
# endif
}
//- refine boundary faces where needed
......@@ -477,8 +478,7 @@ void refineBoundaryLayers::generateNewFaces()
}
//- refine the face
DynList<DynList<label, 4> > newFacesForFace;
DynList<DynList<label, 4>, 128> newFacesForFace;
refineFace(bf, nRefinementInDirection, newFacesForFace);
//- store the refined faces
......@@ -539,7 +539,18 @@ void refineBoundaryLayers::generateNewFaces()
nLayersAtBndFace_[beFaces(beI, 0)];
//- add the data to the list for sending
const label dir = ((eI + 1) % 2);
const label dir = (eI % 2);
# ifdef DEBUGLayer
Pout << "Face " << fI << " owner of proc patch "
<< procBoundaries[patchI].myProcNo()
<< " nei proc "
<< procBoundaries[patchI].neiProcNo()
<< " bnd face patch "
<< facePatches[beFaces(beI, 0)]
<< " direction " << dir
<< " nSplits " << nSplits0 << endl;
# endif
//- add face label, direction
//- and the number of splits
......@@ -582,7 +593,7 @@ void refineBoundaryLayers::generateNewFaces()
while( counter < receivedData.size() )
{
const label fI = receivedData[counter++];
const label dir = receivedData[counter++];
const label dir = ((receivedData[counter++] + 1) % 2);
const label nSplits = receivedData[counter++];
DynList<labelPair, 2>& currentSplits = localSplits[start+fI];
......@@ -622,7 +633,6 @@ void refineBoundaryLayers::generateNewFaces()
}
//- split the face and add the faces to the list
DynList<DynList<label, 4> > facesFromFace;
if( procBoundaries[patchI].owner() )
{
//- this processor owns this patch
......@@ -632,6 +642,16 @@ void refineBoundaryLayers::generateNewFaces()
nLayersInDirection[dirSplits[i].first()] =
dirSplits[i].second();
# ifdef DEBUGLayer
Pout << "Face " << fI << " at owner processor "
<< procBoundaries[patchI].myProcNo()
<< " neighbour processor "
<< procBoundaries[patchI].neiProcNo()
<< " face " << faces[faceI] << " refinement direction "
<< nLayersInDirection << endl;
# endif
DynList<DynList<label, 4>, 128> facesFromFace;
refineFace(faces[faceI], nLayersInDirection, facesFromFace);
//- add faces
......@@ -651,7 +671,17 @@ void refineBoundaryLayers::generateNewFaces()
dirSplits[i].second();
const face rFace = faces[faceI].reverseFace();
Pout << "Refining face " << rFace << endl;
# ifdef DEBUGLayer
Pout << "Face " << fI << " at owner processor "
<< procBoundaries[patchI].myProcNo()
<< " neighbour processor "
<< procBoundaries[patchI].neiProcNo()
<< " face " << rFace << " refinement direction "
<< nLayersInDirection << endl;
# endif
DynList<DynList<label, 4>, 128> facesFromFace;
refineFace(rFace, nLayersInDirection, facesFromFace);
forAll(facesFromFace, i)
......@@ -665,12 +695,72 @@ void refineBoundaryLayers::generateNewFaces()
}
}
}
# ifdef DEBUGLayer
returnReduce(1, sumOp<label>());
for(label procI=0;procI<Pstream::nProcs();++procI)
{
if( procI == Pstream::myProcNo() )
{
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
const label size = procBoundaries[patchI].patchSize();
for(label fI=0;fI<size;++fI)
{
const label faceI = start + fI;
const face& f = faces[faceI];
Pout << "Face " << fI << " in patch "
<< procBoundaries[patchI].patchName()
<< " has nodes " << f
<< " local splits " << localSplits[faceI]
<< " new faces from face " << facesFromFace_[faceI]
<< endl;
Pout << " Face points ";
forAll(f, pI)
Pout << mesh_.points()[f[pI]] << " ";
Pout << endl;
forAllRow(facesFromFace_, faceI, ffI)
{
const label nfI = facesFromFace_(faceI, ffI);
Pout << "New face " << ffI << " with label " << nfI
<< " consists of points ";
forAllRow(newFaces_, nfI, pI)
Pout << mesh_.points()[newFaces_(nfI, pI)]
<< " ";
Pout << endl;
}
}
}
}
returnReduce(1, sumOp<label>());
}
returnReduce(1, sumOp<label>());
//::exit(1);
# endif
}
# ifdef DEBUGLayer
Info << "facesFromFace_ " << facesFromFace_ << endl;
Info << "newFaces_ " << newFaces_ << endl;
returnReduce(1, sumOp<label>());
for(label procI=0;procI<Pstream::nProcs();++procI)
{
if( procI == Pstream::myProcNo() )
{
Pout << "facesFromFace_ " << facesFromFace_ << endl;
Pout << "newFaces_ " << newFaces_ << endl;
}
returnReduce(1, sumOp<label>());
}
# endif
Info << "Finished refining boundary-layer faces " << endl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -39,7 +39,7 @@ Description
#include <omp.h>
# endif
#define DEBUGLayer
//#define DEBUGLayer
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -448,7 +448,7 @@ void refineBoundaryLayers::analyseLayers()
}
# ifdef DEBUGLayer
Info << "nLayersAtPatch " << nLayersAtPatch << endl;
Pout << "nLayersAtPatch " << nLayersAtPatch << endl;
# endif
//- set the number of boundary layers which shall be generated above
......@@ -474,7 +474,9 @@ void refineBoundaryLayers::analyseLayers()
}
# ifdef DEBUGLayer
Info << "nLayersAtBndFace_ " << nLayersAtBndFace_ << endl;
forAll(nLayersAtBndFace_, bfI)
Pout << "Boundary face " << bfI << " in patch "
<< facePatch[bfI] << " num layers " << nLayersAtBndFace_[bfI] << endl;
//::exit(1);
# endif
}
......@@ -638,6 +640,7 @@ bool refineBoundaryLayers::findSplitEdges()
const meshSurfaceEngine& mse = surfaceEngine();
const faceList::subList& bFaces = mse.boundaryFaces();
const labelList& facePatch = mse.boundaryFacePatches();
mse.faceOwners();
const VRWGraph& pFaces = mse.pointFaces();
const labelList& bp = mse.bp();
......@@ -719,7 +722,15 @@ bool refineBoundaryLayers::findSplitEdges()
reduce(validLayer, minOp<bool>());
# ifdef DEBUGLayer
Info << "Generated split edges " << splitEdges_ << endl;
for(label procI=0;procI<Pstream::nProcs();++procI)
{
if( procI == Pstream::myProcNo() )
{
Pout << "Generated split edges " << splitEdges_ << endl;
}
returnReduce(1, sumOp<label>());
}
# endif
return validLayer;
......@@ -912,7 +923,18 @@ void refineBoundaryLayers::generateNewVertices()
{
const labelPair& lp = receivedNumLayers[i];
const label eI = globalToLocal[lp.first()];
nNodesAtEdge[eI] = std::max(nNodesAtEdge[eI], lp.second());
const edge& e = edges[eI];
label seI(-1);
forAllRow(splitEdgesAtPoint_, e.start(), i)
{
const label seJ = splitEdgesAtPoint_(e.start(), i);
if( splitEdges_[seJ] == e )
{
seI = seJ;
break;
}
}
nNodesAtEdge[seI] = std::max(nNodesAtEdge[seI], lp.second());
}
//- exchange thickness ratio
......@@ -923,10 +945,21 @@ void refineBoundaryLayers::generateNewVertices()
{
const labelledScalar& ls = receivedScalar[i];
const label eI = globalToLocal[ls.scalarLabel()];
thicknessRatio[eI] = std::max(thicknessRatio[eI], ls.value());
const edge& e = edges[eI];
label seI(-1);
forAllRow(splitEdgesAtPoint_, e.start(), i)
{
const label seJ = splitEdgesAtPoint_(e.start(), i);
if( splitEdges_[seJ] == e )
{
seI = seJ;
break;
}
}
thicknessRatio[seI] = std::max(thicknessRatio[seI], ls.value());
}
//- exchange ,aximum thickness of the first layer
//- exchange maximum thickness of the first layer
receivedScalar.clear();
help::exchangeMap(exchangeThickness, receivedScalar);
......@@ -934,8 +967,19 @@ void refineBoundaryLayers::generateNewVertices()
{
const labelledScalar& ls = receivedScalar[i];
const label eI = globalToLocal[ls.scalarLabel()];
firstLayerThickness[eI] =
std::min(firstLayerThickness[eI], ls.value());
const edge& e = edges[eI];
label seI(-1);
forAllRow(splitEdgesAtPoint_, e.start(), i)
{
const label seJ = splitEdgesAtPoint_(e.start(), i);
if( splitEdges_[seJ] == e )
{
seI = seJ;
break;
}
}
firstLayerThickness[seI] =
std::min(firstLayerThickness[seI], ls.value());
}
}
......@@ -1033,6 +1077,27 @@ void refineBoundaryLayers::generateNewVertices()
}
# ifdef DEBUGLayer
for(label procI=0;procI<Pstream::nProcs();++procI)
{
if( procI == Pstream::myProcNo() )
{
forAll(splitEdges_, seI)
{
Pout << "\nSplit edge " << seI << " nodes " << splitEdges_[seI]
<< " coordinates " << points[splitEdges_[seI][0]]
<< " " << points[splitEdges_[seI][1]]
<< " has new points "
<< newVerticesForSplitEdge_[seI] << endl;
forAllRow(newVerticesForSplitEdge_, seI, i)
Pout << "Point " << i << " on edge ha coordinates "
<< points[newVerticesForSplitEdge_(seI, i)] << endl;
}
}
returnReduce(1, sumOp<label>());
}
Info << "Finished generating vertices at split edges" << endl;
//::exit(1);
# endif
......
......@@ -170,12 +170,61 @@ void checkMeshDict::checkBoundaryLayers() const
{
if( meshDict_.found("boundaryLayers") )
{
if( meshDict_.isDict("boundaryLayers") )
const dictionary& bndLayers = meshDict_.subDict("boundaryLayers");
//- read global properties
if( bndLayers.found("nLayers") )
{
readLabel(bndLayers.lookup("nLayers"));
}
if( bndLayers.found("thicknessRatio") )