blockMeshMergeFast.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) 2015-2021 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 Functions * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32 
33 // faces
34 // 6
35 // (
36 // 4(0 4 7 3) // x-min
37 // 4(1 2 6 5) // x-max
38 // 4(0 1 5 4) // y-min
39 // 4(3 7 6 2) // y-max
40 // 4(0 3 2 1) // z-min
41 // 4(4 5 6 7) // z-max
42 // );
43 
44 // Face-edge directions
45 static const int faceEdgeDirs[6][4] =
46 {
47  {2, 1, -2, -1},
48  {1, 2, -1, -2},
49  {1, 2, -1, -2},
50  {2, 1, -2, -1},
51  {2, 1, -2, -1},
52  {1, 2, -1, -2}
53 };
54 
55 // The face-face-rotation direction correspondence map
56 static Pair<int> faceFaceRotMap[6][6][4];
57 
58 // Generate the face-face-rotation direction correspondence map
60 {
61  for(int facePi=0; facePi<6; facePi++)
62  {
63  for(int faceNi=0; faceNi<6; faceNi++)
64  {
65  for(int rot=0; rot<4; rot++)
66  {
67  Pair<int>& map = faceFaceRotMap[facePi][faceNi][rot];
68 
69  for(int Pp=0; Pp<2; Pp++)
70  {
71  int Pdir = faceEdgeDirs[facePi][Pp];
72  int Np = (3 - Pp + rot)%4;
73  int Ndir = faceEdgeDirs[faceNi][Np];
74  map[Pdir-1] = -Ndir;
75  }
76 
77  // Handle sign change due to match-face transpose
78  if (mag(map[0]) == 2 && map[0]*map[1] < 0)
79  {
80  map[0] = -map[0];
81  map[1] = -map[1];
82  }
83  }
84  }
85  }
86 }
87 
88 // Return the direction map for the merge-faces
90 (
91  const label facePi,
92  const face& faceP,
93  const label faceNi,
94  const face& faceN
95 )
96 {
97  // Search for the point on faceN corresponding to the 0-point on faceP
98  for(int rot=0; rot<4; rot++)
99  {
100  if (faceN[rot] == faceP[0])
101  {
102  return faceFaceRotMap[facePi][faceNi][rot];
103  }
104  }
105 
107  << "Cannot find point correspondence for faces "
108  << faceP << " and " << faceN
109  << exit(FatalError);
110 
111  return Pair<int>(0, 0);
112 }
113 
114 // Set the block and face indices for all the merge faces
116 (
117  const cellList& topoCells,
118  const faceList::subList& topoInternalFaces,
119  const labelList& topoFaceCell,
120  List<Pair<label>>& mergeBlock
121 )
122 {
123  forAll(topoInternalFaces, topoFacei)
124  {
125  label topoPi = topoFaceCell[topoFacei];
126  const labelList& topoPfaces = topoCells[topoPi];
127 
128  bool foundFace = false;
129  label topoPfacei;
130  for
131  (
132  topoPfacei = 0;
133  topoPfacei < topoPfaces.size();
134  topoPfacei++
135  )
136  {
137  if (topoPfaces[topoPfacei] == topoFacei)
138  {
139  foundFace = true;
140  break;
141  }
142  }
143 
144  if (!foundFace)
145  {
147  << "Cannot find merge face for block " << topoPi
148  << exit(FatalError);
149  }
150 
151  mergeBlock[topoFacei].first() = topoPi;
152  mergeBlock[topoFacei].second() = topoPfacei;
153  }
154 }
155 
156 // Return the number of divisions in each direction for the face
157 Pair<label> faceNij(const label facei, const block& block)
158 {
159  Pair<label> fnij;
160 
161  int i = facei/2;
162 
163  if (i == 0)
164  {
165  fnij.first() = block.density().y() + 1;
166  fnij.second() = block.density().z() + 1;
167  }
168  else if (i == 1)
169  {
170  fnij.first() = block.density().x() + 1;
171  fnij.second() = block.density().z() + 1;
172  }
173  else if (i == 2)
174  {
175  fnij.first() = block.density().x() + 1;
176  fnij.second() = block.density().y() + 1;
177  }
178 
179  return fnij;
180 }
181 
182 // Sign the index corresponding to the map
183 inline label signIndex(const int map, const label i)
184 {
185  return map < 0 ? -i-1 : i;
186 }
187 
188 // Reverse a signed index with the number of divisions
189 inline label unsignIndex(const label i, const label ni)
190 {
191  return i >= 0 ? i : ni + i + 1;
192 }
193 
194 // Return the mapped index
195 inline label mapij(const int map, const label i, const label j)
196 {
197  return signIndex(map, mag(map) == 1 ? i : j);
198 }
199 
200 // Return the face point index
201 inline label facePoint
202 (
203  const int facei,
204  const block& block,
205  const label i,
206  const label j
207 )
208 {
209  switch (facei)
210  {
211  case 0:
212  return block.pointLabel(0, i, j);
213  case 1:
214  return block.pointLabel(block.density().x(), i, j);
215  case 2:
216  return block.pointLabel(i, 0, j);
217  case 3:
218  return block.pointLabel(i, block.density().y(), j);
219  case 4:
220  return block.pointLabel(i, j, 0);
221  case 5:
222  return block.pointLabel(i, j, block.density().z());
223  default:
224  return -1;
225  }
226 }
227 
228 // Return the neighbour face point from the signed indices
229 inline label facePointN
230 (
231  const block& block,
232  const label i,
233  const label j,
234  const label k
235 )
236 {
237  return block.pointLabel
238  (
239  unsignIndex(i, block.density().x()),
240  unsignIndex(j, block.density().y()),
241  unsignIndex(k, block.density().z())
242  );
243 }
244 
245 // Return the neighbour face point from the mapped indices
246 inline label facePointN
247 (
248  const int facei,
249  const block& block,
250  const label i,
251  const label j
252 )
253 {
254  switch (facei)
255  {
256  case 0:
257  return facePointN(block, 0, i, j);
258  case 1:
259  return facePointN(block, block.density().x(), i, j);
260  case 2:
261  return facePointN(block, i, 0, j);
262  case 3:
263  return facePointN(block, i, block.density().y(), j);
264  case 4:
265  return facePointN(block, i, j, 0);
266  case 5:
267  return facePointN(block, i, j, block.density().z());
268  default:
269  return -1;
270  }
271 }
272 
273 // Return the neighbour face point using the map
274 inline label facePointN
275 (
276  const int facei,
277  const Pair<int>& fmap,
278  const block& block,
279  const label i,
280  const label j
281 )
282 {
283  return facePointN(facei, block, mapij(fmap[0], i, j), mapij(fmap[1], i, j));
284 }
285 
286 }
287 
288 
289 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
290 
291 void Foam::blockMesh::calcMergeInfoFast()
292 {
293  // Generate the static face-face map
295 
296  const blockList& blocks = *this;
297 
298  if (verboseOutput)
299  {
300  Info<< "Creating block offsets" << endl;
301  }
302 
303  blockOffsets_.setSize(blocks.size());
304 
305  nPoints_ = 0;
306  nCells_ = 0;
307 
308  forAll(blocks, blocki)
309  {
310  blockOffsets_[blocki] = nPoints_;
311 
312  nPoints_ += blocks[blocki].nPoints();
313  nCells_ += blocks[blocki].nCells();
314  }
315 
316  if (verboseOutput)
317  {
318  Info<< "Creating merge list using the fast topological search"
319  << flush;
320  }
321 
322  // Size merge list and initialise to -1
323  mergeList_.setSize(nPoints_, -1);
324 
325  // Block mesh topology
326  const pointField& topoPoints = topology().points();
327  const cellList& topoCells = topology().cells();
328  const faceList& topoFaces = topology().faces();
329  const labelList& topoFaceOwn = topology().faceOwner();
330  const labelList& topoFaceNei = topology().faceNeighbour();
331 
332  // Topological merging is only necessary for internal block faces
333  // Note edge and face collapse may apply to boundary faces
334  // but is not yet supported in the "fast" algorithm
335  const faceList::subList topoInternalFaces
336  (
337  topoFaces,
338  topology().nInternalFaces()
339  );
340 
341  List<Pair<label>> mergeBlockP(topoInternalFaces.size());
343  (
344  topoCells,
345  topoInternalFaces,
346  topoFaceOwn,
347  mergeBlockP
348  );
349 
350  List<Pair<label>> mergeBlockN(topoInternalFaces.size());
352  (
353  topoCells,
354  topoInternalFaces,
355  topoFaceNei,
356  mergeBlockN
357  );
358 
359  if (debug)
360  {
361  Info<< endl;
362  }
363 
364  forAll(topoInternalFaces, topoFacei)
365  {
366  if (debug)
367  {
368  Info<< "Processing face " << topoFacei << endl;
369  }
370 
371  label blockPi = mergeBlockP[topoFacei].first();
372  label blockPfacei = mergeBlockP[topoFacei].second();
373 
374  label blockNi = mergeBlockN[topoFacei].first();
375  label blockNfacei = mergeBlockN[topoFacei].second();
376 
377  Pair<int> fmap
378  (
379  faceMap
380  (
381  blockPfacei,
382  blocks[blockPi].blockShape().faces()[blockPfacei],
383  blockNfacei,
384  blocks[blockNi].blockShape().faces()[blockNfacei]
385  )
386  );
387 
388  if (debug)
389  {
390  Info<< " Face map for faces "
391  << blocks[blockPi].blockShape().faces()[blockPfacei] << " "
392  << blocks[blockNi].blockShape().faces()[blockNfacei] << ": "
393  << fmap << endl;
394  }
395 
396  const pointField& blockPpoints = blocks[blockPi].points();
397  const pointField& blockNpoints = blocks[blockNi].points();
398 
399  Pair<label> Pnij(faceNij(blockPfacei, blocks[blockPi]));
400 
401  // Check block subdivision correspondence
402  {
403  Pair<label> Nnij(faceNij(blockNfacei, blocks[blockNi]));
404  Pair<label> NPnij;
405  NPnij[0] = Nnij[mag(fmap[0]) - 1];
406  NPnij[1] = Nnij[mag(fmap[1]) - 1];
407 
408  if (Pnij != NPnij)
409  {
411  << "Sub-division mismatch between face "
412  << blockPfacei << " of block " << blockPi << Pnij
413  << " and face "
414  << blockNfacei << " of block " << blockNi << Nnij
415  << exit(FatalError);
416  }
417  }
418 
419  // Calculate a suitable test distance from the bounding box of the face.
420  // Note this is used only as a sanity check and for diagnostics in
421  // case there is a grading inconsistency.
422  const boundBox bb(topoCells[blockPi].points(topoFaces, topoPoints));
423  const scalar testSqrDist = magSqr(1e-6*bb.span());
424 
425  // Accumulate the maximum merge distance for diagnostics
426  scalar maxSqrDist = 0;
427 
428  for (label j=0; j<Pnij.second(); j++)
429  {
430  for (label i=0; i<Pnij.first(); i++)
431  {
432  label blockPpointi =
433  facePoint(blockPfacei, blocks[blockPi], i, j);
434 
435  label blockNpointi =
436  facePointN(blockNfacei, fmap, blocks[blockNi], i, j);
437 
438  scalar sqrDist
439  (
440  magSqr
441  (
442  blockPpoints[blockPpointi]
443  - blockNpoints[blockNpointi]
444  )
445  );
446 
447  if (sqrDist > testSqrDist)
448  {
450  << "Point merge failure between face "
451  << blockPfacei << " of block " << blockPi
452  << " and face "
453  << blockNfacei << " of block " << blockNi
454  << endl
455  << " Points: " << blockPpoints[blockPpointi]
456  << " " << blockNpoints[blockNpointi]
457  << endl
458  << " This may be due to inconsistent grading."
459  << exit(FatalError);
460  }
461 
462  maxSqrDist = max(maxSqrDist, sqrDist);
463 
464  label Ppointi = blockPpointi + blockOffsets_[blockPi];
465  label Npointi = blockNpointi + blockOffsets_[blockNi];
466 
467  label minPNi = min(Ppointi, Npointi);
468 
469  if (mergeList_[Ppointi] != -1)
470  {
471  minPNi = min(minPNi, mergeList_[Ppointi]);
472  }
473 
474  if (mergeList_[Npointi] != -1)
475  {
476  minPNi = min(minPNi, mergeList_[Npointi]);
477  }
478 
479  mergeList_[Ppointi] = mergeList_[Npointi] = minPNi;
480  }
481  }
482 
483  if (debug)
484  {
485  Info<< " Max distance between merge points: "
486  << sqrt(maxSqrDist) << endl;
487  }
488  }
489 
490 
491  bool changedPointMerge = false;
492  label nPasses = 0;
493 
494  do
495  {
496  changedPointMerge = false;
497  nPasses++;
498 
499  forAll(topoInternalFaces, topoFacei)
500  {
501  label blockPi = mergeBlockP[topoFacei].first();
502  label blockPfacei = mergeBlockP[topoFacei].second();
503 
504  label blockNi = mergeBlockN[topoFacei].first();
505  label blockNfacei = mergeBlockN[topoFacei].second();
506 
507  Pair<int> fmap
508  (
509  faceMap
510  (
511  blockPfacei,
512  blocks[blockPi].blockShape().faces()[blockPfacei],
513  blockNfacei,
514  blocks[blockNi].blockShape().faces()[blockNfacei]
515  )
516  );
517 
518  Pair<label> Pnij(faceNij(blockPfacei, blocks[blockPi]));
519 
520  for (label j=0; j<Pnij.second(); j++)
521  {
522  for (label i=0; i<Pnij.first(); i++)
523  {
524  label blockPpointi =
525  facePoint(blockPfacei, blocks[blockPi], i, j);
526 
527  label blockNpointi =
528  facePointN(blockNfacei, fmap, blocks[blockNi], i, j);
529 
530  label Ppointi =
531  blockPpointi + blockOffsets_[blockPi];
532 
533  label Npointi =
534  blockNpointi + blockOffsets_[blockNi];
535 
536  if (mergeList_[Ppointi] != mergeList_[Npointi])
537  {
538  changedPointMerge = true;
539 
540  mergeList_[Ppointi]
541  = mergeList_[Npointi]
542  = min(mergeList_[Ppointi], mergeList_[Npointi]);
543  }
544  }
545  }
546  }
547 
548  if (verboseOutput)
549  {
550  Info<< "." << flush;
551  }
552 
553  if (nPasses > 100)
554  {
556  << "Point merging failed after 100 passes."
557  << exit(FatalError);
558  }
559 
560  } while (changedPointMerge);
561 
562 
563  // Sort merge list and count number of unique points
564  label nUniqPoints = 0;
565 
566  forAll(mergeList_, pointi)
567  {
568  if (mergeList_[pointi] > pointi)
569  {
571  << "Merge list contains point index out of range"
572  << exit(FatalError);
573  }
574 
575  if (mergeList_[pointi] == -1 || mergeList_[pointi] == pointi)
576  {
577  mergeList_[pointi] = nUniqPoints++;
578  }
579  else
580  {
581  mergeList_[pointi] = mergeList_[mergeList_[pointi]];
582  }
583  }
584 
585  // Correct number of points in mesh
586  nPoints_ = nUniqPoints;
587 }
588 
589 
590 // ************************************************************************* //
void genFaceFaceRotMap()
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
label facePointN(const block &block, const label i, const label j, const label k)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Vector< label > & density() const
Return the mesh density (number of cells) in the i,j,k directions.
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
label unsignIndex(const label i, const label ni)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Holds information (coordinate and normal) regarding nearest wall point.
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
A List obtained as a section of another List.
Definition: SubList.H:53
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
const pointField & points
label signIndex(const int map, const label i)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Pair< label > faceNij(const label facei, const block &block)
label mapij(const int map, const label i, const label j)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
label facePoint(const int facei, const block &block, const label i, const label j)
const Type & second() const
Return second.
Definition: Pair.H:110
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
Ostream & flush(Ostream &os)
Flush stream.
Definition: Ostream.H:243
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:84
label pointLabel(const label i, const label j, const label k) const
Vertex label offset for a particular i,j,k position.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
const Type & first() const
Return first.
Definition: Pair.H:98
void setBlockFaceCorrespondence(const cellList &topoCells, const faceList::subList &topoInternalFaces, const labelList &topoFaceCell, List< Pair< label >> &mergeBlock)
Namespace for OpenFOAM.