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

A stable mpi version of boundary layer refinement


git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@38 fdcce57e-7e00-11e2-b579-49867b4cea03
parent 9139db4a
......@@ -414,22 +414,22 @@ void refineBoundaryLayers::refineEdgeHexCell::determineFacesInDirections()
//- check the orientation of faces
const labelList& owner = mesh.owner();
//- checking face at direction 0
//- checking face at direction k = 0
faceOrientation_[0] = owner[c[faceInDirection_[0]]] == cellI_?true:false;
//- checking face in direction 1
//- checking face in direction k = 1
faceOrientation_[1] = owner[c[faceInDirection_[1]]] == cellI_?false:true;
//- set orientation flag for face in direction 2
//- set orientation flag for face in direction j = 0
faceOrientation_[2] = true;
//- checking face in direction 3
//- checking face in direction j = nLayersJ_
faceOrientation_[3] = owner[c[faceInDirection_[3]]] == cellI_?false:true;
//- set orientation flag for face in direction 4
//- set orientation flag for face in direction i = 0
faceOrientation_[4] = true;
//- checking face in direction 5
//- checking face in direction i = nLayersI_
faceOrientation_[5] = owner[c[faceInDirection_[5]]] == cellI_?false:true;
# ifdef DEBUGLayer
......@@ -790,10 +790,13 @@ void refineBoundaryLayers::refineCornerHexCell::determineFacesInDirections()
//- determine the directions of cell faces
//- store boundary faces first. Their normals point in the wrong direction
//- face at k = 0
faceInDirection_[0] = dirFace[permutation[0]];
faceOrientation_[0] = true;
//- face at j = 0
faceInDirection_[2] = dirFace[permutation[1]];
faceOrientation_[2] = true;
//- face at i = 0
faceInDirection_[4] = dirFace[permutation[2]];
faceOrientation_[4] = true;
......@@ -808,6 +811,7 @@ void refineBoundaryLayers::refineCornerHexCell::determineFacesInDirections()
if( !help::shareAnEdge(faces[c[fI]], faces[c[faceInDirection_[0]]]) )
{
//- face at k = nLayersK_
faceInDirection_[1] = fI;
faceOrientation_[1] = orientation;
}
......@@ -816,6 +820,7 @@ void refineBoundaryLayers::refineCornerHexCell::determineFacesInDirections()
!help::shareAnEdge(faces[c[fI]], faces[c[faceInDirection_[2]]])
)
{
//- face at j = nLayersJ_
faceInDirection_[3] = fI;
faceOrientation_[3] = orientation;
}
......@@ -824,6 +829,7 @@ void refineBoundaryLayers::refineCornerHexCell::determineFacesInDirections()
!help::shareAnEdge(faces[c[fI]], faces[c[faceInDirection_[4]]])
)
{
//- face at i = nLayersI_
faceInDirection_[5] = fI;
faceOrientation_[5] = orientation;
}
......@@ -1343,9 +1349,9 @@ void refineBoundaryLayers::generateNewCells()
# ifdef DEBUGLayer
forAll(nCellsFromCell, cellI)
{
Info << "\nCell " << cellI << endl;
Info << "nCellsFromCell " << nCellsFromCell[cellI] << endl;
Info << "Ref type " << refType[cellI] << endl;
Pout << "\nCell " << cellI << endl;
Pout << "nCellsFromCell " << nCellsFromCell[cellI] << endl;
Pout << "Ref type " << refType[cellI] << endl;
}
# endif
......@@ -1444,7 +1450,7 @@ void refineBoundaryLayers::generateNewCells()
const DynList<DynList<label, 4>, 6>& nc = cellsFromCell[cI];
# ifdef DEBUGLayer
Info << "Adding cell " << (cI==0?cellI:nCells)
Pout << "Adding cell " << (cI==0?cellI:nCells)
<< " originating from cell " << cellI << endl;
# endif
......@@ -1627,8 +1633,8 @@ void refineBoundaryLayers::generateNewCells()
# ifdef DEBUGLayer
Info << "Copying internal faces " << endl;
Info << "Original number of internal faces " << nOrigInternalFaces << endl;
Pout << "Copying internal faces " << endl;
Pout << "Original number of internal faces " << nOrigInternalFaces << endl;
# endif
//- store internal faces originating from existing faces
......@@ -1655,10 +1661,10 @@ void refineBoundaryLayers::generateNewCells()
//- store newly-generated internal faces
# ifdef DEBUGLayer
Info << "Copying newly generated internal faces" << endl;
Info << "nNewInternalFaces " << currFace << endl;
Info << "numFacesBefore " << numFacesBefore << endl;
Info << "Total number of faces " << newFaces_.size() << endl;
Pout << "Copying newly generated internal faces" << endl;
Pout << "nNewInternalFaces " << currFace << endl;
Pout << "numFacesBefore " << numFacesBefore << endl;
Pout << "Total number of faces " << newFaces_.size() << endl;
# endif
for(label faceI=numFacesBefore;faceI<newFaces_.size();++faceI)
......@@ -1675,10 +1681,10 @@ void refineBoundaryLayers::generateNewCells()
//- store new boundary faces
# ifdef DEBUGLayer
Info << "Copying boundary faces " << endl;
Info << "currFace " << currFace << endl;
Info << "Faces size " << faces.size() << endl;
Info << "Initial number of faces " << facesFromFace_.size() << endl;
Pout << "Copying boundary faces " << endl;
Pout << "currFace " << currFace << endl;
Pout << "Faces size " << faces.size() << endl;
Pout << "Initial number of faces " << facesFromFace_.size() << endl;
# endif
PtrList<writePatch>& boundaries = meshModifier.boundariesAccess();
......@@ -1720,7 +1726,7 @@ void refineBoundaryLayers::generateNewCells()
if( Pstream::parRun() )
{
# ifdef DEBUGLayer
Info << "Copying processor faces" << endl;
Pout << "Copying processor faces" << endl;
# endif
//- copy faces at inter-processor boundaries
......@@ -1763,8 +1769,8 @@ void refineBoundaryLayers::generateNewCells()
}
# ifdef DEBUGLayer
Info << "Faces after refinement " << faces << endl;
Info << "newFaceLabel " << newFaceLabel << endl;
Pout << "Faces after refinement " << faces << endl;
Pout << "newFaceLabel " << newFaceLabel << endl;
# endif
//- update face subsets
......@@ -1774,7 +1780,7 @@ void refineBoundaryLayers::generateNewCells()
//- update cells to match the faces
# ifdef DEBUGLayer
Info << "Updating cells to match new faces" << endl;
Pout << "Updating cells to match new faces" << endl;
# endif
forAll(cells, cellI)
......@@ -1786,7 +1792,7 @@ void refineBoundaryLayers::generateNewCells()
}
# ifdef DEBUGLayer
Info << "Cleaning mesh " << endl;
Pout << "Cleaning mesh " << endl;
# endif
//- delete all adressing which is no longer up-to-date
......
......@@ -126,13 +126,13 @@ void refineBoundaryLayers::refineFace
}
# ifdef DEBUGLayer
Info << "Refining face " << f << endl;
Info << "Splits in direction " << nLayersInDirection << endl;
Info << "Here " << endl;
Info << "Dir0 " << dir0 << endl;
Info << "dir0Edges " << dir0Edges << endl;
Info << "Dir1 " << dir1 << endl;
Info << "dir1Edges " << dir1Edges << endl;
Pout << "Refining face " << f << endl;
Pout << "Splits in direction " << nLayersInDirection << endl;
Pout << "Here " << endl;
Pout << "Dir0 " << dir0 << endl;
Pout << "dir0Edges " << dir0Edges << endl;
Pout << "Dir1 " << dir1 << endl;
Pout << "dir1Edges " << dir1Edges << endl;
# endif
//- in case of only one refinement direction, it must direction 0
......@@ -159,29 +159,29 @@ void refineBoundaryLayers::refineFace
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;
Pout << "Face has points " << f << endl;
Pout << "dirEdges0 " << dir0Edges << endl;
Pout << "dir1Edges " << dir1Edges << endl;
if( dir0 >= 0 )
{
Info << "Points on edge " << dir0Edges.first() << " with nodes "
Pout << "Points on edge " << dir0Edges.first() << " with nodes "
<< splitEdges_[dir0Edges.first()]
<< " are " << newVerticesForSplitEdge_[dir0Edges.first()] << endl;
Info << "Points on edge " << dir0Edges.second() << " with nodes "
Pout << "Points on edge " << dir0Edges.second() << " with nodes "
<< splitEdges_[dir0Edges.second()]
<< " are " << newVerticesForSplitEdge_[dir0Edges.second()] << endl;
}
if( dir1 >= 0 )
{
Info << "Points on edge " << dir1Edges.first() << " with nodes "
Pout << "Points on edge " << dir1Edges.first() << " with nodes "
<< splitEdges_[dir1Edges.first()]
<< " are " << newVerticesForSplitEdge_[dir1Edges.first()] << endl;
Info << "Points on edge " << dir1Edges.second() << " with nodes "
Pout << "Points on edge " << dir1Edges.second() << " with nodes "
<< splitEdges_[dir1Edges.second()]
<< " are " << newVerticesForSplitEdge_[dir1Edges.second()] << endl;
}
Info << "nLayersDir0 " << nLayersDir0 << endl;
Info << "nLayersDir1 " << nLayersDir1 << endl;
Pout << "nLayersDir0 " << nLayersDir0 << endl;
Pout << "nLayersDir1 " << nLayersDir1 << endl;
# endif
//- map the face onto a matrix for easier orientation
......@@ -211,7 +211,20 @@ void refineBoundaryLayers::refineFace
}
//- create missing vertices if there are any
const pointFieldPMG& points = mesh_.points();
pointFieldPMG& points = mesh_.points();
const point v00 = points[facePoints[0][0]];
const point v10 = points[facePoints[nLayersDir0][0]];
const point v01 = points[facePoints[0][nLayersDir1]];
const point v11 = points[facePoints[nLayersDir0][nLayersDir1]];
# ifdef DEBUGLayer
forAll(points, pointI)
Pout << "Point " << pointI << " coordinates " << points[pointI] << endl;
Pout << "v00 = " << v00 << endl;
Pout << "v10 = " << v10 << endl;
Pout << "v11 = " << v11 << endl;
Pout << "v01 = " << v01 << endl;
# endif
forAll(facePoints, i)
{
......@@ -220,60 +233,64 @@ 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;
Pout << "Determining u " << facePoints[0][0]
<< " coordinates " << points[facePoints[0][0]] << endl;
Pout << "Other point " << facePoints[i][0]
<< " coordinates " << points[facePoints[i][0]] << endl;
Pout << "Points at aplit edge "
<< newVerticesForSplitEdge_[dir0Edges.second()] << endl;}
# endif
const scalar u
(
mag(points[facePoints[i][0]] - points[facePoints[0][0]]) /
Foam::mag(points[facePoints[i][0]] - v00) /
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;
Pout << "Determining v " << facePoints[0][0]
<< " coordinates " << points[facePoints[0][0]] << endl;
Pout << "Other point " << facePoints[0][j]
<< " coordinates " << points[facePoints[0][j]] << endl;
Pout << "Points at aplit edge "
<< newVerticesForSplitEdge_[dir1Edges.first()] << endl;}
# endif
const scalar v
(
mag(points[facePoints[0][j]] - points[facePoints[0][0]]) /
Foam::mag(points[facePoints[0][j]] - v00) /
splitEdges_[dir1Edges.first()].mag(points)
);
# ifdef DEBUGLayer
Info << "Generating point of face " << endl;
Info << "u = " << u << endl;
Info << "v = " << v << endl;
Pout << "Generating point of face " << endl;
Pout << "u = " << u << endl;
Pout << "v = " << v << endl;}
# endif
//- calculate the coordinates of the missing point via
//- transfinite interpolation
const point newP
(
(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]]
(1.0 - u) * (1.0 - v) * v00 +
u * (1.0 - v) * v10 +
u * v * v11 +
(1.0 - u) * v * v01
);
# ifdef DEBUGLayer
Info << "Point coordinate " << newP << endl;
Pout << "Point coordinate " << newP << endl;
# endif
//- add the vertex to the mesh
facePoints[i][j] = points.size();
mesh_.appendVertex(newP);
points.append(newP);
}
}
}
# ifdef DEBUGLayer
Info << "Face points after creating vertices " << facePoints << endl;
Pout << "Face points after creating vertices " << facePoints << endl;
# endif
//- Finally, create the faces
......@@ -289,7 +306,7 @@ void refineBoundaryLayers::refineFace
if( (i == (nLayersDir0 - 1)) && (j == 0) )
{
# ifdef DEBUGLayer
Info << "1. Adding additional points on edge " << endl;
Pout << "1. Adding additional points on edge " << endl;
# endif
//- add additional points on edge
......@@ -310,7 +327,7 @@ void refineBoundaryLayers::refineFace
)
{
# ifdef DEBUGLayer
Info << "2. Adding additional points on edge " << endl;
Pout << "2. Adding additional points on edge " << endl;
# endif
//- add additional points on edge
......@@ -327,7 +344,7 @@ void refineBoundaryLayers::refineFace
if( (i == (nLayersDir0 - 1)) && (j == (nLayersDir1 - 1)) )
{
# ifdef DEBUGLayer
Info << "3. Adding additional points on edge " << endl;
Pout << "3. Adding additional points on edge " << endl;
# endif
const label eLabel = dir0Edges.first();
......@@ -342,7 +359,7 @@ void refineBoundaryLayers::refineFace
if( (dir1 != -1) && (i == 0) && (j == (nLayersDir1 - 1)) )
{
# ifdef DEBUGLayer
Info << "4. Adding additional points on edge " << endl;
Pout << "4. Adding additional points on edge " << endl;
# endif
const label eLabel = dir1Edges.first();
......@@ -357,8 +374,8 @@ void refineBoundaryLayers::refineFace
}
# ifdef DEBUGLayer
Info << "Input face " << f << endl;
Info << "Decomposed faces are " << newFaces << endl;
Pout << "Input face " << f << endl;
Pout << "Decomposed faces are " << newFaces << endl;
//if( (nLayersInDirection[0] > 1) && (nLayersInDirection[1] > 1) )
//::exit(1);
# endif
......@@ -375,9 +392,9 @@ void refineBoundaryLayers::sortFacePoints
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;
Pout << "Creating matrix of points on a split face " << faceI << endl;
Pout << "Face comprises of points " << f << endl;
Pout << "New faces from face " << facesFromFace_.sizeOfRow(faceI) << endl;
# endif
label procStart = mesh_.faces().size();
......@@ -500,6 +517,7 @@ void refineBoundaryLayers::sortFacePoints
const label numSplitsI = (facesFromFace_.sizeOfRow(faceI) / numSplitsJ);
# ifdef DEBUGLayer
Pout << "2. Face comprises of points " << f << endl;
Pout << "2. Num splits in direction 0 " << numSplitsI << endl;
Pout << "2. Num splits in direction 1 " << numSplitsJ << endl;
# endif
......@@ -525,9 +543,9 @@ void refineBoundaryLayers::sortFacePoints
if( newFaces_.sizeOfRow(nfI) == 4 )
{
facePoints[i][j] = newFaces_(nfI, 0);
facePoints[i][j+1] = newFaces_(nfI, 1);
facePoints[i+1][j] = newFaces_(nfI, 1);
facePoints[i+1][j+1] = newFaces_(nfI, 2);
facePoints[i+1][j] = newFaces_(nfI, 3);
facePoints[i][j+1] = newFaces_(nfI, 3);
}
else
{
......@@ -597,8 +615,8 @@ void refineBoundaryLayers::sortFaceFaces
const face& f = faces[faceI];
# ifdef DEBUGLayer
Info << "Creating matrix of faces on a split face " << faceI << endl;
Info << "Face comprises of points " << f << endl;
Pout << "Creating matrix of faces on a split face " << faceI << endl;
Pout << "Face comprises of points " << f << endl;
# endif
label procStart = mesh_.faces().size();
......@@ -814,10 +832,10 @@ void refineBoundaryLayers::generateNewFaces()
}
# ifdef DEBUGLayer
Info << "Internal face " << faceI << " with points " << f
Pout << "Internal face " << faceI << " with points " << f
<< " is refined " << endl;
forAllRow(facesFromFace_, faceI, i)
Info << "New face " << i << " is "
Pout << "New face " << i << " is "
<< newFaces_[facesFromFace_(faceI, i)] << endl;
DynList<DynList<label> > tralala;
sortFacePoints(faceI, tralala);
......@@ -879,10 +897,10 @@ void refineBoundaryLayers::generateNewFaces()
}
# ifdef DEBUGLayer
Info << "Boundary face " << faceI << " with points " << bf
Pout << "Boundary face " << faceI << " with points " << bf
<< " owner cell " << mesh_.owner()[faceI] << " is refined " << endl;
forAllRow(facesFromFace_, faceI, i)
Info << "New face " << i << " is "
Pout << "New face " << i << " is "
<< newFaces_[facesFromFace_(faceI, i)] << endl;
# endif
}
......@@ -998,7 +1016,7 @@ void refineBoundaryLayers::generateNewFaces()
# ifdef DEBUGLayer
returnReduce(1, sumOp<label>());
Info << "Starting splitting processor boundaries" << endl;
Pout << "Starting splitting processor boundaries" << endl;
# endif
//- perform splitting
......@@ -1026,7 +1044,7 @@ void refineBoundaryLayers::generateNewFaces()
if( procBoundaries[patchI].owner() )
{
//- this processor owns this patch
FixedList<label, 2> nLayersInDirection;
FixedList<label, 2> nLayersInDirection(1);
const DynList<labelPair, 2>& dirSplits = it->second;
forAll(dirSplits, i)
nLayersInDirection[dirSplits[i].first()] =
......@@ -1054,7 +1072,7 @@ void refineBoundaryLayers::generateNewFaces()
else
{
//- reverse the face before splitting
FixedList<label, 2> nLayersInDirection;
FixedList<label, 2> nLayersInDirection(1);
const DynList<labelPair, 2>& dirSplits = it->second;
forAll(dirSplits, i)
nLayersInDirection[(dirSplits[i].first()+1)%2] =
......
......@@ -759,24 +759,19 @@ void refineBoundaryLayers::generateNewVertices()
const label nThreads = 1;
# endif
DynList<label> numPointsAtThread;
numPointsAtThread.setSize(nThreads);
# ifdef USE_OMP
# pragma omp parallel num_threads(nThreads)
# endif
{
# ifdef USE_OMP
const label threadI = omp_get_thread_num();
# else
const label threadI(0);
# endif
label& nPoints = numPointsAtThread[threadI];
nPoints = 0;
// # ifdef USE_OMP
// const label threadI = omp_get_thread_num();
// # else
// const label threadI(0);
// # endif
//- start counting vertices at each thread
# ifdef USE_OMP
# pragma omp for schedule(static)
# pragma omp for schedule(static, 1)
# endif
forAll(splitEdges_, seI)
{
......@@ -833,7 +828,6 @@ void refineBoundaryLayers::generateNewVertices()
firstLayerThickness[seI] = thickness;
thicknessRatio[seI] = ratio;
nNodesAtEdge[seI] = nLayers + 1;
nPoints += nLayers - 1;
}
}
......@@ -983,9 +977,31 @@ void refineBoundaryLayers::generateNewVertices()
}
}
//- allocate the memory
//- calculate the number of additional vertices which will be generated
//- on edges of the mesh
DynList<label> numPointsAtThread;
numPointsAtThread.setSize(nThreads);
numPointsAtThread = 0;
# ifdef USE_OMP
# pragma omp parallel for num_threads(nThreads) schedule(static, 1)
# endif
forAll(nNodesAtEdge, seI)
{
# ifdef USE_OMP
const label threadI = omp_get_thread_num();
# else
const label threadI(0);
# endif
numPointsAtThread[threadI] += nNodesAtEdge[seI] - 2;
}
//- allocate the space in a graph storing ids of points on a split edge
newVerticesForSplitEdge_.setSizeAndRowSize(nNodesAtEdge);
//- calculate the number of points which will be generated
//- on split edges
label numPoints = points.size();
forAll(numPointsAtThread, threadI)
{
......@@ -1000,7 +1016,7 @@ void refineBoundaryLayers::generateNewVertices()
Info << "Generating split vertices" << endl;
# endif
//- generate vertices at split edges
//- generate vertices on split edges
# ifdef USE_OMP
# pragma omp parallel num_threads(nThreads)
# endif
......@@ -1014,7 +1030,7 @@ void refineBoundaryLayers::generateNewVertices()
label& nPoints = numPointsAtThread[threadI];
# ifdef USE_OMP
# pragma omp for schedule(static)
# pragma omp for schedule(static, 1)
# endif
forAll(splitEdges_, seI)
{
......@@ -1036,8 +1052,12 @@ void refineBoundaryLayers::generateNewVertices()
);
# ifdef DEBUGLayer
Info << "Edge length " << magv << endl;
Info << "Thickness of the first layer "
Pout << "Thread " << threadI << endl;
Pout << "Generating vertices at split edge "
<< " start point " << points[e.start()]
<< " end point " << points[e.end()] << endl;
Pout << "Edge length " << magv << endl;
Pout << "Thickness of the first layer "
<< firstThickness << endl;
# endif
}
......@@ -1061,9 +1081,13 @@ void refineBoundaryLayers::generateNewVertices()
const point newP = points[e.start()] + param * vec;
# ifdef DEBUGLayer
Info << "Split edge " << seI << " points " << e
Pout << "Split edge " << seI << " edge points " << e
<< " start point " << points[e.start()]
<< " end point " << points[e.end()]
<< " param " << param
<< " first thickness " << firstThickness << endl;
<< " new point " << nPoints
<< " has coordinates " << newP << endl;
}
# endif
param += firstThickness * Foam::pow(thicknessRatio[seI], pI);
......