PrimitivePatchCheck.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) 2011-2016 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 Description
27  Checks topology of the patch.
28 
29 \*---------------------------------------------------------------------------*/
30 
31 #include "PrimitivePatch.H"
32 #include "Map.H"
33 #include "ListOps.H"
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
37 template
38 <
39  class Face,
40  template<class> class FaceList,
41  class PointField,
42  class PointType
43 >
44 void
47 (
48  const label pointi,
49  const labelList& pFaces,
50  const label startFacei,
51  const label startEdgeI,
52  boolList& pFacesHad
53 ) const
54 {
55  label index = pFaces.find(startFacei);
56 
57  if (!pFacesHad[index])
58  {
59  // Mark face as been visited.
60  pFacesHad[index] = true;
61 
62  // Step to next edge on face which is still using pointi
63  const labelList& fEdges = faceEdges()[startFacei];
64 
65  label nextEdgeI = -1;
66 
67  forAll(fEdges, i)
68  {
69  label edgeI = fEdges[i];
70 
71  const edge& e = edges()[edgeI];
72 
73  if (edgeI != startEdgeI && (e[0] == pointi || e[1] == pointi))
74  {
75  nextEdgeI = edgeI;
76 
77  break;
78  }
79  }
80 
81  if (nextEdgeI == -1)
82  {
84  << "Problem: cannot find edge out of " << fEdges
85  << "on face " << startFacei << " that uses point " << pointi
86  << " and is not edge " << startEdgeI << abort(FatalError);
87  }
88 
89  // Walk to next face(s) across edge.
90  const labelList& eFaces = edgeFaces()[nextEdgeI];
91 
92  forAll(eFaces, i)
93  {
94  if (eFaces[i] != startFacei)
95  {
96  visitPointRegion
97  (
98  pointi,
99  pFaces,
100  eFaces[i],
101  nextEdgeI,
102  pFacesHad
103  );
104  }
105  }
106  }
107 }
108 
109 
110 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
111 
112 template
113 <
114  class Face,
115  template<class> class FaceList,
116  class PointField,
117  class PointType
118 >
119 typename
122 surfaceType() const
123 {
124  if (debug)
125  {
126  InfoInFunction << "Calculating patch topology" << endl;
127  }
128 
129  const labelListList& edgeFcs = edgeFaces();
130 
131  surfaceTopo pType = MANIFOLD;
132 
133  forAll(edgeFcs, edgeI)
134  {
135  label nNbrs = edgeFcs[edgeI].size();
136 
137  if (nNbrs < 1 || nNbrs > 2)
138  {
139  pType = ILLEGAL;
140 
141  // Can exit now. Surface is illegal.
142  return pType;
143  }
144  else if (nNbrs == 1)
145  {
146  // Surface might be open or illegal so keep looping.
147  pType = OPEN;
148  }
149  }
150 
151  if (debug)
152  {
153  Info<< " Finished." << endl;
154  }
155 
156  return pType;
157 }
158 
159 
160 template
161 <
162  class Face,
163  template<class> class FaceList,
164  class PointField,
165  class PointType
166 >
167 bool
170 (
171  const bool report,
172  labelHashSet* setPtr
173 ) const
174 {
175  if (debug)
176  {
177  InfoInFunction << "Checking patch topology" << endl;
178  }
179 
180  // Check edgeFaces
181 
182  const labelListList& edgeFcs = edgeFaces();
183 
184  bool illegalTopo = false;
185 
186  forAll(edgeFcs, edgeI)
187  {
188  label nNbrs = edgeFcs[edgeI].size();
189 
190  if (nNbrs < 1 || nNbrs > 2)
191  {
192  illegalTopo = true;
193 
194  if (report)
195  {
196  Info<< "Edge " << edgeI << " with vertices:" << edges()[edgeI]
197  << " has " << nNbrs << " face neighbours"
198  << endl;
199  }
200 
201  if (setPtr)
202  {
203  const edge& e = edges()[edgeI];
204 
205  setPtr->insert(meshPoints()[e.start()]);
206  setPtr->insert(meshPoints()[e.end()]);
207  }
208  }
209  }
210 
211  if (debug)
212  {
213  Info<< " Finished." << endl;
214  }
215 
216  return illegalTopo;
217 }
218 
219 
220 template
221 <
222  class Face,
223  template<class> class FaceList,
224  class PointField,
225  class PointType
226 >
227 bool
230 (
231  const bool report,
232  labelHashSet* setPtr
233 ) const
234 {
235  const labelListList& pf = pointFaces();
236  const labelListList& pe = pointEdges();
237  const labelListList& ef = edgeFaces();
238  const labelList& mp = meshPoints();
239 
240  bool foundError = false;
241 
242  forAll(pf, pointi)
243  {
244  const labelList& pFaces = pf[pointi];
245 
246  // Visited faces (as indices into pFaces)
247  boolList pFacesHad(pFaces.size(), false);
248 
249  // Starting edge
250  const labelList& pEdges = pe[pointi];
251  label startEdgeI = pEdges[0];
252 
253  const labelList& eFaces = ef[startEdgeI];
254 
255  forAll(eFaces, i)
256  {
257  // Visit all faces using pointi, starting from eFaces[i] and
258  // startEdgeI. Mark off all faces visited in pFacesHad.
259  this->visitPointRegion
260  (
261  pointi,
262  pFaces,
263  eFaces[i], // starting face for walk
264  startEdgeI, // starting edge for walk
265  pFacesHad
266  );
267  }
268 
269  // After this all faces using pointi should have been visited and
270  // marked off in pFacesHad.
271 
272  label unset = pFacesHad.find(false);
273 
274  if (unset != -1)
275  {
276  foundError = true;
277 
278  label meshPointi = mp[pointi];
279 
280  if (setPtr)
281  {
282  setPtr->insert(meshPointi);
283  }
284 
285  if (report)
286  {
287  Info<< "Point " << meshPointi
288  << " uses faces which are not connected through an edge"
289  << nl
290  << "This means that the surface formed by this patched"
291  << " is multiply connected at this point" << nl
292  << "Connected (patch) faces:" << nl;
293 
294  forAll(pFacesHad, i)
295  {
296  if (pFacesHad[i])
297  {
298  Info<< " " << pFaces[i] << endl;
299  }
300  }
301 
302  Info<< nl << "Unconnected (patch) faces:" << nl;
303  forAll(pFacesHad, i)
304  {
305  if (!pFacesHad[i])
306  {
307  Info<< " " << pFaces[i] << endl;
308  }
309  }
310  }
311  }
312  }
313 
314  return foundError;
315 }
316 
317 
318 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:74
Foam::constant::atomic::mp
const dimensionedScalar mp
Proton mass.
InfoInFunction
#define InfoInFunction
Report an information message using Foam::Info.
Definition: messageStream.H:316
Foam::edge
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:63
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:72
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:337
Foam::HashSet< label, Hash< label > >
Foam::PrimitivePatch::surfaceType
surfaceTopo surfaceType() const
Calculate surface type formed by patch.
Definition: PrimitivePatchCheck.C:122
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:290
Map.H
Foam::label
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:62
pFaces
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
PrimitivePatch.H
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
Foam::FatalError
error FatalError
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:137
Foam::PrimitivePatch::surfaceTopo
surfaceTopo
Enumeration defining the surface type. Used in check routines.
Definition: PrimitivePatch.H:108
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:355
Foam::PrimitivePatch::checkPointManifold
bool checkPointManifold(const bool report=false, labelHashSet *setPtr=nullptr) const
Checks primitivePatch for faces sharing point but not edge.
Definition: PrimitivePatchCheck.C:230
Foam::nl
constexpr char nl
Definition: Ostream.H:372
Foam::List< labelList >
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::HashSet::insert
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:182
ListOps.H
Various functions to operate on Lists.
Foam::PrimitivePatch::checkTopology
bool checkTopology(const bool report=false, labelHashSet *setPtr=nullptr) const
Check surface formed by patch for manifoldness (see above).
Definition: PrimitivePatchCheck.C:170
Foam::PrimitivePatch
A list of faces which address into the list of points.
Definition: PrimitivePatch.H:90