FreeFOAM The Cross-Platform CFD Toolkit
createMergeList.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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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 Foam::labelList Foam::blockMesh::createMergeList()
31 {
32  Info<< nl << "Creating merge list " << flush;
33 
34  labelList MergeList(nPoints_, -1);
35 
36  blockMesh& blocks = *this;
37 
38  const pointField& blockPoints = topology().points();
39  const cellList& blockCells = topology().cells();
40  const faceList& blockFaces = topology().faces();
41  const labelList& faceOwnerBlocks = topology().faceOwner();
42 
43  // For efficiency, create merge pairs in the first pass
44  labelListListList glueMergePairs(blockFaces.size());
45 
46  const labelList& faceNeighbourBlocks = topology().faceNeighbour();
47 
48  forAll(blockFaces, blockFaceLabel)
49  {
50  label blockPlabel = faceOwnerBlocks[blockFaceLabel];
51  const pointField& blockPpoints = blocks[blockPlabel].points();
52  const labelList& blockPfaces = blockCells[blockPlabel];
53 
54  bool foundFace = false;
55  label blockPfaceLabel;
56  for
57  (
58  blockPfaceLabel = 0;
59  blockPfaceLabel < blockPfaces.size();
60  blockPfaceLabel++
61  )
62  {
63  if
64  (
65  blockFaces[blockPfaces[blockPfaceLabel]]
66  == blockFaces[blockFaceLabel]
67  )
68  {
69  foundFace = true;
70  break;
71  }
72  }
73 
74  if (!foundFace)
75  {
76  FatalErrorIn("blockMesh::createMergeList()")
77  << "Cannot find merge face for block " << blockPlabel
78  << exit(FatalError);
79  };
80 
81  const labelListList& blockPfaceFaces =
82  blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
83 
84  labelListList& curPairs = glueMergePairs[blockFaceLabel];
85  curPairs.setSize(blockPfaceFaces.size());
86 
87  // Calculate sqr of the merge tolerance as 1/10th of the min sqr
88  // point to point distance on the block face.
89  // At the same time merge collated points on the block's faces
90  // (removes boundary poles etc.)
91  // Co-located points detected by initally taking a constant factor of
92  // the size of the block face:
93  boundBox bb(blockFaces[blockFaceLabel].points(blockPoints));
94  const scalar mergeSqrDist =
95  SMALL*magSqr(bb.span())/blockPfaceFaces.size();
96  scalar sqrMergeTol = GREAT;
97 
98  // This is an N^2 algorithm
99  forAll(blockPfaceFaces, blockPfaceFaceLabel)
100  {
101  const labelList& blockPfaceFacePoints
102  = blockPfaceFaces[blockPfaceFaceLabel];
103 
104  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
105  {
106  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel2)
107  {
108  if (blockPfaceFacePointLabel != blockPfaceFacePointLabel2)
109  {
110  scalar magSqrDist = magSqr
111  (
112  blockPpoints[blockPfaceFacePoints
113  [blockPfaceFacePointLabel]]
114  - blockPpoints[blockPfaceFacePoints
115  [blockPfaceFacePointLabel2]]
116  );
117 
118  if (magSqrDist < mergeSqrDist)
119  {
120  label PpointLabel =
121  blockPfaceFacePoints[blockPfaceFacePointLabel]
122  + blockOffsets_[blockPlabel];
123 
124  label PpointLabel2 =
125  blockPfaceFacePoints[blockPfaceFacePointLabel2]
126  + blockOffsets_[blockPlabel];
127 
128  label minPP2 = min(PpointLabel, PpointLabel2);
129 
130  if (MergeList[PpointLabel] != -1)
131  {
132  minPP2 = min(minPP2, MergeList[PpointLabel]);
133  }
134 
135  if (MergeList[PpointLabel2] != -1)
136  {
137  minPP2 = min(minPP2, MergeList[PpointLabel2]);
138  }
139 
140  MergeList[PpointLabel] = MergeList[PpointLabel2]
141  = minPP2;
142  }
143  else
144  {
145  sqrMergeTol = min(sqrMergeTol, magSqrDist);
146  }
147  }
148  }
149  }
150  }
151 
152  sqrMergeTol /= 10.0;
153 
154 
155  if (topology().isInternalFace(blockFaceLabel))
156  {
157  label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
158  const pointField& blockNpoints = blocks[blockNlabel].points();
159  const labelList& blockNfaces = blockCells[blockNlabel];
160 
161  foundFace = false;
162  label blockNfaceLabel;
163  for
164  (
165  blockNfaceLabel = 0;
166  blockNfaceLabel < blockNfaces.size();
167  blockNfaceLabel++
168  )
169  {
170  if
171  (
172  blockFaces[blockNfaces[blockNfaceLabel]]
173  == blockFaces[blockFaceLabel]
174  )
175  {
176  foundFace = true;
177  break;
178  }
179  }
180 
181  if (!foundFace)
182  {
183  FatalErrorIn("blockMesh::createMergeList()")
184  << "Cannot find merge face for block " << blockNlabel
185  << exit(FatalError);
186  };
187 
188  const labelListList& blockNfaceFaces =
189  blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
190 
191  if (blockPfaceFaces.size() != blockNfaceFaces.size())
192  {
193  FatalErrorIn("blockMesh::createMergeList()")
194  << "Inconsistent number of faces between block pair "
195  << blockPlabel << " and " << blockNlabel
196  << exit(FatalError);
197  }
198 
199 
200  bool found = false;
201 
202  // N-squared point search over all points of all faces of
203  // master block over all point of all faces of slave block
204  forAll(blockPfaceFaces, blockPfaceFaceLabel)
205  {
206  const labelList& blockPfaceFacePoints =
207  blockPfaceFaces[blockPfaceFaceLabel];
208 
209  labelList& cp = curPairs[blockPfaceFaceLabel];
210  cp.setSize(blockPfaceFacePoints.size());
211  cp = -1;
212 
213  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
214  {
215  found = false;
216 
217  forAll(blockNfaceFaces, blockNfaceFaceLabel)
218  {
219  const labelList& blockNfaceFacePoints
220  = blockNfaceFaces[blockNfaceFaceLabel];
221 
222  forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
223  {
224  if
225  (
226  magSqr
227  (
228  blockPpoints
229  [
230  blockPfaceFacePoints
231  [blockPfaceFacePointLabel]
232  ]
233  - blockNpoints
234  [
235  blockNfaceFacePoints
236  [blockNfaceFacePointLabel]
237  ]
238  ) < sqrMergeTol
239  )
240  {
241  // Found a new pair
242  found = true;
243 
244  cp[blockPfaceFacePointLabel] =
245  blockNfaceFacePoints[blockNfaceFacePointLabel];
246 
247  label PpointLabel =
248  blockPfaceFacePoints
249  [
250  blockPfaceFacePointLabel
251  ]
252  + blockOffsets_[blockPlabel];
253 
254  label NpointLabel =
255  blockNfaceFacePoints
256  [
257  blockNfaceFacePointLabel
258  ]
259  + blockOffsets_[blockNlabel];
260 
261  label minPN = min(PpointLabel, NpointLabel);
262 
263  if (MergeList[PpointLabel] != -1)
264  {
265  minPN = min(minPN, MergeList[PpointLabel]);
266  }
267 
268  if (MergeList[NpointLabel] != -1)
269  {
270  minPN = min(minPN, MergeList[NpointLabel]);
271  }
272 
273  MergeList[PpointLabel] =
274  MergeList[NpointLabel] =
275  minPN;
276  }
277  }
278  }
279  }
280  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
281  {
282  if (cp[blockPfaceFacePointLabel] == -1)
283  {
284  FatalErrorIn("blockMesh::createMergeList()")
285  << "Inconsistent point locations between blocks "
286  << blockPlabel << " and " << blockNlabel << nl
287  << " probably due to inconsistent grading."
288  << exit(FatalError);
289  }
290  }
291  }
292  }
293  }
294 
295 
296  const faceList::subList blockInternalFaces
297  (
298  blockFaces,
299  topology().nInternalFaces()
300  );
301 
302  bool changedPointMerge = false;
303  label nPasses = 0;
304 
305  do
306  {
307  changedPointMerge = false;
308  nPasses++;
309 
310  forAll(blockInternalFaces, blockFaceLabel)
311  {
312  label blockPlabel = faceOwnerBlocks[blockFaceLabel];
313  label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
314 
315  const labelList& blockPfaces = blockCells[blockPlabel];
316  const labelList& blockNfaces = blockCells[blockNlabel];
317 
318  const labelListList& curPairs = glueMergePairs[blockFaceLabel];
319 
320  bool foundFace = false;
321  label blockPfaceLabel;
322  for
323  (
324  blockPfaceLabel = 0;
325  blockPfaceLabel < blockPfaces.size();
326  blockPfaceLabel++
327  )
328  {
329  if
330  (
331  blockFaces[blockPfaces[blockPfaceLabel]]
332  == blockInternalFaces[blockFaceLabel]
333  )
334  {
335  foundFace = true;
336  break;
337  }
338  }
339 
340  foundFace = false;
341  label blockNfaceLabel;
342  for
343  (
344  blockNfaceLabel = 0;
345  blockNfaceLabel < blockNfaces.size();
346  blockNfaceLabel++
347  )
348  {
349  if
350  (
351  blockFaces[blockNfaces[blockNfaceLabel]]
352  == blockInternalFaces[blockFaceLabel]
353  )
354  {
355  foundFace = true;
356  break;
357  }
358  }
359 
360  const labelListList& blockPfaceFaces =
361  blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
362 
363  forAll(blockPfaceFaces, blockPfaceFaceLabel)
364  {
365  const labelList& blockPfaceFacePoints
366  = blockPfaceFaces[blockPfaceFaceLabel];
367 
368  const labelList& cp = curPairs[blockPfaceFaceLabel];
369 
370  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
371  {
372  label PpointLabel =
373  blockPfaceFacePoints[blockPfaceFacePointLabel]
374  + blockOffsets_[blockPlabel];
375 
376  label NpointLabel =
377  cp[blockPfaceFacePointLabel]
378  + blockOffsets_[blockNlabel];
379 
380  if
381  (
382  MergeList[PpointLabel]
383  != MergeList[NpointLabel]
384  )
385  {
386  changedPointMerge = true;
387 
388  MergeList[PpointLabel]
389  = MergeList[NpointLabel]
390  = min
391  (
392  MergeList[PpointLabel],
393  MergeList[NpointLabel]
394  );
395  }
396  }
397  }
398  }
399  Info << "." << flush;
400 
401  if (nPasses > 100)
402  {
403  FatalErrorIn("blockMesh::createMergeList()")
404  << "Point merging failed after max number of passes."
405  << abort(FatalError);
406  }
407  }
408  while (changedPointMerge);
409  Info << endl;
410 
411  forAll(blockInternalFaces, blockFaceLabel)
412  {
413  label blockPlabel = faceOwnerBlocks[blockFaceLabel];
414  label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
415 
416  const labelList& blockPfaces = blockCells[blockPlabel];
417  const labelList& blockNfaces = blockCells[blockNlabel];
418 
419  const pointField& blockPpoints = blocks[blockPlabel].points();
420  const pointField& blockNpoints = blocks[blockNlabel].points();
421 
422  bool foundFace = false;
423  label blockPfaceLabel;
424  for
425  (
426  blockPfaceLabel = 0;
427  blockPfaceLabel < blockPfaces.size();
428  blockPfaceLabel++
429  )
430  {
431  if
432  (
433  blockFaces[blockPfaces[blockPfaceLabel]]
434  == blockInternalFaces[blockFaceLabel]
435  )
436  {
437  foundFace = true;
438  break;
439  }
440  }
441 
442  if (!foundFace)
443  {
444  FatalErrorIn("blockMesh::createMergeList()")
445  << "Cannot find merge face for block " << blockPlabel
446  << exit(FatalError);
447  };
448 
449  foundFace = false;
450  label blockNfaceLabel;
451  for
452  (
453  blockNfaceLabel = 0;
454  blockNfaceLabel < blockNfaces.size();
455  blockNfaceLabel++
456  )
457  {
458  if
459  (
460  blockFaces[blockNfaces[blockNfaceLabel]]
461  == blockInternalFaces[blockFaceLabel]
462  )
463  {
464  foundFace = true;
465  break;
466  }
467  }
468 
469  if (!foundFace)
470  {
471  FatalErrorIn("blockMesh::createMergeList()")
472  << "Cannot find merge face for block " << blockNlabel
473  << exit(FatalError);
474  };
475 
476  const labelListList& blockPfaceFaces =
477  blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
478 
479  const labelListList& blockNfaceFaces =
480  blocks[blockNlabel].boundaryPatches()[blockNfaceLabel];
481 
482  forAll(blockPfaceFaces, blockPfaceFaceLabel)
483  {
484  const labelList& blockPfaceFacePoints
485  = blockPfaceFaces[blockPfaceFaceLabel];
486 
487  forAll(blockPfaceFacePoints, blockPfaceFacePointLabel)
488  {
489  label PpointLabel =
490  blockPfaceFacePoints[blockPfaceFacePointLabel]
491  + blockOffsets_[blockPlabel];
492 
493  if (MergeList[PpointLabel] == -1)
494  {
495  FatalErrorIn("blockMesh::createMergeList()")
496  << "Unable to merge point "
497  << blockPfaceFacePointLabel
498  << ' ' << blockPpoints[blockPfaceFacePointLabel]
499  << " of face "
500  << blockPfaceLabel
501  << " of block "
502  << blockPlabel
503  << exit(FatalError);
504  }
505  }
506  }
507 
508  forAll(blockNfaceFaces, blockNfaceFaceLabel)
509  {
510  const labelList& blockNfaceFacePoints
511  = blockNfaceFaces[blockNfaceFaceLabel];
512 
513  forAll(blockNfaceFacePoints, blockNfaceFacePointLabel)
514  {
515  label NpointLabel =
516  blockNfaceFacePoints[blockNfaceFacePointLabel]
517  + blockOffsets_[blockNlabel];
518 
519  if (MergeList[NpointLabel] == -1)
520  {
521  FatalErrorIn("blockMesh::createMergeList()")
522  << "unable to merge point "
523  << blockNfaceFacePointLabel
524  << ' ' << blockNpoints[blockNfaceFacePointLabel]
525  << " of face "
526  << blockNfaceLabel
527  << " of block "
528  << blockNlabel
529  << exit(FatalError);
530  }
531  }
532  }
533  }
534 
535 
536  // sort merge list to return new point label (in new shorter list)
537  // given old point label
538  label newPointLabel = 0;
539 
540  forAll(MergeList, pointLabel)
541  {
542  if (MergeList[pointLabel] > pointLabel)
543  {
544  FatalErrorIn("blockMesh::createMergeList()")
545  << "ouch" << exit(FatalError);
546  }
547 
548  if
549  (
550  (MergeList[pointLabel] == -1)
551  || MergeList[pointLabel] == pointLabel
552  )
553  {
554  MergeList[pointLabel] = newPointLabel;
555  newPointLabel++;
556  }
557  else
558  {
559  MergeList[pointLabel] = MergeList[MergeList[pointLabel]];
560  }
561  }
562 
563  nPoints_ = newPointLabel;
564 
565 
566  return MergeList;
567 }
568 
569 // ************************ vim: set sw=4 sts=4 et: ************************ //