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