meshStructure.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2013-2016 OpenFOAM Foundation
9  Copyright (C) 2019-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "meshStructure.H"
30 #include "FaceCellWave.H"
31 #include "topoDistanceData.H"
32 #include "pointTopoDistanceData.H"
33 #include "PointEdgeWave.H"
34 #include "globalIndex.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40  defineTypeNameAndDebug(meshStructure, 0);
41 }
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 bool Foam::meshStructure::isStructuredCell
47 (
48  const polyMesh& mesh,
49  const label layerI,
50  const label celli
51 ) const
52 {
53  const cell& cFaces = mesh.cells()[celli];
54 
55  // Count number of side faces
56  label nSide = 0;
57  forAll(cFaces, i)
58  {
59  if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
60  {
61  nSide++;
62  }
63  }
64 
65  if (nSide != cFaces.size()-2)
66  {
67  return false;
68  }
69 
70  // Check that side faces have correct point layers
71  forAll(cFaces, i)
72  {
73  if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
74  {
75  const face& f = mesh.faces()[cFaces[i]];
76 
77  label nLayer = 0;
78  label nLayerPlus1 = 0;
79  forAll(f, fp)
80  {
81  label pointi = f[fp];
82  if (pointLayer_[pointi] == layerI)
83  {
84  nLayer++;
85  }
86  else if (pointLayer_[pointi] == layerI+1)
87  {
88  nLayerPlus1++;
89  }
90  }
91 
92  if (f.size() != 4 || (nLayer+nLayerPlus1 != 4))
93  {
94  return false;
95  }
96  }
97  }
98 
99  return true;
100 }
101 
102 
103 void Foam::meshStructure::correct
104 (
105  const polyMesh& mesh,
106  const uindirectPrimitivePatch& pp,
107  const globalIndex& globalFaces,
108  const globalIndex& globalEdges,
109  const globalIndex& globalPoints
110 )
111 {
112  // Field on cells and faces.
113  List<topoDistanceData<label>> cellData(mesh.nCells());
114  List<topoDistanceData<label>> faceData(mesh.nFaces());
115 
116  {
117  if (debug)
118  {
119  Info<< typeName << " : seeding "
120  << returnReduce(pp.size(), sumOp<label>()) << " patch faces"
121  << nl << endl;
122  }
123 
124 
125  // Start of changes
126  labelList patchFaces(pp.size());
127  List<topoDistanceData<label>> patchData(pp.size());
128  forAll(pp, patchFacei)
129  {
130  patchFaces[patchFacei] = pp.addressing()[patchFacei];
131  patchData[patchFacei] = topoDistanceData<label>
132  (
133  0, // distance
134  globalFaces.toGlobal(patchFacei) // passive data
135  );
136  }
137 
138 
139  // Propagate information inwards
140  FaceCellWave<topoDistanceData<label>> distanceCalc
141  (
142  mesh,
143  patchFaces,
144  patchData,
145  faceData,
146  cellData,
148  );
149 
150 
151  // Determine cells from face-cell-walk
152  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
153 
154  cellToPatchFaceAddressing_.setSize(mesh.nCells());
155  cellLayer_.setSize(mesh.nCells());
156  forAll(cellToPatchFaceAddressing_, celli)
157  {
158  cellToPatchFaceAddressing_[celli] = cellData[celli].data();
159  cellLayer_[celli] = cellData[celli].distance();
160  }
161 
162 
163 
164  // Determine faces from face-cell-walk
165  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
166 
167  faceToPatchFaceAddressing_.setSize(mesh.nFaces());
168  faceToPatchEdgeAddressing_.setSize(mesh.nFaces());
169  faceToPatchEdgeAddressing_ = labelMin;
170  faceLayer_.setSize(mesh.nFaces());
171 
172  forAll(faceToPatchFaceAddressing_, facei)
173  {
174  label own = mesh.faceOwner()[facei];
175  label patchFacei = faceData[facei].data();
176  label patchDist = faceData[facei].distance();
177 
178  if (mesh.isInternalFace(facei))
179  {
180  label nei = mesh.faceNeighbour()[facei];
181 
182  if (cellData[own].distance() == cellData[nei].distance())
183  {
184  // side face
185  faceToPatchFaceAddressing_[facei] = 0;
186  faceLayer_[facei] = cellData[own].distance();
187  }
188  else if (cellData[own].distance() < cellData[nei].distance())
189  {
190  // unturned face
191  faceToPatchFaceAddressing_[facei] = patchFacei+1;
192  faceToPatchEdgeAddressing_[facei] = -1;
193  faceLayer_[facei] = patchDist;
194  }
195  else
196  {
197  // turned face
198  faceToPatchFaceAddressing_[facei] = -(patchFacei+1);
199  faceToPatchEdgeAddressing_[facei] = -1;
200  faceLayer_[facei] = patchDist;
201  }
202  }
203  else if (patchDist == cellData[own].distance())
204  {
205  // starting face
206  faceToPatchFaceAddressing_[facei] = -(patchFacei+1);
207  faceToPatchEdgeAddressing_[facei] = -1;
208  faceLayer_[facei] = patchDist;
209  }
210  else
211  {
212  // unturned face or side face. Cannot be determined until
213  // we determine the point layers. Problem is that both are
214  // the same number of steps away from the initial seed face.
215  }
216  }
217  }
218 
219 
220  // Determine points from separate walk on point-edge
221  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
222 
223  {
224  pointToPatchPointAddressing_.setSize(mesh.nPoints());
225  pointLayer_.setSize(mesh.nPoints());
226 
227  if (debug)
228  {
229  Info<< typeName << " : seeding "
230  << returnReduce(pp.nPoints(), sumOp<label>()) << " patch points"
231  << nl << endl;
232  }
233 
234  // Field on edges and points.
235  List<pointTopoDistanceData<label>> edgeData(mesh.nEdges());
236  List<pointTopoDistanceData<label>> pointData(mesh.nPoints());
237 
238  // Start of changes
239  labelList patchPoints(pp.nPoints());
240  List<pointTopoDistanceData<label>> patchData(pp.nPoints());
241  forAll(pp.meshPoints(), patchPointi)
242  {
243  patchPoints[patchPointi] = pp.meshPoints()[patchPointi];
244  patchData[patchPointi] = pointTopoDistanceData<label>
245  (
246  0, // distance
247  globalPoints.toGlobal(patchPointi) // passive data
248  );
249  }
250 
251 
252  // Walk
253  PointEdgeWave<pointTopoDistanceData<label>> distanceCalc
254  (
255  mesh,
256  patchPoints,
257  patchData,
258 
259  pointData,
260  edgeData,
261  mesh.globalData().nTotalPoints() // max iterations
262  );
263 
264  forAll(pointData, pointi)
265  {
266  pointToPatchPointAddressing_[pointi] = pointData[pointi].data();
267  pointLayer_[pointi] = pointData[pointi].distance();
268  }
269 
270 
271  // Derive from originating patch points what the patch edges were.
272  EdgeMap<label> pointsToEdge(pp.nEdges());
273  forAll(pp.edges(), edgeI)
274  {
275  const edge& e = pp.edges()[edgeI];
276  edge globalEdge
277  (
278  globalPoints.toGlobal(e[0]),
279  globalPoints.toGlobal(e[1])
280  );
281  pointsToEdge.insert(globalEdge, globalEdges.toGlobal(edgeI));
282  }
283 
284  // Look up on faces
285  forAll(faceToPatchEdgeAddressing_, facei)
286  {
287  if (faceToPatchEdgeAddressing_[facei] == labelMin)
288  {
289  // Face not yet done. Check if all points on same level
290  // or if not see what edge it originates from
291 
292  const face& f = mesh.faces()[facei];
293 
294  label levelI = pointLayer_[f[0]];
295  for (label fp = 1; fp < f.size(); fp++)
296  {
297  if (pointLayer_[f[fp]] != levelI)
298  {
299  levelI = -1;
300  break;
301  }
302  }
303 
304  if (levelI != -1)
305  {
306  // All same level
307  //Pout<< "Horizontal boundary face " << facei
308  // << " at:" << mesh.faceCentres()[facei]
309  // << " data:" << faceData[facei]
310  // << " pointDatas:"
311  // << UIndirectList<pointTopoDistanceData<label>>
312  // (pointData, f)
313  // << endl;
314 
315  label patchFacei = faceData[facei].data();
316  label patchDist = faceData[facei].distance();
317 
318  faceToPatchEdgeAddressing_[facei] = -1;
319  faceToPatchFaceAddressing_[facei] = patchFacei+1;
320  faceLayer_[facei] = patchDist;
321  }
322  else
323  {
324  // Points of face on different levels
325 
326  // See if there is any edge
327  forAll(f, fp)
328  {
329  label pointi = f[fp];
330  label nextPointi = f.nextLabel(fp);
331 
332  const auto fnd = pointsToEdge.cfind
333  (
334  edge
335  (
336  pointData[pointi].data(),
337  pointData[nextPointi].data()
338  )
339  );
340 
341  if (fnd.found())
342  {
343  faceToPatchEdgeAddressing_[facei] = fnd.val();
344  faceToPatchFaceAddressing_[facei] = 0;
345  label own = mesh.faceOwner()[facei];
346  faceLayer_[facei] = cellData[own].distance();
347 
348  // Note: could test whether the other edges on the
349  // face are consistent
350  break;
351  }
352  }
353  }
354  }
355  }
356  }
357 
358 
359 
360  // Use maps to find out mesh structure.
361  {
362  label nLayers = gMax(cellLayer_)+1;
363  labelListList layerToCells(invertOneToMany(nLayers, cellLayer_));
364 
365  structured_ = true;
366  forAll(layerToCells, layerI)
367  {
368  const labelList& lCells = layerToCells[layerI];
369 
370  forAll(lCells, lCelli)
371  {
372  label celli = lCells[lCelli];
373 
374  structured_ = isStructuredCell
375  (
376  mesh,
377  layerI,
378  celli
379  );
380 
381  if (!structured_)
382  {
383  break;
384  }
385  }
386 
387  if (!structured_)
388  {
389  break;
390  }
391  }
392 
393  reduce(structured_, andOp<bool>());
394  }
395 }
396 
397 
398 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
399 
401 (
402  const polyMesh& mesh,
403  const uindirectPrimitivePatch& pp
404 )
405 {
406  correct
407  (
408  mesh,
409  pp,
410  globalIndex(pp.size()),
411  globalIndex(pp.nEdges()),
412  globalIndex(pp.nPoints())
413  );
414 }
415 
416 
418 (
419  const polyMesh& mesh,
420  const uindirectPrimitivePatch& pp,
421  const globalIndex& globalFaces,
422  const globalIndex& globalEdges,
424 )
425 {
426  correct
427  (
428  mesh,
429  pp,
430  globalFaces,
431  globalEdges,
433  );
434 }
435 
436 
437 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:71
topoDistanceData.H
Foam::globalPoints
Calculates points shared by more than two processor patches or cyclic patches.
Definition: globalPoints.H:102
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
Foam::primitiveMesh::nFaces
label nFaces() const
Number of mesh faces.
Definition: primitiveMeshI.H:90
globalIndex.H
Foam::PrimitivePatch::nEdges
label nEdges() const
Return number of edges in patch.
Definition: PrimitivePatch.H:322
Foam::globalMeshData::nTotalPoints
label nTotalPoints() const
Return total number of points in decomposed mesh. Not.
Definition: globalMeshData.H:381
Foam::meshStructure::meshStructure
meshStructure(const polyMesh &mesh, const uindirectPrimitivePatch &)
Construct from mesh and faces in mesh. Any addressing to.
Definition: meshStructure.C:401
Foam::globalMeshData::nTotalCells
label nTotalCells() const
Return total number of cells in decomposed mesh.
Definition: globalMeshData.H:394
Foam::primitiveMesh::nEdges
label nEdges() const
Number of mesh edges.
Definition: primitiveMeshI.H:67
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
meshStructure.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
Foam::primitiveMesh::nCells
label nCells() const
Number of mesh cells.
Definition: primitiveMeshI.H:96
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
PointEdgeWave.H
correct
fvOptions correct(rho)
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1076
Foam::PrimitivePatch::nPoints
label nPoints() const
Return number of points supporting patch faces.
Definition: PrimitivePatch.H:316
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::globalIndex
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
Foam::distance
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
Foam::polyMesh::faces
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1063
pointTopoDistanceData.H
Foam::nl
constexpr char nl
Definition: Ostream.H:385
f
labelList f(nPoints)
Foam::primitiveMesh::isInternalFace
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
Definition: primitiveMeshI.H:102
FaceCellWave.H
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::primitiveMesh::nPoints
label nPoints() const
Number of mesh points.
Definition: primitiveMeshI.H:37
Foam::labelMin
constexpr label labelMin
Definition: label.H:60
Foam::uindirectPrimitivePatch
PrimitivePatch< UIndirectList< face >, const pointField & > uindirectPrimitivePatch
A PrimitivePatch with UIndirectList for the faces, const reference for the point field.
Definition: uindirectPrimitivePatch.H:49
Foam::polyMesh::globalData
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1234
Foam::List::setSize
void setSize(const label newSize)
Alias for resize(const label)
Definition: ListI.H:146
Foam::invertOneToMany
labelListList invertOneToMany(const label len, const labelUList &map)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:114
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1082
Foam::PrimitivePatch
A list of faces which address into the list of points.
Definition: PrimitivePatch.H:85