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 | Copyright (C) 2011 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Description
25  Construct a cell shape from face information
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "cellShapeRecognition.H"
30 #include "labelList.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 cellShape create3DCellShape
40 (
41  const label cellIndex,
42  const labelList& faceLabels,
43  const faceList& faces,
44  const labelList& owner,
45  const labelList& neighbour,
46  const label fluentCellModelID
47 )
48 {
49  // List of pointers to shape models for 3-D shape recognition
50  static List<const cellModel*> fluentCellModelLookup
51  (
52  7,
53  reinterpret_cast<const cellModel*>(0)
54  );
55 
56  fluentCellModelLookup[2] = cellModeller::lookup("tet");
57  fluentCellModelLookup[4] = cellModeller::lookup("hex");
58  fluentCellModelLookup[5] = cellModeller::lookup("pyr");
59  fluentCellModelLookup[6] = cellModeller::lookup("prism");
60 
61  static label faceMatchingOrder[7][6] =
62  {
63  {-1, -1, -1, -1, -1, -1},
64  {-1, -1, -1, -1, -1, -1},
65  { 0, 1, 2, 3, -1, -1}, // tet
66  {-1, -1, -1, -1, -1, -1},
67  { 0, 2, 4, 3, 5, 1}, // hex
68  { 0, 1, 2, 3, 4, -1}, // pyr
69  { 0, 2, 3, 4, 1, -1}, // prism
70  };
71 
72  const cellModel& curModel = *fluentCellModelLookup[fluentCellModelID];
73 
74  // Checking
75  if (faceLabels.size() != curModel.nFaces())
76  {
78  (
79  "create3DCellShape(const label cellIndex, "
80  "const labelList& faceLabels, const labelListList& faces, "
81  "const labelList& owner, const labelList& neighbour, "
82  "const label fluentCellModelID)"
83  ) << "Number of face labels not equal to"
84  << "number of face in the model. "
85  << "Number of face labels: " << faceLabels.size()
86  << " number of faces in model: " << curModel.nFaces()
87  << abort(FatalError);
88  }
89 
90  // make a list of outward-pointing faces
91  labelListList localFaces(faceLabels.size());
92 
93  forAll(faceLabels, faceI)
94  {
95  const label curFaceLabel = faceLabels[faceI];
96 
97  const labelList& curFace = faces[curFaceLabel];
98 
99  if (owner[curFaceLabel] == cellIndex)
100  {
101  localFaces[faceI] = curFace;
102  }
103  else if (neighbour[curFaceLabel] == cellIndex)
104  {
105  // Reverse the face
106  localFaces[faceI].setSize(curFace.size());
107 
108  forAllReverse(curFace, i)
109  {
110  localFaces[faceI][curFace.size() - i - 1] =
111  curFace[i];
112  }
113  }
114  else
115  {
117  (
118  "create3DCellShape(const label cellIndex, "
119  "const labelList& faceLabels, const labelListList& faces, "
120  "const labelList& owner, const labelList& neighbour, "
121  "const label fluentCellModelID)"
122  ) << "face " << curFaceLabel
123  << " does not belong to cell " << cellIndex
124  << ". Face owner: " << owner[curFaceLabel] << " neighbour: "
125  << neighbour[curFaceLabel]
126  << abort(FatalError);
127  }
128  }
129 
130  // Algorithm:
131  // Make an empty list of pointLabels and initialise it with -1. Pick the
132  // first face from modelFaces and look through the faces to find one with
133  // the same number of labels. Insert face by copying its labels into
134  // pointLabels. Mark the face as used. Loop through all model faces.
135  // For each model face loop through faces. If the face is unused and the
136  // numbers of labels fit, try to match the face onto the point labels. If
137  // at least one edge is matched, insert the face into pointLabels. If at
138  // any stage the matching algorithm reaches the end of faces, the matching
139  // algorithm has failed. Once all the faces are matched, the list of
140  // pointLabels defines the model.
141 
142  // Make a list of empty pointLabels
143  labelList pointLabels(curModel.nPoints(), -1);
144 
145  // Follow the used mesh faces
146  List<bool> meshFaceUsed(localFaces.size(), false);
147 
148  // Get the raw model faces
149  const faceList& modelFaces = curModel.modelFaces();
150 
151  // Insert the first face into the list
152  const labelList& firstModelFace =
153  modelFaces[faceMatchingOrder[fluentCellModelID][0]];
154 
155  bool found = false;
156 
157  forAll(localFaces, meshFaceI)
158  {
159  if (localFaces[meshFaceI].size() == firstModelFace.size())
160  {
161  // Match. Insert points into the pointLabels
162  found = true;
163 
164  const labelList& curMeshFace = localFaces[meshFaceI];
165 
166  meshFaceUsed[meshFaceI] = true;
167 
168  forAll(curMeshFace, pointI)
169  {
170  pointLabels[firstModelFace[pointI]] = curMeshFace[pointI];
171  }
172 
173  break;
174  }
175  }
176 
177  if (!found)
178  {
180  (
181  "create3DCellShape(const label cellIndex, "
182  "const labelList& faceLabels, const labelListList& faces, "
183  "const labelList& owner, const labelList& neighbour, "
184  "const label fluentCellModelID)"
185  ) << "Cannot find match for first face. "
186  << "cell model: " << curModel.name() << " first model face: "
187  << firstModelFace << " Mesh faces: " << localFaces
188  << abort(FatalError);
189  }
190 
191  for (label modelFaceI = 1; modelFaceI < modelFaces.size(); modelFaceI++)
192  {
193  // get the next model face
194  const labelList& curModelFace =
195  modelFaces
196  [faceMatchingOrder[fluentCellModelID][modelFaceI]];
197 
198  found = false;
199 
200  // Loop through mesh faces until a match is found
201  forAll(localFaces, meshFaceI)
202  {
203  if
204  (
205  !meshFaceUsed[meshFaceI]
206  && localFaces[meshFaceI].size() == curModelFace.size()
207  )
208  {
209  // A possible match. A mesh face will be rotated, so make a copy
210  labelList meshFaceLabels = localFaces[meshFaceI];
211 
212  for
213  (
214  label rotation = 0;
215  rotation < meshFaceLabels.size();
216  rotation++
217  )
218  {
219  // try matching the face
220  label nMatchedLabels = 0;
221 
222  forAll(meshFaceLabels, pointI)
223  {
224  if
225  (
226  pointLabels[curModelFace[pointI]]
227  == meshFaceLabels[pointI]
228  )
229  {
230  nMatchedLabels++;
231  }
232  }
233 
234  if (nMatchedLabels >= 2)
235  {
236  // match!
237  found = true;
238  }
239 
240  if (found)
241  {
242  // match found. Insert mesh face
243  forAll(meshFaceLabels, pointI)
244  {
245  pointLabels[curModelFace[pointI]] =
246  meshFaceLabels[pointI];
247  }
248 
249  meshFaceUsed[meshFaceI] = true;
250 
251  break;
252  }
253  else
254  {
255  // No match found. Rotate face
256  label firstLabel = meshFaceLabels[0];
257 
258  for (label i = 1; i < meshFaceLabels.size(); i++)
259  {
260  meshFaceLabels[i - 1] = meshFaceLabels[i];
261  }
262 
263  meshFaceLabels.last() = firstLabel;
264  }
265  }
266 
267  if (found) break;
268  }
269  }
270 
271  if (!found)
272  {
273  // A model face is not matched. Shape detection failed
275  (
276  "create3DCellShape(const label cellIndex, "
277  "const labelList& faceLabels, const labelListList& faces, "
278  "const labelList& owner, const labelList& neighbour, "
279  "const label fluentCellModelID)"
280  ) << "Cannot find match for face "
281  << modelFaceI
282  << ".\nModel: " << curModel.name() << " model face: "
283  << curModelFace << " Mesh faces: " << localFaces
284  << "Matched points: " << pointLabels
285  << abort(FatalError);
286  }
287  }
288 
289  return cellShape(curModel, pointLabels);
290 }
291 
292 
293 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
294 
295 } // End namespace Foam
296 
297 // ************************************************************************* //
labelList pointLabels(nPoints,-1)
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:59
List< face > faceList
Definition: faceListFwd.H:43
Namespace for OpenFOAM.
void setSize(const label)
Reset size of List.
Definition: List.C:318
cellShape create3DCellShape(const label cellIndex, const labelList &faceLabels, const faceList &faces, const labelList &owner, const labelList &neighbour, const label fluentCellModelID)
#define forAll(list, i)
Definition: UList.H:421
errorManip< error > abort(error &err)
Definition: errorManip.H:131
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
#define forAllReverse(list, i)
Definition: UList.H:424
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const cellModel * lookup(const word &)
Look up a model by name and return a pointer to the model or NULL.
Definition: cellModeller.C:91
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57