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

Refinement of bnd layer (mostly finished) + templating of helper

functions + bug fixes


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@35 fdcce57e-7e00-11e2-b579-49867b4cea03
parent 9211bebd
......@@ -53,7 +53,7 @@ refineBoundaryLayers::refineBoundaryLayers(polyMeshGen& mesh)
msePtr_(NULL),
globalNumLayers_(1),
globalThicknessRatio_(1.0),
globalMaxThicknessFirstLayer_(-1.0),
globalMaxThicknessFirstLayer_(VGREAT),
numLayersForPatch_(),
thicknessRatioForPatch_(),
maxThicknessForPatch_(),
......
......@@ -32,6 +32,8 @@ Description
#include "FixedList.H"
#include "helperFunctions.H"
#define DEBUGLayer
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
......@@ -123,7 +125,7 @@ void refineBoundaryLayers::refineFace
}
}
/*
# ifdef DEBUGLayer
Info << "Refining face " << f << endl;
Info << "Splits in direction " << nLayersInDirection << endl;
Info << "Here " << endl;
......@@ -131,7 +133,7 @@ void refineBoundaryLayers::refineFace
Info << "dir0Edges " << dir0Edges << endl;
Info << "Dir1 " << dir1 << endl;
Info << "dir1Edges " << dir1Edges << endl;
*/
# endif
//- in case of only one refinement direction, it must direction 0
if( (dir1 != -1) && (dir0 == -1) )
......@@ -140,7 +142,7 @@ void refineBoundaryLayers::refineFace
dir0Edges = dir1Edges;
dir1 = -1;
}
else if( (dir0 != -1) && (dir1 != -1) && (dir1 != f.rcIndex(dir0)) )
else if( (dir0 != -1) && (dir1 != -1) && (dir1 != f.fcIndex(dir0)) )
{
//- alternate value to preserve correct face orientation
const label add = dir0;
......@@ -156,27 +158,31 @@ void refineBoundaryLayers::refineFace
const label nLayersDir0 = dir0>=0?nLayersInDirection[dir0%2]:1;
const label nLayersDir1 = dir1>=0?nLayersInDirection[dir1%2]:1;
/*
# ifdef DEBUGLayer
Info << "Face has points " << f << endl;
Info << "dirEdges0 " << dir0Edges << endl;
Info << "dir1Edges " << dir1Edges << endl;
if( dir0 >= 0 )
{
Info << "Points on edge " << dir0Edges.first()
Info << "Points on edge " << dir0Edges.first() << " with nodes "
<< splitEdges_[dir0Edges.first()]
<< " are " << newVerticesForSplitEdge_[dir0Edges.first()] << endl;
Info << "Points on edge " << dir0Edges.second()
Info << "Points on edge " << dir0Edges.second() << " with nodes "
<< splitEdges_[dir0Edges.second()]
<< " are " << newVerticesForSplitEdge_[dir0Edges.second()] << endl;
}
if( dir1 >= 0 )
{
Info << "Points on edge " << dir1Edges.first()
Info << "Points on edge " << dir1Edges.first() << " with nodes "
<< splitEdges_[dir1Edges.first()]
<< " are " << newVerticesForSplitEdge_[dir1Edges.first()] << endl;
Info << "Points on edge " << dir1Edges.second()
Info << "Points on edge " << dir1Edges.second() << " with nodes "
<< splitEdges_[dir1Edges.second()]
<< " are " << newVerticesForSplitEdge_[dir1Edges.second()] << endl;
}
Info << "nLayersDir0 " << nLayersDir0 << endl;
Info << "nLayersDir1 " << nLayersDir1 << endl;
*/
# endif
//- map the face onto a matrix for easier orientation
DynList<DynList<label> > facePoints;
......@@ -213,33 +219,51 @@ void refineBoundaryLayers::refineFace
{
if( facePoints[i][j] < 0 )
{
# ifdef DEBUGLayer
Info << "Determining u " << facePoints[0][0] << endl;
Info << "Other point " << facePoints[i][0] << endl;
Info << "Points at aplit edge "
<< newVerticesForSplitEdge_[dir0Edges.second()] << endl;
# endif
const scalar u
(
mag(points[facePoints[i][0]] - points[facePoints[0][0]]) /
splitEdges_[dir0Edges.second()].mag(points)
);
# ifdef DEBUGLayer
Info << "Determining v " << facePoints[0][0] << endl;
Info << "Other point " << facePoints[0][j] << endl;
Info << "Points at aplit edge "
<< newVerticesForSplitEdge_[dir1Edges.first()] << endl;
# endif
const scalar v
(
mag(points[facePoints[0][i]] - points[facePoints[0][0]]) /
mag(points[facePoints[0][j]] - points[facePoints[0][0]]) /
splitEdges_[dir1Edges.first()].mag(points)
);
# ifdef DEBUGLayer
Info << "Generating point of face " << endl;
Info << "u = " << u << endl;
Info << "v = " << v << endl;
# endif
//- calculate the coordinates of the missing point via
//- transfinite interpolation
const point newP
(
(1.0 - v) * points[facePoints[i][0]] +
v * points[facePoints[i][nLayersDir1]] +
(1.0 - u) * points[facePoints[0][j]] +
u * points[facePoints[nLayersDir0][j]] -
(1.0 - u) * (1.0 - v) * points[facePoints[0][0]] -
u * v * points[facePoints[nLayersDir0][0]] -
u * (1.0 - v) * points[facePoints[nLayersDir0][nLayersDir1]] -
(1.0 - u) * (1.0 - v) * points[facePoints[0][0]] +
u * (1.0 - v) * points[facePoints[nLayersDir0][0]] +
u * v * points[facePoints[nLayersDir0][nLayersDir1]] +
(1.0 - u) * v * points[facePoints[0][nLayersDir1]]
);
# ifdef DEBUGLayer
Info << "Point coordinate " << newP << endl;
# endif
//- add the vertex to the mesh
facePoints[i][j] = points.size();
......@@ -248,12 +272,14 @@ void refineBoundaryLayers::refineFace
}
}
//Info << "Face points after creating vertices " << facePoints << endl;
# ifdef DEBUGLayer
Info << "Face points after creating vertices " << facePoints << endl;
# endif
//- Finally, create the faces
for(label i=0;i<nLayersDir0;++i)
for(label j=0;j<nLayersDir1;++j)
{
for(label j=0;j<nLayersDir1;++j)
for(label i=0;i<nLayersDir0;++i)
{
if( !((i == (nLayersDir0 - 1)) && (j == (nLayersDir1 - 1))) )
{
......@@ -274,7 +300,7 @@ void refineBoundaryLayers::refineFace
//- create the last face which may not be a quad
DynList<label, 4> newF;
if( dir0 != -1 && dir1 == -1 )
if( (dir0 != -1) && (dir1 == -1) )
{
//- face is split in one direction, only
label eLabel = dir0Edges.second();
......@@ -312,10 +338,12 @@ void refineBoundaryLayers::refineFace
newFaces.append(newF);
//Info << "Input face " << f << endl;
//Info << "Decomposed faces are " << newFaces << endl;
# ifdef DEBUGLayer
Info << "Input face " << f << endl;
Info << "Decomposed faces are " << newFaces << endl;
//if( (nLayersInDirection[0] > 1) && (nLayersInDirection[1] > 1) )
//::exit(1);
# endif
}
void refineBoundaryLayers::generateNewFaces()
......@@ -397,6 +425,12 @@ void refineBoundaryLayers::generateNewFaces()
facesFromFace_.append(faceI, newFaces_.size());
newFaces_.appendList(newFacesForFace[fI]);
}
Info << "Internal face " << faceI << " with points " << f
<< " is refined " << endl;
forAllRow(facesFromFace_, faceI, i)
Info << "New face " << i << " is "
<< newFaces_[facesFromFace_(faceI, i)] << endl;
}
//- refine boundary faces where needed
......@@ -453,6 +487,14 @@ void refineBoundaryLayers::generateNewFaces()
facesFromFace_.append(faceI, newFaces_.size());
newFaces_.appendList(newFacesForFace[fI]);
}
# ifdef DEBUGLayer
Info << "Boundary face " << faceI << " with points " << bf
<< " owner cell " << mesh_.owner()[faceI] << " is refined " << endl;
forAllRow(facesFromFace_, faceI, i)
Info << "New face " << i << " is "
<< newFaces_[facesFromFace_(faceI, i)] << endl;
# endif
}
if( Pstream::parRun() )
......@@ -497,7 +539,7 @@ void refineBoundaryLayers::generateNewFaces()
nLayersAtBndFace_[beFaces(beI, 0)];
//- add the data to the list for sending
const label dir = !(eI % 2);
const label dir = ((eI + 1) % 2);
//- add face label, direction
//- and the number of splits
......@@ -553,6 +595,11 @@ void refineBoundaryLayers::generateNewFaces()
}
}
# ifdef DEBUGLayer
returnReduce(1, sumOp<label>());
Info << "Starting splitting processor boundaries" << endl;
# endif
//- perform splitting
forAll(procBoundaries, patchI)
{
......@@ -578,7 +625,7 @@ void refineBoundaryLayers::generateNewFaces()
DynList<DynList<label, 4> > facesFromFace;
if( procBoundaries[patchI].owner() )
{
//- this processor own this patch
//- this processor owns this patch
FixedList<label, 2> nLayersInDirection;
const DynList<labelPair, 2>& dirSplits = it->second;
forAll(dirSplits, i)
......@@ -604,6 +651,7 @@ void refineBoundaryLayers::generateNewFaces()
dirSplits[i].second();
const face rFace = faces[faceI].reverseFace();
Pout << "Refining face " << rFace << endl;
refineFace(rFace, nLayersInDirection, facesFromFace);
forAll(facesFromFace, i)
......@@ -619,8 +667,10 @@ void refineBoundaryLayers::generateNewFaces()
}
}
# ifdef DEBUGLayer
Info << "facesFromFace_ " << facesFromFace_ << endl;
Info << "newFaces_ " << newFaces_ << endl;
# endif
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -39,6 +39,8 @@ Description
#include <omp.h>
# endif
#define DEBUGLayer
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
......@@ -46,79 +48,76 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace bndLayerOps
{
class meshBndLayerNeighbourOperator
{
const polyMeshGen& mesh_;
const meshSurfaceEngine& mse_;
const label size_;
public:
meshBndLayerNeighbourOperator
(
const polyMeshGen& mesh,
const meshSurfaceEngine& mse
)
meshBndLayerNeighbourOperator(const meshSurfaceEngine& mse)
:
mesh_(mesh),
mse_(mse)
mse_(mse),
size_(mse.boundaryFaces().size())
{}
label size() const
{
return mesh_.cells().size();
return size_;
}
void operator()(const label cellI, DynList<label>& neighbourCells) const
void operator()(const label bfI, DynList<label>& neighbourFaces) const
{
neighbourCells.clear();
neighbourFaces.clear();
const edgeList& edges = mse_.edges();
const VRWGraph& bpEdges = mse_.boundaryPointEdges();
const labelList& bp = mse_.bp();
const cellListPMG& cells = mse_.cells();
const labelList& owner = mesh_.owner();
const labelList& neighbour = mesh_.neighbour();
const labelList& faceOwner = mse_.faceOwners();
const label cellI = faceOwner[bfI];
const cell& c = cells[cellI];
const faceListPMG& faces = mesh_.faces();
const cell& c = mesh_.cells()[cellI];
const VRWGraph& faceEdges = mse_.faceEdges();
const VRWGraph& edgeFaces = mse_.edgeFaces();
forAll(c, fI)
forAllRow(faceEdges, bfI, feI)
{
//- find the neighbour cell over this face
label nei = owner[c[fI]];
if( nei == cellI )
nei = neighbour[c[fI]];
if( nei < 0 )
continue;
//- face must be a quad
const face& f = faces[c[fI]];
if( f.size() != 4 )
continue;
const label edgeI = faceEdges(bfI, feI);
forAll(f, eI)
if( edgeFaces.sizeOfRow(edgeI) == 2 )
{
const edge e = f.faceEdge(eI);
label nei = edgeFaces(edgeI, 0);
const label bpI = bp[e.start()];
const label bpJ = bp[e.end()];
if( nei == bfI )
nei = edgeFaces(edgeI, 1);
if( (bpI < 0) || (bpJ < 0) )
//- faces must not be part of the same cell
if( faceOwner[nei] == cellI )
continue;
forAllRow(bpEdges, bpI, peI)
//- owner cell of the other face must
//- have cellI as its neighbour
const cell& neiC = cells[faceOwner[nei]];
bool sharedFace(false);
forAll(c, fI)
{
const label beI = bpEdges(bpI, peI);
if( edges[beI] == e )
forAll(neiC, fJ)
{
neighbourCells.append(nei);
break;
if( c[fI] == neiC[fJ] )
{
sharedFace = true;
break;
}
}
if( sharedFace )
break;
}
if( sharedFace )
neighbourFaces.append(nei);
}
}
}
......@@ -131,91 +130,130 @@ public:
const DynList<label>& localGroupLabel
) const
{
const PtrList<writeProcessorPatch>& procBoundaries =
mesh_.procBoundaries();
const labelList& owner = mesh_.owner();
const polyMeshGen& mesh = mse_.mesh();
const faceListPMG& faces = mesh.faces();
const cellListPMG& cells = mesh.cells();
const edgeList& edges = mse_.edges();
const labelList& faceOwner = mse_.faceOwners();
const VRWGraph& edgeFaces = mse_.edgeFaces();
const Map<label>& otherProc = mse_.otherEdgeFaceAtProc();
const Map<label>& globalToLocal = mse_.globalToLocalBndEdgeAddressing();
//- send the data to other processors
forAll(procBoundaries, patchI)
std::map<label, LongList<labelPair> > exchangeData;
forAll(mse_.beNeiProcs(), procI)
exchangeData[mse_.beNeiProcs()[procI]].clear();
forAllConstIter(Map<label>, globalToLocal, it)
{
const label start = procBoundaries[patchI].patchStart();
const label size = procBoundaries[patchI].patchSize();
const label beI = it();
//- combine data if the cell attached to this face has a face
//- attached to the inter-processor boundary
//- this must hold for boundary layer cells
const cell& c = cells[faceOwner[edgeFaces(beI, 0)]];
labelList groupOwner(procBoundaries[patchI].patchSize());
for(label faceI=0;faceI<size;++faceI)
bool validCell(false);
forAll(c, fI)
{
const label groupI = elementInGroup[owner[start+faceI]];
const face& f = faces[c[fI]];
if( groupI < 0 )
forAll(f, eI)
{
groupOwner[faceI] = -1;
continue;
const edge fe = f.faceEdge(eI);
if( fe == edges[beI] && mesh.faceIsInProcPatch(c[fI]) >= 0 )
{
validCell = true;
break;
}
}
groupOwner[faceI] = localGroupLabel[groupI];
if( validCell )
break;
}
OPstream toOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo(),
groupOwner.byteSize()
);
if( !validCell )
continue;
const label groupI = elementInGroup[edgeFaces(beI, 0)];
if( groupI < 0 )
continue;
toOtherProc << groupOwner;
const label lgI = localGroupLabel[groupI];
exchangeData[otherProc[beI]].append(labelPair(it.key(), lgI));
}
//- receive data from other processors
forAll(procBoundaries, patchI)
{
const label start = procBoundaries[patchI].patchStart();
LongList<labelPair> receivedData;
help::exchangeMap(exchangeData, receivedData);
labelList receivedData;
forAll(receivedData, i)
{
const labelPair& lp = receivedData[i];
IPstream fromOtherProc
(
Pstream::blocking,
procBoundaries[patchI].neiProcNo()
);
const label beI = globalToLocal[lp.first()];
fromOtherProc >> receivedData;
const cell& c = cells[faceOwner[edgeFaces(beI, 0)]];
forAll(receivedData, faceI)
//- combine data if the cell attached to this face has a face
//- attached to the inter-processor boundary
//- this must hold for boundary layer cells
bool validCell(false);
forAll(c, fI)
{
if( receivedData[faceI] < 0 )
continue;
const label groupI = elementInGroup[owner[start+faceI]];
const face& f = faces[c[fI]];
if( groupI < 0 )
continue;
forAll(f, eI)
{
const edge fe = f.faceEdge(eI);
DynList<label>& ng = neiGroups[localGroupLabel[groupI]];
if( fe == edges[beI] && mesh.faceIsInProcPatch(c[fI]) >= 0 )
{
validCell = true;
break;
}
}
//- store the connection over the inter-processor boundary
ng.appendIfNotIn(receivedData[faceI]);
if( validCell )
break;
}
if( !validCell )
continue;
const label groupI = elementInGroup[edgeFaces(beI, 0)];
if( groupI < 0 )
continue;
DynList<label>& ng = neiGroups[localGroupLabel[groupI]];
//- store the connection over the inter-processor boundary
ng.appendIfNotIn(lp.second());
}
}
};
class meshBndLayerSelectorOperator
{
const polyMeshGen& mesh_;
const meshSurfaceEngine& mse_;
public:
meshBndLayerSelectorOperator(const polyMeshGen& mesh)
meshBndLayerSelectorOperator(const meshSurfaceEngine& mse)
:
mesh_(mesh)
mse_(mse)
{}
bool operator()(const label cellI) const
bool operator()(const label bfI) const
{
const faceListPMG& faces = mesh_.faces();
const labelList& faceOwner = mse_.faceOwners();
const polyMeshGen& mesh = mse_.mesh();
const faceListPMG& faces = mesh.faces();
const cell& c = mesh_.cells()[cellI];
const PtrList<writePatch>& boundaries = mesh_.boundaries();
const cell& c = mesh.cells()[faceOwner[bfI]];
const PtrList<writePatch>& boundaries = mesh.boundaries();
const label start = boundaries[0].patchStart();
label nBndFaces(0), baseFace(-1), otherBase(-1), nQuads(0);
......@@ -224,59 +262,85 @@ public:
if( faces[c[fI]].size() == 4 )
++nQuads;
if( (c[fI] - start) >= 0 )
if( (c[fI] - start) == bfI )
{
baseFace = fI;
++nBndFaces;
}
else if( faces[c[fI]].size() != 4 )
{