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

A stable version of boundary layer refinement


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@37 fdcce57e-7e00-11e2-b579-49867b4cea03
parent 965e43d5
......@@ -140,11 +140,11 @@ class refineBoundaryLayers
//- find edges which shall be split due to refinement
bool findSplitEdges();
//- generate new vertices on edges, faces and in cells
//- generate new points on edges, faces and in cells
void generateNewVertices();
//- refine a given face and return the new faces
//- generates new vertices at cross-split faces
//- generates new points at cross-split faces
void refineFace
(
const face& f,
......@@ -152,6 +152,24 @@ class refineBoundaryLayers
DynList<DynList<label, 4>, 128>& newFaces
);
//- generate a matrix of points generated by splitting a face
//- and return them in the local i, j system of the face
void sortFacePoints
(
const label faceI,
DynList<DynList<label> >& facePoints,
const label transpose = false
) const;
//- generate a matrix of faces generated by splitting a face
//- and return them in the local i, j, system of the face
void sortFaceFaces
(
const label faceI,
DynList<DynList<label> >& faceFaces,
const label transpose = false
) const;
//- map split edges onto a cell
void generateNewFaces();
......@@ -162,23 +180,150 @@ class refineBoundaryLayers
DynList<DynList<DynList<label, 8>, 10> >& cellsFromCell
);
//- generate new cells from a hex at a feature edge
void generateNewCellsEdgeHex
//- a helper function which stores faces generated from
//- an existing face into new cells
void storeFacesIntoCells
(
const label cellI,
DynList<DynList<DynList<label, 4>, 6>, 64>& cellsFromCell
);
//- generate new cells from a hex at a corner
void generateNewCellsCornerHex
(
const label cellI,
const label faceI,
const bool reverseOrientation,
const label normalDirection,
const bool maxCoordinate,
const label nLayersI,
const label nLayersJ,
const label nLayersK,
DynList<DynList<DynList<label, 4>, 6>, 256>& cellsFromCell
);
) const;
//- generate new cells and add them to the mesh
void generateNewCells();
// Nested classes
class refineEdgeHexCell
{
// Private data
//- label of cell
const label cellI_;
//- number of cells in local direction i
label nLayersI_;
//- number of cells in locatiol direction j
label nLayersJ_;
//- container for new cells
DynList<DynList<DynList<label, 4>, 6>, 256> cellsFromCell_;
//- reference to the boundary layer class
refineBoundaryLayers& bndLayers_;
//- faces sorted into directions of a hex shape
FixedList<label, 6> faceInDirection_;
//- information about orientation of faces
//- false means the orientation as expected
//- true means wrong orientation
FixedList<bool, 6> faceOrientation_;
//- points on cross-split faces
FixedList<DynList<DynList<label> >, 2> cellPoints_;
// Private member functions
//- populate faceInDirection_nad wrongFaceOrientation_
void determineFacesInDirections();
//- populate new cells with new faces generated from already
//- existing faces
void populateExistingFaces();
//- generate new internal faces and tore them to new cells
void generateMissingFaces();
public:
// Constructor
//- construct from cell label and the refineBoundaryLayers
refineEdgeHexCell(const label cellI, refineBoundaryLayers& ref);
// Public member functions
inline const DynList<DynList<DynList<label, 4>, 6>, 256>&
newCells() const
{
return cellsFromCell_;
}
};
class refineCornerHexCell
{
// Private data
//- label of cell
const label cellI_;
//- number of cells in local direction i
label nLayersI_;
//- number of cells in local direction j
label nLayersJ_;
//- number of cells in local direction k
label nLayersK_;
//- split edge in directions
FixedList<label, 3> splitEdgeInDirection_;
//- container for new cells
DynList<DynList<DynList<label, 4>, 6>, 256> cellsFromCell_;
//- reference to the boundary layer class
refineBoundaryLayers& bndLayers_;
//- faces sorted into directions of a hex shape
FixedList<label, 6> faceInDirection_;
//- information about orientation of faces
//- false means the orientation as expected
//- true means wrong orientation
FixedList<bool, 6> faceOrientation_;
//- points on cross-split faces
FixedList<DynList<DynList<label> >, 6> facePoints_;
//- points inside the cell
DynList<DynList<DynList<label> > > cellPoints_;
// Private member functions
//- populate faceInDirection_nad wrongFaceOrientation_
void determineFacesInDirections();
//- populate new cells with new faces generated from already
//- existing faces
void populateExistingFaces();
//- generate missing points inside the cell
void generateNewPoints();
//- generate new internal faces and tore them to new cells
void generateMissingFaces();
public:
// Constructor
//- construct from cell label and the refineBoundaryLayers
refineCornerHexCell
(
const label cellI,
refineBoundaryLayers& ref
);
// Public member functions
inline const DynList<DynList<DynList<label, 4>, 6>, 256>&
newCells() const
{
return cellsFromCell_;
}
};
// Private member functions
//- Disallow bitwise copy construct
refineBoundaryLayers(const refineBoundaryLayers&);
......
......@@ -281,69 +281,457 @@ void refineBoundaryLayers::refineFace
{
for(label i=0;i<nLayersDir0;++i)
{
if( !((i == (nLayersDir0 - 1)) && (j == (nLayersDir1 - 1))) )
//- create quad face
DynList<label, 4> f;
f.append(facePoints[i][j]);
if( (i == (nLayersDir0 - 1)) && (j == 0) )
{
//- create quad face
DynList<label, 4> f;
f.setSize(4);
# ifdef DEBUGLayer
Info << "1. Adding additional points on edge " << endl;
# endif
f[0] = facePoints[i][j];
f[1] = facePoints[i+1][j];
f[2] = facePoints[i+1][j+1];
f[3] = facePoints[i][j+1];
//- add additional points on edge
const label eLabel = dir0Edges.second();
const label size =
newVerticesForSplitEdge_.sizeOfRow(eLabel) - 1;
newFaces.append(f);
for(label index=i+1;index<size;++index)
f.append(newVerticesForSplitEdge_(eLabel, index));
}
f.append(facePoints[i+1][j]);
if(
(dir1 != -1) &&
(i == (nLayersDir0 - 1)) &&
(j == (nLayersDir1 - 1))
)
{
# ifdef DEBUGLayer
Info << "2. Adding additional points on edge " << endl;
# endif
//- add additional points on edge
const label eLabel = dir1Edges.second();
const label size =
newVerticesForSplitEdge_.sizeOfRow(eLabel) - 1;
for(label index=j+1;index<size;++index)
f.append(newVerticesForSplitEdge_(eLabel, index));
}
f.append(facePoints[i+1][j+1]);
if( (i == (nLayersDir0 - 1)) && (j == (nLayersDir1 - 1)) )
{
# ifdef DEBUGLayer
Info << "3. Adding additional points on edge " << endl;
# endif
const label eLabel = dir0Edges.first();
const label size =
newVerticesForSplitEdge_.sizeOfRow(eLabel) - 2;
for(label index=size;index>i;--index)
f.append(newVerticesForSplitEdge_(eLabel, index));
}
f.append(facePoints[i][j+1]);
if( (dir1 != -1) && (i == 0) && (j == (nLayersDir1 - 1)) )
{
# ifdef DEBUGLayer
Info << "4. Adding additional points on edge " << endl;
# endif
const label eLabel = dir1Edges.first();
const label size =
newVerticesForSplitEdge_.sizeOfRow(eLabel) - 2;
for(label index=size;index>j;--index)
f.append(newVerticesForSplitEdge_(eLabel, index));
}
newFaces.append(f);
}
}
//- create the last face which may not be a quad
DynList<label, 4> newF;
# 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::sortFacePoints
(
const label faceI,
DynList<DynList<label> >& facePoints,
const label transpose
) const
{
const faceListPMG& faces = mesh_.faces();
const face& f = faces[faceI];
# ifdef DEBUGLayer
Info << "Creating matrix of points on a split face " << faceI << endl;
Info << "Face comprises of points " << f << endl;
Info << "New faces from face " << facesFromFace_.sizeOfRow(faceI) << endl;
# endif
if( (dir0 != -1) && (dir1 == -1) )
label procStart = mesh_.faces().size();
const PtrList<writeProcessorPatch>& procBoundaries = mesh_.procBoundaries();
if( Pstream::parRun() )
procStart = procBoundaries[0].patchStart();
if(
(faceI < procStart) ||
procBoundaries[mesh_.faceIsInProcPatch(faceI)].owner()
)
{
//- face is split in one direction, only
label eLabel = dir0Edges.second();
label size = newVerticesForSplitEdge_.sizeOfRow(eLabel);
//- orientation of new faces is the same as the face itself
//- start the procedure by finding the number of splits in
//- both i and j direction
label numSplitsI(1);
const label pos = f.which(newFaces_(facesFromFace_(faceI, 0), 0));
forAllRow(facesFromFace_, faceI, i)
{
const label nfI = facesFromFace_(faceI, i);
for(label i=nLayersDir0-1;i<size;++i)
newF.append(newVerticesForSplitEdge_(eLabel, i));
if( (numSplitsI == 1) && newFaces_.contains(nfI, f.nextLabel(pos)) )
{
numSplitsI = i + 1;
break;
}
}
const label numSplitsJ = (facesFromFace_.sizeOfRow(faceI) / numSplitsI);
# ifdef DEBUGLayer
Pout << "Pos " << pos << endl;
Pout << "Num splits in direction 0 " << numSplitsI << endl;
Pout << "Num splits in direction 1 " << numSplitsJ << endl;
# endif
facePoints.setSize(numSplitsI+1);
forAll(facePoints, i)
facePoints[i].setSize(numSplitsJ+1);
//- start filling in the matrix
forAllRow(facesFromFace_, faceI, fI)
{
const label nfI = facesFromFace_(faceI, fI);
const label i = fI % numSplitsI;
const label j = fI / numSplitsI;
# ifdef DEBUGLayer
Pout << "New face " << fI << " is " << newFaces_[nfI] << endl;
Pout << " i = " << i << endl;
Pout << " j = " << j << endl;
# endif
eLabel = dir0Edges.first();
size = newVerticesForSplitEdge_.sizeOfRow(eLabel);
for(label i=size;i>=nLayersDir0;--i)
newF.append(newVerticesForSplitEdge_(eLabel, i-1));
if( newFaces_.sizeOfRow(nfI) == 4 )
{
facePoints[i][j] = newFaces_(nfI, 0);
facePoints[i+1][j] = newFaces_(nfI, 1);
facePoints[i+1][j+1] = newFaces_(nfI, 2);
facePoints[i][j+1] = newFaces_(nfI, 3);
}
else
{
if( j == 0 )
{
forAllRow(newFaces_, nfI, pI)
if( f.which(newFaces_(nfI, pI)) >= 0 )
{
facePoints[i+1][0] = newFaces_(nfI, pI);
break;
}
}
else if( i == 0 )
{
forAllRow(newFaces_, nfI, pI)
if( f.which(newFaces_(nfI, pI)) >= 0 )
{
facePoints[0][j+1] = newFaces_(nfI, pI);
break;
}
}
else
{
forAllRow(newFaces_, nfI, pI)
if( f.which(newFaces_(nfI, pI)) >= 0 )
{
facePoints[i+1][j+1] = newFaces_(nfI, pI);
break;
}
}
}
}
# ifdef DEBUGLayer
Pout << "Generated matrix of points on face " << faceI
<< " is " << facePoints << endl;
# endif
}
else if( dir0 != -1 && dir1 != -1 )
else
{
//- face is split in both directions
newF.append(facePoints[nLayersDir0-1][nLayersDir1-1]);
//- this situation exists on inter-processor boundaries
//- on neighbour processor. i and j coordinates are reversed
label numSplitsJ(1);
const label pos = f.which(newFaces_(facesFromFace_(faceI, 0), 0));
forAllRow(facesFromFace_, faceI, j)
{
const label nfI = facesFromFace_(faceI, j);
if( (numSplitsJ == 1) && newFaces_.contains(nfI, f.prevLabel(pos)) )
{
numSplitsJ = j + 1;
break;
}
}
const label numSplitsI = (facesFromFace_.sizeOfRow(faceI) / numSplitsJ);
# ifdef DEBUGLayer
Pout << "2. Num splits in direction 0 " << numSplitsI << endl;
Pout << "2. Num splits in direction 1 " << numSplitsJ << endl;
# endif
facePoints.setSize(numSplitsI+1);
forAll(facePoints, i)
facePoints[i].setSize(numSplitsJ+1);
//- start filling in the matrix
forAllRow(facesFromFace_, faceI, fI)
{
const label nfI = facesFromFace_(faceI, fI);
//- add additional points on edge
label eLabel = dir1Edges.second();
label size = newVerticesForSplitEdge_.sizeOfRow(eLabel) - 1;
const label i = fI / numSplitsJ;
const label j = fI % numSplitsJ;
for(label i=nLayersDir1-1;i<size;++i)
newF.append(newVerticesForSplitEdge_(eLabel, i));
# ifdef DEBUGLayer
Pout << "2. New face " << fI << " is " << newFaces_[nfI] << endl;
Pout << "2. i = " << i << endl;
Pout << "2. j = " << j << endl;
# endif
//- add other corner
newF.append(facePoints[nLayersDir0][nLayersDir1]);
if( newFaces_.sizeOfRow(nfI) == 4 )
{
facePoints[i][j] = newFaces_(nfI, 0);
facePoints[i][j+1] = newFaces_(nfI, 1);
facePoints[i+1][j+1] = newFaces_(nfI, 2);
facePoints[i+1][j] = newFaces_(nfI, 3);
}
else
{
if( i == 0 )
{
forAllRow(newFaces_, nfI, pI)
if( f.which(newFaces_(nfI, pI)) >= 0 )
{
facePoints[0][j+1] = newFaces_(nfI, pI);
break;
}
}
else if( j == 0 )
{
forAllRow(newFaces_, nfI, pI)
if( f.which(newFaces_(nfI, pI)) >= 0 )
{
facePoints[i+1][0] = newFaces_(nfI, pI);
break;
}
}
else
{
forAllRow(newFaces_, nfI, pI)
if( f.which(newFaces_(nfI, pI)) >= 0 )
{
facePoints[i+1][j+1] = newFaces_(nfI, pI);
break;
}
}
}
}
//- add additional points on edge
eLabel = dir0Edges.first();
size = newVerticesForSplitEdge_.sizeOfRow(eLabel) - 1;
for(label i=size;i>=nLayersDir0;--i)
newF.append(newVerticesForSplitEdge_(eLabel, i-1));
# ifdef DEBUGLayer
Pout << "Generated matrix of points on processor face " << faceI
<< " is " << facePoints << endl;
# endif
}
newFaces.append(newF);
if( transpose )
{
DynList<DynList<label> > transposedFacePoints;
transposedFacePoints.setSize(facePoints[0].size());
forAll(transposedFacePoints, j)
transposedFacePoints[j].setSize(facePoints.size());
forAll(facePoints, i)
forAll(facePoints[i], j)
transposedFacePoints[j][i] = facePoints[i][j];
facePoints = transposedFacePoints;
# ifdef DEBUGLayer
Pout << "Transposed face points " << facePoints << endl;
# endif
}
}
void refineBoundaryLayers::sortFaceFaces
(
const label faceI,
DynList<DynList<label> >& faceFaces,
const label transpose
) const
{
const faceListPMG& faces = mesh_.faces();
const face& f = faces[faceI];
# ifdef DEBUGLayer
Info << "Input face " << f << endl;
Info << "Decomposed faces are " << newFaces << endl;
//if( (nLayersInDirection[0] > 1) && (nLayersInDirection[1] > 1) )
//::exit(1);
Info << "Creating matrix of faces on a split face " << faceI << endl;
Info << "Face comprises of points " << f << endl;
# endif
label procStart = mesh_.faces().size();
const PtrList<writeProcessorPatch>& procBoundaries = mesh_.procBoundaries();
if( Pstream::parRun() )
procStart = procBoundaries[0].patchStart();
if(
(faceI < procStart) ||
procBoundaries[mesh_.faceIsInProcPatch(faceI)].owner()
)
{
//- orientation of new faces is the same as the face itself
//- start the procedure by finding the number of splits in
//- both i and j direction