Loading...
Searching...
No Matches
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-2023 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 "DynamicList.H"
35
36// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37
38template<class FaceList, class PointField>
39void
40Foam::PrimitivePatch<FaceList, PointField>::visitPointRegion
41(
42 const label pointi,
43 const labelUList& pFaces,
44 const label startFacei,
45 const label startEdgei,
46 UList<bool>& pFacesVisited
47) const
48{
49 const label index = pFaces.find(startFacei);
50
51 if (index >= 0 && !pFacesVisited[index])
52 {
53 // Mark face as visited
54 pFacesVisited[index] = true;
55
56 // Step to next edge on face which is still using pointi
57 label nextEdgei = -1;
58
59 for (const label faceEdgei : faceEdges()[startFacei])
60 {
61 const edge& e = edges()[faceEdgei];
62
63 if (faceEdgei != startEdgei && e.contains(pointi))
64 {
65 nextEdgei = faceEdgei;
66 break;
67 }
68 }
69
70 if (nextEdgei == -1)
71 {
73 << "Problem: cannot find edge out of "
74 << faceEdges()[startFacei]
75 << "on face " << startFacei << " that uses point " << pointi
76 << " and is not edge " << startEdgei << abort(FatalError);
77 }
78
79 // Walk to next face(s) across edge.
80 for (const label edgeFacei : edgeFaces()[nextEdgei])
81 {
82 if (edgeFacei != startFacei)
83 {
84 visitPointRegion
85 (
86 pointi,
87 pFaces,
88 edgeFacei,
89 nextEdgei,
90 pFacesVisited
91 );
92 }
93 }
94 }
95}
97
98// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
99
100template<class FaceList, class PointField>
103(
105) const
106{
107 if (badEdgesPtr)
108 {
109 badEdgesPtr->clear();
110 }
111
112 bool foundError = false;
113
115
116 const labelListList& edgeFcs = edgeFaces();
117
118 forAll(edgeFcs, edgei)
119 {
120 const label nNbrs = edgeFcs[edgei].size();
121
122 if (nNbrs < 1 || nNbrs > 2)
123 {
124 // Surface is illegal
125 foundError = true;
126
127 if (badEdgesPtr)
128 {
129 // Record and continue
130 if (nNbrs > 2)
131 {
132 badEdgesPtr->insert(edgei);
133 }
134 }
135 else
136 {
137 // Early exit when not recording bad edges
138 break;
139 }
140 }
141 else if (nNbrs == 1)
142 {
143 // Surface might be open or illegal so keep looping.
144 pType = surfaceTopo::OPEN;
145 }
146 }
147
148 if (foundError)
149 {
151 }
152
153 return pType;
154}
155
156
157template<class FaceList, class PointField>
158bool
160(
161 const bool report,
163) const
164{
165 // Check edgeFaces
166
167 bool foundError = false;
168
169 const labelListList& edgeFcs = edgeFaces();
170
171 forAll(edgeFcs, edgei)
172 {
173 const label nNbrs = edgeFcs[edgei].size();
174
175 if (nNbrs < 1 || nNbrs > 2)
176 {
177 foundError = true;
178
179 if (report)
180 {
181 Info<< "Edge " << edgei << " with vertices:" << edges()[edgei]
182 << " has " << nNbrs << " face neighbours" << endl;
183 }
184
185 if (pointSetPtr)
186 {
187 const edge& e = edges()[edgei];
188
189 pointSetPtr->insert(meshPoints()[e.first()]);
190 pointSetPtr->insert(meshPoints()[e.second()]);
191 }
192 }
193 }
194
196}
197
198
199template<class FaceList, class PointField>
200bool
202(
203 const bool report,
205) const
206{
207 const labelListList& pf = pointFaces();
208 const labelListList& pe = pointEdges();
209 const labelListList& ef = edgeFaces();
210 const labelList& mp = meshPoints();
211
212 bool foundError = false;
213
214 // Visited faces (as indices into pFaces)
215 DynamicList<bool> pFacesVisited;
216
217 forAll(pf, pointi)
218 {
219 const labelList& pFaces = pf[pointi];
220
221 pFacesVisited.resize_nocopy(pFaces.size());
222 pFacesVisited = false;
223
224 // Starting edge
225 const labelList& pEdges = pe[pointi];
226 const label startEdgei = pEdges[0];
227
228 const labelList& eFaces = ef[startEdgei];
229
230 for (const label edgeFacei : eFaces)
231 {
232 // Visit all faces using pointi, starting from eFaces[i] and
233 // startEdgei. Mark off all faces visited in pFacesVisited.
234 this->visitPointRegion
235 (
236 pointi,
237 pFaces,
238 edgeFacei, // starting face for walk
239 startEdgei, // starting edge for walk
240 pFacesVisited
241 );
242 }
243
244 // After this all faces using pointi should have been visited and
245 // marked off in pFacesVisited.
246
247 if (pFacesVisited.contains(false))
248 {
249 foundError = true;
250
251 const label meshPointi = mp[pointi];
252
253 if (pointSetPtr)
254 {
255 pointSetPtr->insert(meshPointi);
256 }
257
258 if (report)
259 {
260 Info<< "Point " << meshPointi
261 << " uses faces which are not connected through an edge"
262 << nl
263 << "This means that the surface formed by this patched"
264 << " is multiply connected at this point" << nl
265 << "Connected (patch) faces:" << nl;
266
267 forAll(pFacesVisited, i)
268 {
269 if (pFacesVisited[i])
270 {
271 Info<< " " << pFaces[i] << endl;
272 }
273 }
274
275 Info<< nl << "Unconnected (patch) faces:" << nl;
276 forAll(pFacesVisited, i)
277 {
278 if (!pFacesVisited[i])
279 {
280 Info<< " " << pFaces[i] << endl;
281 }
282 }
283 }
284 }
285 }
286
287 return foundError;
288}
289
290
291// ************************************************************************* //
labelHashSet * pointSetPtr
bool foundError
labelHashSet * badEdgesPtr
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition DynamicList.H:68
void resize_nocopy(const label len)
Alter addressable list size, allocating new space if required without necessarily recovering old cont...
const labelListList & pointEdges() const
Return point-edge addressing.
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
bool checkTopology(const bool report=false, labelHashSet *pointSetPtr=nullptr) const
Check surface formed by patch for manifoldness (see above).
const labelList & meshPoints() const
Return labelList of mesh points in patch.
bool checkPointManifold(const bool report=false, labelHashSet *pointSetPtr=nullptr) const
Checks primitivePatch for faces sharing point but not edge.
const labelListList & pointFaces() const
Return point-face addressing.
surfaceTopo surfaceType(labelHashSet *badEdgesPtr=nullptr) const
Calculate surface type formed by patch, optionally recording the indices of illegal edges.
const labelListList & edgeFaces() const
Return edge-face addressing.
surfaceTopo
Enumeration defining the surface type. Used in check routines.
bool contains(const T &val) const
True if the value is contained in the list.
Definition UListI.H:302
void size(const label n)
Older name for setAddressableSize.
Definition UList.H:118
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
Definition edge.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition error.H:600
List< labelList > labelListList
List of labelList.
Definition labelList.H:38
List< label > labelList
A List of labels.
Definition List.H:62
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition HashSet.H:85
messageStream Info
Information stream (stdout output on master, null elsewhere).
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition Ostream.H:519
errorManip< error > abort(error &err)
Definition errorManip.H:139
constexpr char nl
The newline '\n' character (0x0a).
Definition Ostream.H:50
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];} pointMap[start]=pointMap[end]=Foam::min(start, end);} } 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};constexpr 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< DynamicList< face > > pFaces[nBCs]
volScalarField & e
#define forAll(list, i)
Loop across all elements in list.
Definition stdFoam.H:299