createPolyCells.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-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 Description
25  create cellPolys
26  - use pointCells when searching for connectivity
27  - initialise the cell connectivity with '-1'
28  - find both cell faces corresponding to the baffles and mark them
29  to prevent a connection
30  - standard connectivity checks
31 
32  - added baffle support
33 
34 \*---------------------------------------------------------------------------*/
35 
36 #include "meshReader.H"
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 void Foam::meshReader::createPolyCells()
41 {
42  // loop through all cell faces and create connectivity. This will produce
43  // a global face list and will describe all cells as lists of face labels
44 
45  const faceListList& cFaces = cellFaces();
46 
47  // count the maximum number of faces and set the size of the cellPolys_
48  cellPolys_.setSize(cFaces.size());
49 
50  label maxFaces = 0;
51 
52  forAll(cellPolys_, celli)
53  {
54  cellPolys_[celli].setSize(cFaces[celli].size(), -1);
55 
56  maxFaces += cFaces[celli].size();
57  }
58 
59  Info<< "Maximum possible number of faces in mesh: " << maxFaces << endl;
60 
61  meshFaces_.setSize(maxFaces);
62 
63  // set reference to point-cell addressing
64  const labelListList& ptCells = pointCells();
65 
66  // size the baffle lists and initialise to -1
67  baffleIds_.setSize(baffleFaces_.size());
68  forAll(baffleIds_, baffleI)
69  {
70  baffleIds_[baffleI].setSize(2);
71  }
72 
73  // block off baffles first
74  //
75  // To prevent internal faces, we'll mark the cell faces
76  // with negative cell ids (offset by nCells).
77  // eg,
78  // celli = -(nCells + baffleI)
79  //
80  // To distinguish these from the normal '-1' marker, we require
81  // celli = -(nCells + baffleI) < -1
82  //
83  // This condition is met provided that nCells > 1.
84  // ie., baffles require at least 2 volume cells
85 
86  label baffleOffset = cFaces.size();
87  forAll(baffleFaces_, baffleI)
88  {
89  label celli = -(baffleOffset + baffleI);
90  const face& curFace = baffleFaces_[baffleI];
91 
92  // get the list of labels
93  const labelList& curPoints = curFace;
94 
95  // a baffle is a single face - only need to match one face
96  // get the list of cells sharing this point
97  const labelList& curNeighbours = ptCells[curPoints[0]];
98 
99  label nNeighbours = 0;
100 
101  // For all neighbours
102  forAll(curNeighbours, neiI)
103  {
104  label curNei = curNeighbours[neiI];
105 
106  // get the list of search faces
107  const faceList& searchFaces = cFaces[curNei];
108 
109  forAll(searchFaces, neiFacei)
110  {
111  int cmp = face::compare(curFace, searchFaces[neiFacei]);
112 
113  if (cmp)
114  {
115  // maintain baffle orientation
116  // side0: baffle normal same as attached face
117  // side1: baffle normal opposite from attached face
118  //
119  label side = 0;
120  if (cmp < 0)
121  {
122  side = 1;
123  }
124 
125 #ifdef DEBUG_FACE_ORDERING
126  Info<< "cmp " << cmp << " matched " << curFace
127  << " with " << searchFaces[neiFacei]
128  << endl;
129 
130 
131  Info<< "match " << baffleI
132  << " (" << origCellId_[baffleOffset+baffleI] << ")"
133  << " side " << side
134  << " against cell " << curNei
135  << " face " << neiFacei
136  << " curFace " << curFace[1]
137  << " neiFace " << searchFaces[neiFacei][1]
138  << endl;
139 #endif
140 
141  if (baffleIds_[baffleI][side].unused())
142  {
143  baffleIds_[baffleI][side] = cellFaceIdentifier
144  (
145  curNei,
146  neiFacei
147  );
148 
149  nNeighbours++;
150  }
151  else
152  {
153  Info<< "multiple matches for side " << side
154  << " of baffle " << baffleI
155  << " (original cell "
156  << origCellId_[baffleOffset+baffleI] << ")"
157  << endl;
158  }
159  break;
160  }
161  }
162  if (nNeighbours >= 2) break;
163  }
164 
165  if (nNeighbours == 2)
166  {
167  for (label side = 0; side < nNeighbours; ++side)
168  {
169  label neiCell = baffleIds_[baffleI][side].cell;
170  label neiFace = baffleIds_[baffleI][side].face;
171 
172  if (baffleIds_[baffleI][side].used())
173  {
174  cellPolys_[neiCell][neiFace] = celli;
175  }
176  }
177  }
178  else
179  {
180  Info<< "drop baffle " << baffleI
181  << " (original cell "
182  << origCellId_[baffleOffset+baffleI] << ")"
183  << " with " << nNeighbours << " neighbours" << endl;
184 
185  baffleFaces_[baffleI].clear();
186  baffleIds_[baffleI].clear();
187  }
188  }
189 
190 #ifdef DEBUG_CELLPOLY
191  Info<< "cellPolys_" << cellPolys_ << endl;
192  Info<< "baffleFaces_" << baffleFaces_ << endl;
193  Info<< "baffleIds_" << baffleIds_ << endl;
194 #endif
195 
196  bool found = false;
197 
198  nInternalFaces_ = 0;
199 
200  forAll(cFaces, celli)
201  {
202  // Note:
203  // Insertion cannot be done in one go as the faces need to be
204  // added into the list in the increasing order of neighbour
205  // cells. Therefore, all neighbours will be detected first
206  // and then added in the correct order.
207 
208  const faceList& curFaces = cFaces[celli];
209 
210  // Record the neighbour cell
211  labelList neiCells(curFaces.size(), -1);
212 
213  // Record the face of neighbour cell
214  labelList faceOfNeiCell(curFaces.size(), -1);
215 
216  label nNeighbours = 0;
217 
218  // For all faces ...
219  forAll(curFaces, facei)
220  {
221  // Skip already matched faces or those tagged by baffles
222  if (cellPolys_[celli][facei] != -1) continue;
223 
224  found = false;
225 
226  const face& curFace = curFaces[facei];
227 
228  // get the list of labels
229  const labelList& curPoints = curFace;
230 
231  // For all points
232  forAll(curPoints, pointi)
233  {
234  // get the list of cells sharing this point
235  const labelList& curNeighbours = ptCells[curPoints[pointi]];
236 
237  // For all neighbours
238  forAll(curNeighbours, neiI)
239  {
240  label curNei = curNeighbours[neiI];
241 
242  // reject neighbours with the lower label. This should
243  // also reject current cell.
244  if (curNei > celli)
245  {
246  // get the list of search faces
247  const faceList& searchFaces = cFaces[curNei];
248 
249  forAll(searchFaces, neiFacei)
250  {
251  if (searchFaces[neiFacei] == curFace)
252  {
253  // Record the neighbour cell and face
254  neiCells[facei] = curNei;
255  faceOfNeiCell[facei] = neiFacei;
256  nNeighbours++;
257 #ifdef DEBUG_FACE_ORDERING
258  Info<< " cell " << celli
259  << " face " << facei
260  << " point " << pointi
261  << " nei " << curNei
262  << " neiFace " << neiFacei
263  << endl;
264 #endif
265  found = true;
266  break;
267  }
268  }
269  if (found) break;
270  }
271  if (found) break;
272  }
273  if (found) break;
274  } // End of current points
275  } // End of current faces
276 
277  // Add the faces in the increasing order of neighbours
278  for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
279  {
280  // Find the lowest neighbour which is still valid
281  label nextNei = -1;
282  label minNei = cellPolys_.size();
283 
284  forAll(neiCells, ncI)
285  {
286  if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
287  {
288  nextNei = ncI;
289  minNei = neiCells[ncI];
290  }
291  }
292 
293  if (nextNei > -1)
294  {
295  // Add the face to the list of faces
296  meshFaces_[nInternalFaces_] = curFaces[nextNei];
297 
298  // Mark for owner
299  cellPolys_[celli][nextNei] = nInternalFaces_;
300 
301  // Mark for neighbour
302  cellPolys_[neiCells[nextNei]][faceOfNeiCell[nextNei]] =
303  nInternalFaces_;
304 
305  // Stop the neighbour from being used again
306  neiCells[nextNei] = -1;
307 
308  // Increment number of faces counter
309  nInternalFaces_++;
310  }
311  else
312  {
314  << "Error in internal face insertion"
315  << abort(FatalError);
316  }
317  }
318  }
319 
320 #ifdef DEBUG_CELLPOLY
321  Info<< "cellPolys = " << cellPolys_ << endl;
322 #endif
323 
324  // don't reset the size of internal faces, because more faces will be
325  // added in createPolyBoundary()
326 }
327 
328 
329 // ************************************************************************* //
bool found
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
void setSize(const label)
Reset size of List.
Definition: List.C:281
static int compare(const face &, const face &)
Compare faces.
Definition: face.C:48
labelList origCellId_
Lookup original Cell number for a given cell.
Definition: meshReader.H:252
faceList baffleFaces_
List of each baffle face.
Definition: meshReader.H:271
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
errorManip< error > abort(error &err)
Definition: errorManip.H:131
messageStream Info
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
List< faceList > faceListList
Definition: faceListFwd.H:43
error FatalError
List< face > faceList
Definition: faceListFwd.H:41