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 Copyright (C) 2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
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
27Description
28 Checks topology of the patch.
29
30\*---------------------------------------------------------------------------*/
31
32#include "PrimitivePatch.H"
33#include "Map.H"
34#include "ListOps.H"
35
36// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37
38template<class FaceList, class PointField>
39void
41(
42 const label pointi,
43 const labelList& pFaces,
44 const label startFacei,
45 const label startEdgeI,
46 boolList& pFacesHad
47) const
48{
49 label index = pFaces.find(startFacei);
50
51 if (!pFacesHad[index])
52 {
53 // Mark face as been visited.
54 pFacesHad[index] = true;
55
56 // Step to next edge on face which is still using pointi
57 const labelList& fEdges = faceEdges()[startFacei];
58
59 label nextEdgeI = -1;
60
61 forAll(fEdges, i)
62 {
63 label edgeI = fEdges[i];
64
65 const edge& e = edges()[edgeI];
66
67 if (edgeI != startEdgeI && (e[0] == pointi || e[1] == pointi))
68 {
69 nextEdgeI = edgeI;
70
71 break;
72 }
73 }
74
75 if (nextEdgeI == -1)
76 {
78 << "Problem: cannot find edge out of " << fEdges
79 << "on face " << startFacei << " that uses point " << pointi
80 << " and is not edge " << startEdgeI << abort(FatalError);
81 }
82
83 // Walk to next face(s) across edge.
84 const labelList& eFaces = edgeFaces()[nextEdgeI];
85
86 forAll(eFaces, i)
87 {
88 if (eFaces[i] != startFacei)
89 {
90 visitPointRegion
91 (
92 pointi,
93 pFaces,
94 eFaces[i],
95 nextEdgeI,
96 pFacesHad
97 );
98 }
99 }
100 }
101}
102
103
104// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
105
106template<class FaceList, class PointField>
109{
110 DebugInFunction << "Calculating patch topology" << nl;
111
112 const labelListList& edgeFcs = edgeFaces();
113
114 surfaceTopo pType = MANIFOLD;
115
116 forAll(edgeFcs, edgeI)
117 {
118 label nNbrs = edgeFcs[edgeI].size();
119
120 if (nNbrs < 1 || nNbrs > 2)
121 {
122 pType = ILLEGAL;
123
124 // Can exit now. Surface is illegal.
125 return pType;
126 }
127 else if (nNbrs == 1)
128 {
129 // Surface might be open or illegal so keep looping.
130 pType = OPEN;
131 }
132 }
133
134 DebugInFunction << "Calculated patch topology" << nl;
135
136 return pType;
137}
138
139
140template<class FaceList, class PointField>
141bool
143(
144 const bool report,
145 labelHashSet* setPtr
146) const
147{
148 DebugInFunction << "Checking patch topology" << nl;
149
150 // Check edgeFaces
151
152 const labelListList& edgeFcs = edgeFaces();
153
154 bool illegalTopo = false;
155
156 forAll(edgeFcs, edgeI)
157 {
158 label nNbrs = edgeFcs[edgeI].size();
159
160 if (nNbrs < 1 || nNbrs > 2)
161 {
162 illegalTopo = true;
163
164 if (report)
165 {
166 Info<< "Edge " << edgeI << " with vertices:" << edges()[edgeI]
167 << " has " << nNbrs << " face neighbours"
168 << endl;
169 }
170
171 if (setPtr)
172 {
173 const edge& e = edges()[edgeI];
174
175 setPtr->insert(meshPoints()[e.start()]);
176 setPtr->insert(meshPoints()[e.end()]);
177 }
178 }
179 }
180
181 DebugInFunction << "Checked patch topology" << nl;
182
183 return illegalTopo;
184}
185
186
187template<class FaceList, class PointField>
188bool
190(
191 const bool report,
192 labelHashSet* setPtr
193) const
194{
195 const labelListList& pf = pointFaces();
196 const labelListList& pe = pointEdges();
197 const labelListList& ef = edgeFaces();
198 const labelList& mp = meshPoints();
199
200 bool foundError = false;
201
202 forAll(pf, pointi)
203 {
204 const labelList& pFaces = pf[pointi];
205
206 // Visited faces (as indices into pFaces)
207 boolList pFacesHad(pFaces.size(), false);
208
209 // Starting edge
210 const labelList& pEdges = pe[pointi];
211 label startEdgeI = pEdges[0];
212
213 const labelList& eFaces = ef[startEdgeI];
214
215 forAll(eFaces, i)
216 {
217 // Visit all faces using pointi, starting from eFaces[i] and
218 // startEdgeI. Mark off all faces visited in pFacesHad.
219 this->visitPointRegion
220 (
221 pointi,
222 pFaces,
223 eFaces[i], // starting face for walk
224 startEdgeI, // starting edge for walk
225 pFacesHad
226 );
227 }
228
229 // After this all faces using pointi should have been visited and
230 // marked off in pFacesHad.
231
232 label unset = pFacesHad.find(false);
233
234 if (unset != -1)
235 {
236 foundError = true;
237
238 label meshPointi = mp[pointi];
239
240 if (setPtr)
241 {
242 setPtr->insert(meshPointi);
243 }
244
245 if (report)
246 {
247 Info<< "Point " << meshPointi
248 << " uses faces which are not connected through an edge"
249 << nl
250 << "This means that the surface formed by this patched"
251 << " is multiply connected at this point" << nl
252 << "Connected (patch) faces:" << nl;
253
254 forAll(pFacesHad, i)
255 {
256 if (pFacesHad[i])
257 {
258 Info<< " " << pFaces[i] << endl;
259 }
260 }
261
262 Info<< nl << "Unconnected (patch) faces:" << nl;
263 forAll(pFacesHad, i)
264 {
265 if (!pFacesHad[i])
266 {
267 Info<< " " << pFaces[i] << endl;
268 }
269 }
270 }
271 }
272 }
273
274 return foundError;
275}
276
277
278// ************************************************************************* //
Various functions to operate on Lists.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
A list of faces which address into the list of points.
bool checkPointManifold(const bool report=false, labelHashSet *setPtr=nullptr) const
Checks primitivePatch for faces sharing point but not edge.
surfaceTopo surfaceType() const
Calculate surface type formed by patch.
bool checkTopology(const bool report=false, labelHashSet *setPtr=nullptr) const
Check surface formed by patch for manifoldness (see above).
surfaceTopo
Enumeration defining the surface type. Used in check routines.
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
#define DebugInFunction
Report an information message using Foam::Info.
List< label > labelList
A List of labels.
Definition: List.H:66
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
errorManip< error > abort(error &err)
Definition: errorManip.H:144
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
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].reset(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
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333