meshStructure.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) 2013-2018 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 "meshStructure.H"
27 #include "FaceCellWave.H"
28 #include "topoDistanceData.H"
29 #include "pointTopoDistanceData.H"
30 #include "PointEdgeWave.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 defineTypeNameAndDebug(meshStructure, 0);
37 }
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 bool Foam::meshStructure::isStructuredCell
43 (
44  const polyMesh& mesh,
45  const label layerI,
46  const label celli
47 ) const
48 {
49  const cell& cFaces = mesh.cells()[celli];
50 
51  // Count number of side faces
52  label nSide = 0;
53  forAll(cFaces, i)
54  {
55  if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
56  {
57  nSide++;
58  }
59  }
60 
61  if (nSide != cFaces.size()-2)
62  {
63  return false;
64  }
65 
66  // Check that side faces have correct point layers
67  forAll(cFaces, i)
68  {
69  if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
70  {
71  const face& f = mesh.faces()[cFaces[i]];
72 
73  label nLayer = 0;
74  label nLayerPlus1 = 0;
75  forAll(f, fp)
76  {
77  label pointi = f[fp];
78  if (pointLayer_[pointi] == layerI)
79  {
80  nLayer++;
81  }
82  else if (pointLayer_[pointi] == layerI+1)
83  {
84  nLayerPlus1++;
85  }
86  }
87 
88  if (f.size() != 4 || (nLayer+nLayerPlus1 != 4))
89  {
90  return false;
91  }
92  }
93  }
94 
95  return true;
96 }
97 
98 
99 void Foam::meshStructure::correct
100 (
101  const polyMesh& mesh,
102  const uindirectPrimitivePatch& pp
103 )
104 {
105  // Field on cells and faces.
106  List<topoDistanceData> cellData(mesh.nCells());
107  List<topoDistanceData> faceData(mesh.nFaces());
108 
109  {
110  if (debug)
111  {
112  Info<< typeName << " : seeding "
113  << returnReduce(pp.size(), sumOp<label>()) << " patch faces"
114  << nl << endl;
115  }
116 
117 
118  // Start of changes
119  labelList patchFaces(pp.size());
120  List<topoDistanceData> patchData(pp.size());
121  forAll(pp, patchFacei)
122  {
123  patchFaces[patchFacei] = pp.addressing()[patchFacei];
124  patchData[patchFacei] = topoDistanceData(patchFacei, 0);
125  }
126 
127 
128  // Propagate information inwards
129  FaceCellWave<topoDistanceData> distanceCalc
130  (
131  mesh,
132  patchFaces,
133  patchData,
134  faceData,
135  cellData,
136  mesh.globalData().nTotalCells()+1
137  );
138 
139 
140  // Determine cells from face-cell-walk
141  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
142 
143  cellToPatchFaceAddressing_.setSize(mesh.nCells());
144  cellLayer_.setSize(mesh.nCells());
145  forAll(cellToPatchFaceAddressing_, celli)
146  {
147  cellToPatchFaceAddressing_[celli] = cellData[celli].data();
148  cellLayer_[celli] = cellData[celli].distance();
149  }
150 
151 
152 
153  // Determine faces from face-cell-walk
154  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
155 
156  faceToPatchFaceAddressing_.setSize(mesh.nFaces());
157  faceToPatchEdgeAddressing_.setSize(mesh.nFaces());
158  faceToPatchEdgeAddressing_ = labelMin;
159  faceLayer_.setSize(mesh.nFaces());
160 
161  forAll(faceToPatchFaceAddressing_, facei)
162  {
163  label own = mesh.faceOwner()[facei];
164  label patchFacei = faceData[facei].data();
165  label patchDist = faceData[facei].distance();
166 
167  if (mesh.isInternalFace(facei))
168  {
169  label nei = mesh.faceNeighbour()[facei];
170 
171  if (cellData[own].distance() == cellData[nei].distance())
172  {
173  // side face
174  faceToPatchFaceAddressing_[facei] = 0;
175  faceLayer_[facei] = cellData[own].distance();
176  }
177  else if (cellData[own].distance() < cellData[nei].distance())
178  {
179  // unturned face
180  faceToPatchFaceAddressing_[facei] = patchFacei+1;
181  faceToPatchEdgeAddressing_[facei] = -1;
182  faceLayer_[facei] = patchDist;
183  }
184  else
185  {
186  // turned face
187  faceToPatchFaceAddressing_[facei] = -(patchFacei+1);
188  faceToPatchEdgeAddressing_[facei] = -1;
189  faceLayer_[facei] = patchDist;
190  }
191  }
192  else if (patchDist == cellData[own].distance())
193  {
194  // starting face
195  faceToPatchFaceAddressing_[facei] = -(patchFacei+1);
196  faceToPatchEdgeAddressing_[facei] = -1;
197  faceLayer_[facei] = patchDist;
198  }
199  else
200  {
201  // unturned face or side face. Cannot be determined until
202  // we determine the point layers. Problem is that both are
203  // the same number of steps away from the initial seed face.
204  }
205  }
206  }
207 
208 
209  // Determine points from separate walk on point-edge
210  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211 
212  {
213  pointToPatchPointAddressing_.setSize(mesh.nPoints());
214  pointLayer_.setSize(mesh.nPoints());
215 
216  if (debug)
217  {
218  Info<< typeName << " : seeding "
219  << returnReduce(pp.nPoints(), sumOp<label>()) << " patch points"
220  << nl << endl;
221  }
222 
223  // Field on edges and points.
224  List<pointTopoDistanceData> edgeData(mesh.nEdges());
225  List<pointTopoDistanceData> pointData(mesh.nPoints());
226 
227  // Start of changes
228  labelList patchPoints(pp.nPoints());
229  List<pointTopoDistanceData> patchData(pp.nPoints());
230  forAll(pp.meshPoints(), patchPointi)
231  {
232  patchPoints[patchPointi] = pp.meshPoints()[patchPointi];
233  patchData[patchPointi] = pointTopoDistanceData(patchPointi, 0);
234  }
235 
236 
237  // Walk
238  PointEdgeWave<pointTopoDistanceData> distanceCalc
239  (
240  mesh,
241  patchPoints,
242  patchData,
243 
244  pointData,
245  edgeData,
246  mesh.globalData().nTotalPoints() // max iterations
247  );
248 
249  forAll(pointData, pointi)
250  {
251  pointToPatchPointAddressing_[pointi] = pointData[pointi].data();
252  pointLayer_[pointi] = pointData[pointi].distance();
253  }
254 
255 
256  // Derive from originating patch points what the patch edges were.
257  EdgeMap<label> pointsToEdge(pp.nEdges());
258  forAll(pp.edges(), edgeI)
259  {
260  pointsToEdge.insert(pp.edges()[edgeI], edgeI);
261  }
262 
263  // Look up on faces
264  forAll(faceToPatchEdgeAddressing_, facei)
265  {
266  if (faceToPatchEdgeAddressing_[facei] == labelMin)
267  {
268  // Face not yet done. Check if all points on same level
269  // or if not see what edge it originates from
270 
271  const face& f = mesh.faces()[facei];
272 
273  label levelI = pointLayer_[f[0]];
274  for (label fp = 1; fp < f.size(); fp++)
275  {
276  if (pointLayer_[f[fp]] != levelI)
277  {
278  levelI = -1;
279  break;
280  }
281  }
282 
283  if (levelI != -1)
284  {
285  // All same level
286  // Pout<< "Horizontal boundary face " << facei
287  // << " at:" << mesh.faceCentres()[facei]
288  // << " data:" << faceData[facei]
289  // << " pointDatas:"
290  // << UIndirectList<pointTopoDistanceData>(pointData, f)
291  // << endl;
292 
293  label patchFacei = faceData[facei].data();
294  label patchDist = faceData[facei].distance();
295 
296  faceToPatchEdgeAddressing_[facei] = -1;
297  faceToPatchFaceAddressing_[facei] = patchFacei+1;
298  faceLayer_[facei] = patchDist;
299  }
300  else
301  {
302  // Points of face on different levels
303 
304  // See if there is any edge
305  forAll(f, fp)
306  {
307  label pointi = f[fp];
308  label nextPointi = f.nextLabel(fp);
309 
310  EdgeMap<label>::const_iterator fnd = pointsToEdge.find
311  (
312  edge
313  (
314  pointData[pointi].data(),
315  pointData[nextPointi].data()
316  )
317  );
318  if (fnd != pointsToEdge.end())
319  {
320  faceToPatchEdgeAddressing_[facei] = fnd();
321  faceToPatchFaceAddressing_[facei] = 0;
322  label own = mesh.faceOwner()[facei];
323  faceLayer_[facei] = cellData[own].distance();
324 
325  // Note: could test whether the other edges on the
326  // face are consistent
327  break;
328  }
329  }
330  }
331  }
332  }
333  }
334 
335 
336 
337  // Use maps to find out mesh structure.
338  {
339  label nLayers = gMax(cellLayer_)+1;
340  labelListList layerToCells(invertOneToMany(nLayers, cellLayer_));
341 
342  structured_ = true;
343  forAll(layerToCells, layerI)
344  {
345  const labelList& lCells = layerToCells[layerI];
346 
347  forAll(lCells, lCelli)
348  {
349  label celli = lCells[lCelli];
350 
351  structured_ = isStructuredCell
352  (
353  mesh,
354  layerI,
355  celli
356  );
357 
358  if (!structured_)
359  {
360  break;
361  }
362  }
363 
364  if (!structured_)
365  {
366  break;
367  }
368  }
369 
370  reduce(structured_, andOp<bool>());
371  }
372 }
373 
374 
375 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
376 
378 (
379  const polyMesh& mesh,
380  const uindirectPrimitivePatch& pp
381 )
382 {
383  correct(mesh, pp);
384 }
385 
386 
387 // ************************************************************************* //
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
meshStructure(const polyMesh &mesh, const uindirectPrimitivePatch &)
Construct null.
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
A list of faces which address into the list of points.
Info<< "Predicted p max-min : "<< max(p).value()<< " "<< min(p).value()<< endl;rho==max(rho0+psi *p, rhoMin);# 1 "/home/ubuntu/OpenFOAM-9/applications/solvers/multiphase/cavitatingFoam/alphavPsi.H" 1{ alphav=max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));alphal=1.0 - alphav;Info<< "max-min alphav: "<< max(alphav).value()<< " "<< min(alphav).value()<< endl;psiModel-> correct()
Definition: pEqn.H:68
friend class const_iterator
Declare friendship with the const_iterator.
Definition: HashTable.H:197
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const char nl
Definition: Ostream.H:260
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
Foam::uindirectPrimitivePatch.
messageStream Info
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static const label labelMin
Definition: label.H:61
Namespace for OpenFOAM.