33void Foam::blockMesh::calcGeometricalMerge()
39 Info<<
"Creating block offsets" <<
endl;
42 blockOffsets_.
resize(blocks.size());
49 blockOffsets_[blocki] = nPoints_;
51 nPoints_ += blocks[blocki].nPoints();
52 nCells_ += blocks[blocki].nCells();
58 Info<<
"Creating merge list (geometric search).." <<
flush;
62 mergeList_.
resize(nPoints_);
84 forAll(blockFaces, blockFaceLabel)
86 label blockPlabel = faceOwnerBlocks[blockFaceLabel];
87 const pointField& blockPpoints = blocks[blockPlabel].points();
88 const labelList& blockPfaces = blockCells[blockPlabel];
90 bool foundFace =
false;
91 label blockPfaceLabel;
95 blockPfaceLabel < blockPfaces.size();
99 if (blockPfaces[blockPfaceLabel] == blockFaceLabel)
109 <<
"Cannot find merge face for block " << blockPlabel
113 const List<FixedList<label, 4>>& blockPfaceFaces =
114 blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
117 curPairs.
setSize(blockPfaceFaces.size());
126 boundBox bb(blockCells[blockPlabel].
points(blockFaces, blockPoints));
127 const scalar mergeSqrDist =
magSqr(10*SMALL*bb.span());
131 scalar sqrMergeTol = GREAT;
133 forAll(blockPfaceFaces, blockPfaceFaceLabel)
135 const FixedList<label, 4>& blockPfaceFacePoints
136 = blockPfaceFaces[blockPfaceFaceLabel];
138 forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
140 forAll(blockPfaceFacePoints, blockPfaceFacePointLabel2)
142 if (blockPfaceFacePointLabel != blockPfaceFacePointLabel2)
144 scalar magSqrDist =
magSqr
146 blockPpoints[blockPfaceFacePoints
147 [blockPfaceFacePointLabel]]
148 - blockPpoints[blockPfaceFacePoints
149 [blockPfaceFacePointLabel2]]
152 if (magSqrDist < mergeSqrDist)
155 blockPfaceFacePoints[blockPfaceFacePointLabel]
156 + blockOffsets_[blockPlabel];
159 blockPfaceFacePoints[blockPfaceFacePointLabel2]
160 + blockOffsets_[blockPlabel];
162 label minPP2 =
min(PpointLabel, PpointLabel2);
164 if (mergeList_[PpointLabel] != -1)
166 minPP2 =
min(minPP2, mergeList_[PpointLabel]);
169 if (mergeList_[PpointLabel2] != -1)
171 minPP2 =
min(minPP2, mergeList_[PpointLabel2]);
174 mergeList_[PpointLabel] = mergeList_[PpointLabel2]
179 sqrMergeTol =
min(sqrMergeTol, magSqrDist);
191 label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
192 const pointField& blockNpoints = blocks[blockNlabel].points();
193 const labelList& blockNfaces = blockCells[blockNlabel];
196 label blockNfaceLabel;
200 blockNfaceLabel < blockNfaces.size();
206 blockFaces[blockNfaces[blockNfaceLabel]]
207 == blockFaces[blockFaceLabel]
218 <<
"Cannot find merge face for block " << blockNlabel
222 const List<FixedList<label, 4>>& blockNfaceFaces =
223 blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
227 checkFaceCorrespondence_
228 && blockPfaceFaces.size() != blockNfaceFaces.size()
232 <<
"Inconsistent number of faces between block pair "
233 << blockPlabel <<
" and " << blockNlabel
240 forAll(blockPfaceFaces, blockPfaceFaceLabel)
242 const FixedList<label, 4>& blockPfaceFacePoints
243 = blockPfaceFaces[blockPfaceFaceLabel];
249 forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
251 forAll(blockNfaceFaces, blockNfaceFaceLabel)
253 const FixedList<label, 4>& blockNfaceFacePoints
254 = blockNfaceFaces[blockNfaceFaceLabel];
256 forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
263 [blockPfaceFacePoints[blockPfaceFacePointLabel]]
265 [blockNfaceFacePoints[blockNfaceFacePointLabel]]
271 cp[blockPfaceFacePointLabel] =
272 blockNfaceFacePoints[blockNfaceFacePointLabel];
275 blockPfaceFacePoints[blockPfaceFacePointLabel]
276 + blockOffsets_[blockPlabel];
279 blockNfaceFacePoints[blockNfaceFacePointLabel]
280 + blockOffsets_[blockNlabel];
282 label minPN =
min(PpointLabel, NpointLabel);
284 if (mergeList_[PpointLabel] != -1)
286 minPN =
min(minPN, mergeList_[PpointLabel]);
289 if (mergeList_[NpointLabel] != -1)
291 minPN =
min(minPN, mergeList_[NpointLabel]);
294 mergeList_[PpointLabel] = mergeList_[NpointLabel]
300 forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
302 if (
cp[blockPfaceFacePointLabel] == -1)
305 <<
"Inconsistent point locations between block pair "
306 << blockPlabel <<
" and " << blockNlabel <<
nl
307 <<
" probably due to inconsistent grading."
316 bool changedPointMerge =
false;
321 changedPointMerge =
false;
324 forAll(blockinternalFaces, blockFaceLabel)
326 label blockPlabel = faceOwnerBlocks[blockFaceLabel];
327 label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
329 const labelList& blockPfaces = blockCells[blockPlabel];
330 const labelList& blockNfaces = blockCells[blockNlabel];
332 const labelListList& curPairs = glueMergePairs[blockFaceLabel];
334 label blockPfaceLabel;
338 blockPfaceLabel < blockPfaces.size();
344 blockFaces[blockPfaces[blockPfaceLabel]]
345 == blockinternalFaces[blockFaceLabel]
353 label blockNfaceLabel;
357 blockNfaceLabel < blockNfaces.size();
363 blockFaces[blockNfaces[blockNfaceLabel]]
364 == blockinternalFaces[blockFaceLabel]
372 const List<FixedList<label, 4>>& blockPfaceFaces =
373 blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
375 forAll(blockPfaceFaces, blockPfaceFaceLabel)
377 const FixedList<label, 4>& blockPfaceFacePoints
378 = blockPfaceFaces[blockPfaceFaceLabel];
380 const labelList&
cp = curPairs[blockPfaceFaceLabel];
382 forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
385 blockPfaceFacePoints[blockPfaceFacePointLabel]
386 + blockOffsets_[blockPlabel];
389 cp[blockPfaceFacePointLabel]
390 + blockOffsets_[blockNlabel];
394 mergeList_[PpointLabel]
395 != mergeList_[NpointLabel]
398 changedPointMerge =
true;
400 mergeList_[PpointLabel]
401 = mergeList_[NpointLabel]
404 mergeList_[PpointLabel],
405 mergeList_[NpointLabel]
419 <<
"Point merging failed after max number of passes."
423 while (changedPointMerge);
430 forAll(blockinternalFaces, blockFaceLabel)
432 label blockPlabel = faceOwnerBlocks[blockFaceLabel];
433 label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
435 const labelList& blockPfaces = blockCells[blockPlabel];
436 const labelList& blockNfaces = blockCells[blockNlabel];
438 const pointField& blockPpoints = blocks[blockPlabel].points();
439 const pointField& blockNpoints = blocks[blockNlabel].points();
441 bool foundFace =
false;
442 label blockPfaceLabel;
446 blockPfaceLabel < blockPfaces.size();
452 blockFaces[blockPfaces[blockPfaceLabel]]
453 == blockinternalFaces[blockFaceLabel]
464 <<
"Cannot find merge face for block " << blockPlabel
469 label blockNfaceLabel;
473 blockNfaceLabel < blockNfaces.size();
479 blockFaces[blockNfaces[blockNfaceLabel]]
480 == blockinternalFaces[blockFaceLabel]
491 <<
"Cannot find merge face for block " << blockNlabel
495 const List<FixedList<label, 4>>& blockPfaceFaces =
496 blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
498 const List<FixedList<label, 4>>& blockNfaceFaces =
499 blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
501 forAll(blockPfaceFaces, blockPfaceFaceLabel)
503 const FixedList<label, 4>& blockPfaceFacePoints
504 = blockPfaceFaces[blockPfaceFaceLabel];
506 forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
509 blockPfaceFacePoints[blockPfaceFacePointLabel]
510 + blockOffsets_[blockPlabel];
512 if (mergeList_[PpointLabel] == -1)
515 <<
"Unable to merge point "
516 << blockPfaceFacePointLabel
517 <<
' ' << blockPpoints[blockPfaceFacePointLabel]
527 forAll(blockNfaceFaces, blockNfaceFaceLabel)
529 const FixedList<label, 4>& blockNfaceFacePoints
530 = blockNfaceFaces[blockNfaceFaceLabel];
532 forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
535 blockNfaceFacePoints[blockNfaceFacePointLabel]
536 + blockOffsets_[blockNlabel];
538 if (mergeList_[NpointLabel] == -1)
541 <<
"unable to merge point "
542 << blockNfaceFacePointLabel
543 <<
' ' << blockNpoints[blockNfaceFacePointLabel]
557 label newPointLabel = 0;
559 forAll(mergeList_, pointLabel)
561 if (mergeList_[pointLabel] > pointLabel)
564 <<
"Merge list contains point index out of range"
570 mergeList_[pointLabel] == -1
571 || mergeList_[pointLabel] == pointLabel
574 mergeList_[pointLabel] = newPointLabel;
579 mergeList_[pointLabel] = mergeList_[mergeList_[pointLabel]];
583 nPoints_ = newPointLabel;
const polyMesh & topoMesh
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.
const pointField & points() const
const polyMesh & topology() const
PtrList< block > blockList
The list of blocks is stored as a PtrList.
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.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
label nInternalFaces() const noexcept
Number of internal faces.
const cellList & cells() const
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
List< label > labelList
A List of labels.
List< cell > cellList
A List of cells.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Ostream & endl(Ostream &os)
Add newline and flush stream.
List< labelList > labelListList
A List of labelList.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
List< labelListList > labelListListList
A List of labelListList.
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.
constexpr char nl
The newline '\n' character (0x0a)
const volScalarField & cp
#define forAll(list, i)
Loop across all elements in list.