blockMeshMergeGeometrical.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "blockMesh.H"
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 void Foam::blockMesh::calcGeometricalMerge()
34 {
35  const blockList& blocks = *this;
36 
37  if (verbose_)
38  {
39  Info<< "Creating block offsets" << endl;
40  }
41 
42  blockOffsets_.resize(blocks.size());
43 
44  nPoints_ = 0;
45  nCells_ = 0;
46 
47  forAll(blocks, blocki)
48  {
49  blockOffsets_[blocki] = nPoints_;
50 
51  nPoints_ += blocks[blocki].nPoints();
52  nCells_ += blocks[blocki].nCells();
53  }
54 
55 
56  if (verbose_)
57  {
58  Info<< "Creating merge list (geometric search).." << flush;
59  }
60 
61  // Set unused to -1
62  mergeList_.resize(nPoints_);
63  mergeList_ = -1;
64 
65 
66  // Block mesh topology
67  const polyMesh& topoMesh = topology();
68 
69  const pointField& blockPoints = topoMesh.points();
70  const cellList& blockCells = topoMesh.cells();
71  const faceList& blockFaces = topoMesh.faces();
72  const labelList& faceOwnerBlocks = topoMesh.faceOwner();
73  const labelList& faceNeighbourBlocks = topoMesh.faceNeighbour();
74 
75  const faceList::subList blockinternalFaces
76  (
77  blockFaces,
78  topoMesh.nInternalFaces()
79  );
80 
81  // For efficiency, create merge pairs in the first pass
82  labelListListList glueMergePairs(blockFaces.size());
83 
84  forAll(blockFaces, blockFaceLabel)
85  {
86  label blockPlabel = faceOwnerBlocks[blockFaceLabel];
87  const pointField& blockPpoints = blocks[blockPlabel].points();
88  const labelList& blockPfaces = blockCells[blockPlabel];
89 
90  bool foundFace = false;
91  label blockPfaceLabel;
92  for
93  (
94  blockPfaceLabel = 0;
95  blockPfaceLabel < blockPfaces.size();
96  blockPfaceLabel++
97  )
98  {
99  if (blockPfaces[blockPfaceLabel] == blockFaceLabel)
100  {
101  foundFace = true;
102  break;
103  }
104  }
105 
106  if (!foundFace)
107  {
108  FatalErrorInFunction
109  << "Cannot find merge face for block " << blockPlabel
110  << exit(FatalError);
111  }
112 
113  const List<FixedList<label, 4>>& blockPfaceFaces =
114  blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
115 
116  labelListList& curPairs = glueMergePairs[blockFaceLabel];
117  curPairs.setSize(blockPfaceFaces.size());
118 
119  // Calculate sqr of the merge tolerance as 1/10th of the min sqr
120  // point to point distance on the block face.
121  // At the same time merge collated points on the block's faces
122  // (removes boundary poles etc.)
123  // Collated points detected by initially taking a constant factor of
124  // the size of the block.
125 
126  boundBox bb(blockCells[blockPlabel].points(blockFaces, blockPoints));
127  const scalar mergeSqrDist = magSqr(10*SMALL*bb.span());
128 
129  // This is an N^2 algorithm
130 
131  scalar sqrMergeTol = GREAT;
132 
133  forAll(blockPfaceFaces, blockPfaceFaceLabel)
134  {
135  const FixedList<label, 4>& blockPfaceFacePoints
136  = blockPfaceFaces[blockPfaceFaceLabel];
137 
138  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
139  {
140  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel2)
141  {
142  if (blockPfaceFacePointLabel != blockPfaceFacePointLabel2)
143  {
144  scalar magSqrDist = magSqr
145  (
146  blockPpoints[blockPfaceFacePoints
147  [blockPfaceFacePointLabel]]
148  - blockPpoints[blockPfaceFacePoints
149  [blockPfaceFacePointLabel2]]
150  );
151 
152  if (magSqrDist < mergeSqrDist)
153  {
154  label PpointLabel =
155  blockPfaceFacePoints[blockPfaceFacePointLabel]
156  + blockOffsets_[blockPlabel];
157 
158  label PpointLabel2 =
159  blockPfaceFacePoints[blockPfaceFacePointLabel2]
160  + blockOffsets_[blockPlabel];
161 
162  label minPP2 = min(PpointLabel, PpointLabel2);
163 
164  if (mergeList_[PpointLabel] != -1)
165  {
166  minPP2 = min(minPP2, mergeList_[PpointLabel]);
167  }
168 
169  if (mergeList_[PpointLabel2] != -1)
170  {
171  minPP2 = min(minPP2, mergeList_[PpointLabel2]);
172  }
173 
174  mergeList_[PpointLabel] = mergeList_[PpointLabel2]
175  = minPP2;
176  }
177  else
178  {
179  sqrMergeTol = min(sqrMergeTol, magSqrDist);
180  }
181  }
182  }
183  }
184  }
185 
186  sqrMergeTol /= 10.0;
187 
188 
189  if (topoMesh.isInternalFace(blockFaceLabel))
190  {
191  label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
192  const pointField& blockNpoints = blocks[blockNlabel].points();
193  const labelList& blockNfaces = blockCells[blockNlabel];
194 
195  foundFace = false;
196  label blockNfaceLabel;
197  for
198  (
199  blockNfaceLabel = 0;
200  blockNfaceLabel < blockNfaces.size();
201  blockNfaceLabel++
202  )
203  {
204  if
205  (
206  blockFaces[blockNfaces[blockNfaceLabel]]
207  == blockFaces[blockFaceLabel]
208  )
209  {
210  foundFace = true;
211  break;
212  }
213  }
214 
215  if (!foundFace)
216  {
217  FatalErrorInFunction
218  << "Cannot find merge face for block " << blockNlabel
219  << exit(FatalError);
220  }
221 
222  const List<FixedList<label, 4>>& blockNfaceFaces =
223  blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
224 
225  if
226  (
227  checkFaceCorrespondence_
228  && blockPfaceFaces.size() != blockNfaceFaces.size()
229  )
230  {
231  FatalErrorInFunction
232  << "Inconsistent number of faces between block pair "
233  << blockPlabel << " and " << blockNlabel
234  << exit(FatalError);
235  }
236 
237 
238  // N-squared point search over all points of all faces of
239  // master block over all point of all faces of slave block
240  forAll(blockPfaceFaces, blockPfaceFaceLabel)
241  {
242  const FixedList<label, 4>& blockPfaceFacePoints
243  = blockPfaceFaces[blockPfaceFaceLabel];
244 
245  labelList& cp = curPairs[blockPfaceFaceLabel];
246  cp.setSize(blockPfaceFacePoints.size());
247  cp = -1;
248 
249  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
250  {
251  forAll(blockNfaceFaces, blockNfaceFaceLabel)
252  {
253  const FixedList<label, 4>& blockNfaceFacePoints
254  = blockNfaceFaces[blockNfaceFaceLabel];
255 
256  forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
257  {
258  if
259  (
260  magSqr
261  (
262  blockPpoints
263  [blockPfaceFacePoints[blockPfaceFacePointLabel]]
264  - blockNpoints
265  [blockNfaceFacePoints[blockNfaceFacePointLabel]]
266  ) < sqrMergeTol
267  )
268  {
269  // Found a new pair
270 
271  cp[blockPfaceFacePointLabel] =
272  blockNfaceFacePoints[blockNfaceFacePointLabel];
273 
274  label PpointLabel =
275  blockPfaceFacePoints[blockPfaceFacePointLabel]
276  + blockOffsets_[blockPlabel];
277 
278  label NpointLabel =
279  blockNfaceFacePoints[blockNfaceFacePointLabel]
280  + blockOffsets_[blockNlabel];
281 
282  label minPN = min(PpointLabel, NpointLabel);
283 
284  if (mergeList_[PpointLabel] != -1)
285  {
286  minPN = min(minPN, mergeList_[PpointLabel]);
287  }
288 
289  if (mergeList_[NpointLabel] != -1)
290  {
291  minPN = min(minPN, mergeList_[NpointLabel]);
292  }
293 
294  mergeList_[PpointLabel] = mergeList_[NpointLabel]
295  = minPN;
296  }
297  }
298  }
299  }
300  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
301  {
302  if (cp[blockPfaceFacePointLabel] == -1)
303  {
304  FatalErrorInFunction
305  << "Inconsistent point locations between block pair "
306  << blockPlabel << " and " << blockNlabel << nl
307  << " probably due to inconsistent grading."
308  << exit(FatalError);
309  }
310  }
311  }
312  }
313  }
314 
315 
316  bool changedPointMerge = false;
317  label nPasses = 0;
318 
319  do
320  {
321  changedPointMerge = false;
322  nPasses++;
323 
324  forAll(blockinternalFaces, blockFaceLabel)
325  {
326  label blockPlabel = faceOwnerBlocks[blockFaceLabel];
327  label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
328 
329  const labelList& blockPfaces = blockCells[blockPlabel];
330  const labelList& blockNfaces = blockCells[blockNlabel];
331 
332  const labelListList& curPairs = glueMergePairs[blockFaceLabel];
333 
334  label blockPfaceLabel;
335  for
336  (
337  blockPfaceLabel = 0;
338  blockPfaceLabel < blockPfaces.size();
339  blockPfaceLabel++
340  )
341  {
342  if
343  (
344  blockFaces[blockPfaces[blockPfaceLabel]]
345  == blockinternalFaces[blockFaceLabel]
346  )
347  {
348  break;
349  }
350  }
351 
352 
353  label blockNfaceLabel;
354  for
355  (
356  blockNfaceLabel = 0;
357  blockNfaceLabel < blockNfaces.size();
358  blockNfaceLabel++
359  )
360  {
361  if
362  (
363  blockFaces[blockNfaces[blockNfaceLabel]]
364  == blockinternalFaces[blockFaceLabel]
365  )
366  {
367  break;
368  }
369  }
370 
371 
372  const List<FixedList<label, 4>>& blockPfaceFaces =
373  blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
374 
375  forAll(blockPfaceFaces, blockPfaceFaceLabel)
376  {
377  const FixedList<label, 4>& blockPfaceFacePoints
378  = blockPfaceFaces[blockPfaceFaceLabel];
379 
380  const labelList& cp = curPairs[blockPfaceFaceLabel];
381 
382  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
383  {
384  label PpointLabel =
385  blockPfaceFacePoints[blockPfaceFacePointLabel]
386  + blockOffsets_[blockPlabel];
387 
388  label NpointLabel =
389  cp[blockPfaceFacePointLabel]
390  + blockOffsets_[blockNlabel];
391 
392  if
393  (
394  mergeList_[PpointLabel]
395  != mergeList_[NpointLabel]
396  )
397  {
398  changedPointMerge = true;
399 
400  mergeList_[PpointLabel]
401  = mergeList_[NpointLabel]
402  = min
403  (
404  mergeList_[PpointLabel],
405  mergeList_[NpointLabel]
406  );
407  }
408  }
409  }
410  }
411  if (verbose_)
412  {
413  Info<< '.' << flush;
414  }
415 
416  if (nPasses > 100)
417  {
418  FatalErrorInFunction
419  << "Point merging failed after max number of passes."
420  << exit(FatalError);
421  }
422  }
423  while (changedPointMerge);
424 
425  if (verbose_)
426  {
427  Info<< endl;
428  }
429 
430  forAll(blockinternalFaces, blockFaceLabel)
431  {
432  label blockPlabel = faceOwnerBlocks[blockFaceLabel];
433  label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
434 
435  const labelList& blockPfaces = blockCells[blockPlabel];
436  const labelList& blockNfaces = blockCells[blockNlabel];
437 
438  const pointField& blockPpoints = blocks[blockPlabel].points();
439  const pointField& blockNpoints = blocks[blockNlabel].points();
440 
441  bool foundFace = false;
442  label blockPfaceLabel;
443  for
444  (
445  blockPfaceLabel = 0;
446  blockPfaceLabel < blockPfaces.size();
447  blockPfaceLabel++
448  )
449  {
450  if
451  (
452  blockFaces[blockPfaces[blockPfaceLabel]]
453  == blockinternalFaces[blockFaceLabel]
454  )
455  {
456  foundFace = true;
457  break;
458  }
459  }
460 
461  if (!foundFace)
462  {
463  FatalErrorInFunction
464  << "Cannot find merge face for block " << blockPlabel
465  << exit(FatalError);
466  }
467 
468  foundFace = false;
469  label blockNfaceLabel;
470  for
471  (
472  blockNfaceLabel = 0;
473  blockNfaceLabel < blockNfaces.size();
474  blockNfaceLabel++
475  )
476  {
477  if
478  (
479  blockFaces[blockNfaces[blockNfaceLabel]]
480  == blockinternalFaces[blockFaceLabel]
481  )
482  {
483  foundFace = true;
484  break;
485  }
486  }
487 
488  if (!foundFace)
489  {
490  FatalErrorInFunction
491  << "Cannot find merge face for block " << blockNlabel
492  << exit(FatalError);
493  }
494 
495  const List<FixedList<label, 4>>& blockPfaceFaces =
496  blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
497 
498  const List<FixedList<label, 4>>& blockNfaceFaces =
499  blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
500 
501  forAll(blockPfaceFaces, blockPfaceFaceLabel)
502  {
503  const FixedList<label, 4>& blockPfaceFacePoints
504  = blockPfaceFaces[blockPfaceFaceLabel];
505 
506  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
507  {
508  label PpointLabel =
509  blockPfaceFacePoints[blockPfaceFacePointLabel]
510  + blockOffsets_[blockPlabel];
511 
512  if (mergeList_[PpointLabel] == -1)
513  {
514  FatalErrorInFunction
515  << "Unable to merge point "
516  << blockPfaceFacePointLabel
517  << ' ' << blockPpoints[blockPfaceFacePointLabel]
518  << " of face "
519  << blockPfaceLabel
520  << " of block "
521  << blockPlabel
522  << exit(FatalError);
523  }
524  }
525  }
526 
527  forAll(blockNfaceFaces, blockNfaceFaceLabel)
528  {
529  const FixedList<label, 4>& blockNfaceFacePoints
530  = blockNfaceFaces[blockNfaceFaceLabel];
531 
532  forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
533  {
534  label NpointLabel =
535  blockNfaceFacePoints[blockNfaceFacePointLabel]
536  + blockOffsets_[blockNlabel];
537 
538  if (mergeList_[NpointLabel] == -1)
539  {
540  FatalErrorInFunction
541  << "unable to merge point "
542  << blockNfaceFacePointLabel
543  << ' ' << blockNpoints[blockNfaceFacePointLabel]
544  << " of face "
545  << blockNfaceLabel
546  << " of block "
547  << blockNlabel
548  << exit(FatalError);
549  }
550  }
551  }
552  }
553 
554 
555  // Sort merge list to return new point label (in new shorter list)
556  // given old point label
557  label newPointLabel = 0;
558 
559  forAll(mergeList_, pointLabel)
560  {
561  if (mergeList_[pointLabel] > pointLabel)
562  {
563  FatalErrorInFunction
564  << "Merge list contains point index out of range"
565  << exit(FatalError);
566  }
567 
568  if
569  (
570  mergeList_[pointLabel] == -1
571  || mergeList_[pointLabel] == pointLabel
572  )
573  {
574  mergeList_[pointLabel] = newPointLabel;
575  newPointLabel++;
576  }
577  else
578  {
579  mergeList_[pointLabel] = mergeList_[mergeList_[pointLabel]];
580  }
581  }
582 
583  nPoints_ = newPointLabel;
584 }
585 
586 
587 // ************************************************************************* //
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::List::resize
void resize(const label len)
Adjust allocated size of list.
Definition: ListI.H:139
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
blockMesh.H
Foam::blockMesh::topology
const polyMesh & topology() const
Definition: blockMeshCreate.C:294
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::flush
Ostream & flush(Ostream &os)
Flush stream.
Definition: Ostream.H:361
topoMesh
const polyMesh & topoMesh
Definition: blockMeshOBJ.H:29
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::cellList
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
Foam::List::subList
SubList< T > subList
Declare type of subList.
Definition: List.H:112
Foam::FatalError
error FatalError
Foam::labelListListList
List< labelListList > labelListListList
A List of labelListList.
Definition: labelList.H:57
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::blockMesh::blockList
PtrList< block > blockList
The list of blocks is stored as a PtrList.
Definition: blockMesh.H:174
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::cp
bool cp(const fileName &src, const fileName &dst, const bool followLink=true)
Copy the source to the destination (recursively if necessary).
Definition: MSwindows.C:802
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::blockMesh::points
const pointField & points() const
Definition: blockMesh.C:386
Foam::faceList
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47