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