create3DCellShape.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-------------------------------------------------------------------------------
10License
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
26Description
27 Construct a cell shape from face information
28
29\*---------------------------------------------------------------------------*/
30
32#include "labelList.H"
33
34// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35
36namespace Foam
37{
38
39// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40
41cellShape create3DCellShape
42(
43 const label cellIndex,
44 const labelList& faceLabels,
45 const faceList& faces,
46 const labelList& owner,
47 const labelList& neighbour,
48 const label fluentCellModelID
49)
50{
51 // List of pointers to shape models for 3-D shape recognition
52 static List<const cellModel*> fluentCellModelLookup
53 (
54 7,
55 nullptr
56 );
57
58 fluentCellModelLookup[2] = cellModel::ptr(cellModel::TET);
59 fluentCellModelLookup[4] = cellModel::ptr(cellModel::HEX);
60 fluentCellModelLookup[5] = cellModel::ptr(cellModel::PYR);
61 fluentCellModelLookup[6] = cellModel::ptr(cellModel::PRISM);
62
63 static label faceMatchingOrder[7][6] =
64 {
65 {-1, -1, -1, -1, -1, -1},
66 {-1, -1, -1, -1, -1, -1},
67 { 0, 1, 2, 3, -1, -1}, // tet
68 {-1, -1, -1, -1, -1, -1},
69 { 0, 2, 4, 3, 5, 1}, // hex
70 { 0, 1, 2, 3, 4, -1}, // pyr
71 { 0, 2, 3, 4, 1, -1}, // prism
72 };
73
74 const cellModel& curModel = *fluentCellModelLookup[fluentCellModelID];
75
76 // Checking
77 if (faceLabels.size() != curModel.nFaces())
78 {
80 << "Number of face labels not equal to"
81 << "number of face in the model. "
82 << "Number of face labels: " << faceLabels.size()
83 << " number of faces in model: " << curModel.nFaces()
84 << abort(FatalError);
85 }
86
87 // make a list of outward-pointing faces
88 labelListList localFaces(faceLabels.size());
89
90 forAll(faceLabels, facei)
91 {
92 const label curFaceLabel = faceLabels[facei];
93
94 const labelList& curFace = faces[curFaceLabel];
95
96 if (owner[curFaceLabel] == cellIndex)
97 {
98 localFaces[facei] = curFace;
99 }
100 else if (neighbour[curFaceLabel] == cellIndex)
101 {
102 // Reverse the face
103 localFaces[facei].setSize(curFace.size());
104
105 forAllReverse(curFace, i)
106 {
107 localFaces[facei][curFace.size() - i - 1] =
108 curFace[i];
109 }
110 }
111 else
112 {
114 << "face " << curFaceLabel
115 << " does not belong to cell " << cellIndex
116 << ". Face owner: " << owner[curFaceLabel] << " neighbour: "
117 << neighbour[curFaceLabel]
118 << abort(FatalError);
119 }
120 }
121
122 // Algorithm:
123 // Make an empty list of pointLabels and initialise it with -1. Pick the
124 // first face from modelFaces and look through the faces to find one with
125 // the same number of labels. Insert face by copying its labels into
126 // pointLabels. Mark the face as used. Loop through all model faces.
127 // For each model face loop through faces. If the face is unused and the
128 // numbers of labels fit, try to match the face onto the point labels. If
129 // at least one edge is matched, insert the face into pointLabels. If at
130 // any stage the matching algorithm reaches the end of faces, the matching
131 // algorithm has failed. Once all the faces are matched, the list of
132 // pointLabels defines the model.
133
134 // Make a list of empty pointLabels
135 labelList pointLabels(curModel.nPoints(), -1);
136
137 // Follow the used mesh faces
138 List<bool> meshFaceUsed(localFaces.size(), false);
139
140 // Get the raw model faces
141 const faceList& modelFaces = curModel.modelFaces();
142
143 // Insert the first face into the list
144 const labelList& firstModelFace =
145 modelFaces[faceMatchingOrder[fluentCellModelID][0]];
146
147 bool found = false;
148
149 forAll(localFaces, meshFacei)
150 {
151 if (localFaces[meshFacei].size() == firstModelFace.size())
152 {
153 // Match. Insert points into the pointLabels
154 found = true;
155
156 const labelList& curMeshFace = localFaces[meshFacei];
157
158 meshFaceUsed[meshFacei] = true;
159
160 forAll(curMeshFace, pointi)
161 {
162 pointLabels[firstModelFace[pointi]] = curMeshFace[pointi];
163 }
164
165 break;
166 }
167 }
168
169 if (!found)
170 {
172 << "Cannot find match for first face. "
173 << "cell model: " << curModel.name() << " first model face: "
174 << firstModelFace << " Mesh faces: " << localFaces
175 << abort(FatalError);
176 }
177
178 for (label modelFacei = 1; modelFacei < modelFaces.size(); modelFacei++)
179 {
180 // get the next model face
181 const labelList& curModelFace =
182 modelFaces
183 [faceMatchingOrder[fluentCellModelID][modelFacei]];
184
185 found = false;
186
187 // Loop through mesh faces until a match is found
188 forAll(localFaces, meshFacei)
189 {
190 if
191 (
192 !meshFaceUsed[meshFacei]
193 && localFaces[meshFacei].size() == curModelFace.size()
194 )
195 {
196 // A possible match. A mesh face will be rotated, so make a copy
197 labelList meshFaceLabels = localFaces[meshFacei];
198
199 for
200 (
201 label rotation = 0;
202 rotation < meshFaceLabels.size();
203 rotation++
204 )
205 {
206 // try matching the face
207 label nMatchedLabels = 0;
208
209 forAll(meshFaceLabels, pointi)
210 {
211 if
212 (
213 pointLabels[curModelFace[pointi]]
214 == meshFaceLabels[pointi]
215 )
216 {
217 nMatchedLabels++;
218 }
219 }
220
221 if (nMatchedLabels >= 2)
222 {
223 // match!
224 found = true;
225 }
226
227 if (found)
228 {
229 // match found. Insert mesh face
230 forAll(meshFaceLabels, pointi)
231 {
232 pointLabels[curModelFace[pointi]] =
233 meshFaceLabels[pointi];
234 }
235
236 meshFaceUsed[meshFacei] = true;
237
238 break;
239 }
240 else
241 {
242 // No match found. Rotate face
243 label firstLabel = meshFaceLabels[0];
244
245 for (label i = 1; i < meshFaceLabels.size(); i++)
246 {
247 meshFaceLabels[i - 1] = meshFaceLabels[i];
248 }
249
250 meshFaceLabels.last() = firstLabel;
251 }
252 }
253
254 if (found) break;
255 }
256 }
257
258 if (!found)
259 {
260 // A model face is not matched. Shape detection failed
262 << "Cannot find match for face "
263 << modelFacei
264 << ".\nModel: " << curModel.name() << " model face: "
265 << curModelFace << " Mesh faces: " << localFaces
266 << "Matched points: " << pointLabels
267 << abort(FatalError);
268 }
269 }
270
271 return cellShape(curModel, pointLabels);
272}
273
274
275// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
276
277} // End namespace Foam
278
279// ************************************************************************* //
bool found
void setSize(const label n)
Alias for resize()
Definition: List.H:218
T & last()
Return the last element of the list.
Definition: UListI.H:216
@ PRISM
prism
Definition: cellModel.H:83
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
cellShape create3DCellShape(const label cellIndex, const labelList &faceLabels, const faceList &faces, const labelList &owner, const labelList &neighbour, const label fluentCellModelID)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:144
error FatalError
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
labelList pointLabels(nPoints, -1)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: stdFoam.H:346