64 for(
int facePi=0; facePi<6; facePi++)
66 for(
int faceNi=0; faceNi<6; faceNi++)
68 for(
int rot=0; rot<4; rot++)
72 for(
int Pp=0; Pp<2; Pp++)
75 int Np = (3 - Pp + rot)%4;
81 if (
mag(map[0]) == 2 && map[0]*map[1] < 0)
102 for(
int rot=0; rot<4; rot++)
104 if (faceN[rot] == faceP[0])
111 <<
"Cannot find point correspondence for faces "
112 << faceP <<
" and " << faceN
128 forAll(topoInternalFaces, topoFacei)
130 label topoPi = topoFaceCell[topoFacei];
131 const labelList& topoPfaces = topoCells[topoPi];
133 bool foundFace =
false;
138 topoPfacei < topoPfaces.
size();
142 if (topoPfaces[topoPfacei] == topoFacei)
152 <<
"Cannot find merge face for block " << topoPi
156 mergeBlock[topoFacei].first() = topoPi;
157 mergeBlock[topoFacei].second() = topoPfacei;
192 return map < 0 ? -i-1 : i;
199 return i >= 0 ? i : ni + i + 1;
204inline label
mapij(
const int map,
const label i,
const label j)
304void Foam::blockMesh::calcTopologicalMerge()
313 Info<<
"Creating block offsets" <<
endl;
316 blockOffsets_.
resize(blocks.size());
323 blockOffsets_[blocki] = nPoints_;
325 nPoints_ += blocks[blocki].nPoints();
326 nCells_ += blocks[blocki].nCells();
331 Info<<
"Creating merge list (topological search).." <<
flush;
335 mergeList_.
setSize(nPoints_, -1);
355 List<Pair<label>> mergeBlockP(topoInternalFaces.size());
364 List<Pair<label>> mergeBlockN(topoInternalFaces.size());
375 forAll(topoInternalFaces, topoFacei)
379 label blockPi = mergeBlockP[topoFacei].first();
380 label blockPfacei = mergeBlockP[topoFacei].second();
382 label blockNi = mergeBlockN[topoFacei].first();
383 label blockNfacei = mergeBlockN[topoFacei].second();
390 blocks[blockPi].blockShape().
faces()[blockPfacei],
392 blocks[blockNi].blockShape().
faces()[blockNfacei]
397 <<
" Face map for faces "
398 << blocks[blockPi].blockShape().faces()[blockPfacei] <<
" "
399 << blocks[blockNi].blockShape().faces()[blockNfacei] <<
": "
402 const pointField& blockPpoints = blocks[blockPi].points();
403 const pointField& blockNpoints = blocks[blockNi].points();
405 Pair<label> Pnij(
faceNij(blockPfacei, blocks[blockPi]));
409 Pair<label> Nnij(
faceNij(blockNfacei, blocks[blockNi]));
411 NPnij[0] = Nnij[
mag(fmap[0]) - 1];
412 NPnij[1] = Nnij[
mag(fmap[1]) - 1];
417 <<
"Sub-division mismatch between face "
418 << blockPfacei <<
" of block " << blockPi << Pnij
420 << blockNfacei <<
" of block " << blockNi << Nnij
428 const boundBox bb(topoCells[blockPi].
points(topoFaces, topoPoints));
429 const scalar testSqrDist =
magSqr(1
e-6*bb.span());
432 scalar maxSqrDist = 0;
434 for (label j=0; j<Pnij.second(); j++)
436 for (label i=0; i<Pnij.first(); i++)
439 facePoint(blockPfacei, blocks[blockPi], i, j);
442 facePointN(blockNfacei, fmap, blocks[blockNi], i, j);
448 blockPpoints[blockPpointi]
449 - blockNpoints[blockNpointi]
453 if (sqrDist > testSqrDist)
456 <<
"Point merge failure between face "
457 << blockPfacei <<
" of block " << blockPi
459 << blockNfacei <<
" of block " << blockNi
461 <<
" Points: " << blockPpoints[blockPpointi]
462 <<
" " << blockNpoints[blockNpointi]
464 <<
" This may be due to inconsistent grading."
468 maxSqrDist =
max(maxSqrDist, sqrDist);
470 label Ppointi = blockPpointi + blockOffsets_[blockPi];
471 label Npointi = blockNpointi + blockOffsets_[blockNi];
473 label minPNi =
min(Ppointi, Npointi);
475 if (mergeList_[Ppointi] != -1)
477 minPNi =
min(minPNi, mergeList_[Ppointi]);
480 if (mergeList_[Npointi] != -1)
482 minPNi =
min(minPNi, mergeList_[Npointi]);
485 mergeList_[Ppointi] = mergeList_[Npointi] = minPNi;
490 <<
" Max distance between merge points: "
495 bool changedPointMerge =
false;
500 changedPointMerge =
false;
503 forAll(topoInternalFaces, topoFacei)
505 label blockPi = mergeBlockP[topoFacei].first();
506 label blockPfacei = mergeBlockP[topoFacei].second();
508 label blockNi = mergeBlockN[topoFacei].first();
509 label blockNfacei = mergeBlockN[topoFacei].second();
516 blocks[blockPi].blockShape().
faces()[blockPfacei],
518 blocks[blockNi].blockShape().
faces()[blockNfacei]
522 Pair<label> Pnij(
faceNij(blockPfacei, blocks[blockPi]));
524 for (label j=0; j<Pnij.second(); j++)
526 for (label i=0; i<Pnij.first(); i++)
529 facePoint(blockPfacei, blocks[blockPi], i, j);
532 facePointN(blockNfacei, fmap, blocks[blockNi], i, j);
535 blockPpointi + blockOffsets_[blockPi];
538 blockNpointi + blockOffsets_[blockNi];
540 if (mergeList_[Ppointi] != mergeList_[Npointi])
542 changedPointMerge =
true;
545 = mergeList_[Npointi]
546 =
min(mergeList_[Ppointi], mergeList_[Npointi]);
560 <<
"Point merging failed after 100 passes."
564 }
while (changedPointMerge);
572 label nUniqPoints = 0;
574 forAll(mergeList_, pointi)
576 if (mergeList_[pointi] > pointi)
579 <<
"Merge list contains point index out of range"
583 if (mergeList_[pointi] == -1 || mergeList_[pointi] == pointi)
585 mergeList_[pointi] = nUniqPoints++;
589 mergeList_[pointi] = mergeList_[mergeList_[pointi]];
594 nPoints_ = nUniqPoints;
const polyMesh & topoMesh
T & first() noexcept
The first element of the list, position [0].
SubList< face > subList
Declare type of subList.
void setSize(const label n)
Alias for resize()
void resize(const label len)
Adjust allocated size of list.
An ordered pair of two objects of type <T> with first() and second() elements.
const T & second() const noexcept
Return second element, which is also the last element.
A List obtained as a section of another List.
void size(const label n)
Older name for setAddressableSize.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
const labelVector & density() const noexcept
The mesh density (number of cells) in the i,j,k directions.
const blockFaceList & faces() const noexcept
The curved faces.
const pointField & points() const
const polyMesh & topology() const
PtrList< block > blockList
The list of blocks is stored as a PtrList.
Creates a single block of cells from point coordinates, numbers of cells in each direction and an exp...
A face is a list of labels corresponding to mesh vertices.
label pointLabel(const label i, const label j, const label k) const
The linear point index for an i-j-k position.
virtual const faceList & faces() const
Return raw faces.
virtual const labelList & faceOwner() const
Return face owner.
virtual const labelList & faceNeighbour() const
Return face neighbour.
virtual const pointField & points() const
Return raw points.
label nInternalFaces() const noexcept
Number of internal faces.
const cellList & cells() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
#define DebugInfo
Report an information message using Foam::Info.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
static const int faceEdgeDirs[6][4]
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
List< label > labelList
A List of labels.
List< cell > cellList
A List of cells.
label unsignIndex(const label i, const label ni)
label facePoint(const int facei, const block &block, const label i, const label j)
label facePointN(const block &block, const label i, const label j, const label k)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
label mapij(const int map, const label i, const label j)
Pair< label > faceNij(const label facei, const block &block)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
label signIndex(const int map, const label i)
void setBlockFaceCorrespondence(const cellList &topoCells, const faceList::subList &topoInternalFaces, const labelList &topoFaceCell, List< Pair< label > > &mergeBlock)
static Pair< int > faceFaceRotMap[6][6][4]
List< face > faceList
A List of faces.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Ostream & flush(Ostream &os)
Flush stream.
#define forAll(list, i)
Loop across all elements in list.