33 void 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);
189 if (
topoMesh.isInternalFace(blockFaceLabel))
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];
246 cp.setSize(blockPfaceFacePoints.size());
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;
List< labelList > labelListList
A List of labelList.
errorManipArg< error, int > exit(error &err, const int errNo=1)
void resize(const label len)
Adjust allocated size of list.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
constexpr char nl
The newline '\n' character (0x0a)
List< face > faceList
A List of faces.
Ostream & endl(Ostream &os)
Add newline and flush stream.
bool cp(const fileName &src, const fileName &dst, const bool followLink=true)
Copy the source to the destination (recursively if necessary).
SubList< face > subList
Declare type of subList.
#define forAll(list, i)
Loop across all elements in list.
PtrList< block > blockList
The list of blocks is stored as a PtrList.
vectorField pointField
pointField is a vectorField.
void setSize(const label n)
Alias for resize()
const pointField & points() const
The points for the entire mesh. These points are scaled and transformed.
const polyMesh & topology() const
The blockMesh topology as a polyMesh unscaled and non-transformed.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Ostream & flush(Ostream &os)
Flush stream.
List< labelListList > labelListListList
A List of labelListList.
const polyMesh & topoMesh
messageStream Info
Information stream (stdout output on master, null elsewhere)
List< label > labelList
A List of labels.
List< cell > cellList
A List of cells.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)