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-2016 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  << "Number of face labels not equal to"
79  << "number of face in the model. "
80  << "Number of face labels: " << faceLabels.size()
81  << " number of faces in model: " << curModel.nFaces()
82  << abort(FatalError);
83  }
84 
85  // make a list of outward-pointing faces
86  labelListList localFaces(faceLabels.size());
87 
88  forAll(faceLabels, facei)
89  {
90  const label curFaceLabel = faceLabels[facei];
91 
92  const labelList& curFace = faces[curFaceLabel];
93 
94  if (owner[curFaceLabel] == cellIndex)
95  {
96  localFaces[facei] = curFace;
97  }
98  else if (neighbour[curFaceLabel] == cellIndex)
99  {
100  // Reverse the face
101  localFaces[facei].setSize(curFace.size());
102 
103  forAllReverse(curFace, i)
104  {
105  localFaces[facei][curFace.size() - i - 1] =
106  curFace[i];
107  }
108  }
109  else
110  {
112  << "face " << curFaceLabel
113  << " does not belong to cell " << cellIndex
114  << ". Face owner: " << owner[curFaceLabel] << " neighbour: "
115  << neighbour[curFaceLabel]
116  << abort(FatalError);
117  }
118  }
119 
120  // Algorithm:
121  // Make an empty list of pointLabels and initialise it with -1. Pick the
122  // first face from modelFaces and look through the faces to find one with
123  // the same number of labels. Insert face by copying its labels into
124  // pointLabels. Mark the face as used. Loop through all model faces.
125  // For each model face loop through faces. If the face is unused and the
126  // numbers of labels fit, try to match the face onto the point labels. If
127  // at least one edge is matched, insert the face into pointLabels. If at
128  // any stage the matching algorithm reaches the end of faces, the matching
129  // algorithm has failed. Once all the faces are matched, the list of
130  // pointLabels defines the model.
131 
132  // Make a list of empty pointLabels
133  labelList pointLabels(curModel.nPoints(), -1);
134 
135  // Follow the used mesh faces
136  List<bool> meshFaceUsed(localFaces.size(), false);
137 
138  // Get the raw model faces
139  const faceList& modelFaces = curModel.modelFaces();
140 
141  // Insert the first face into the list
142  const labelList& firstModelFace =
143  modelFaces[faceMatchingOrder[fluentCellModelID][0]];
144 
145  bool found = false;
146 
147  forAll(localFaces, meshFacei)
148  {
149  if (localFaces[meshFacei].size() == firstModelFace.size())
150  {
151  // Match. Insert points into the pointLabels
152  found = true;
153 
154  const labelList& curMeshFace = localFaces[meshFacei];
155 
156  meshFaceUsed[meshFacei] = true;
157 
158  forAll(curMeshFace, pointi)
159  {
160  pointLabels[firstModelFace[pointi]] = curMeshFace[pointi];
161  }
162 
163  break;
164  }
165  }
166 
167  if (!found)
168  {
170  << "Cannot find match for first face. "
171  << "cell model: " << curModel.name() << " first model face: "
172  << firstModelFace << " Mesh faces: " << localFaces
173  << abort(FatalError);
174  }
175 
176  for (label modelFacei = 1; modelFacei < modelFaces.size(); modelFacei++)
177  {
178  // get the next model face
179  const labelList& curModelFace =
180  modelFaces
181  [faceMatchingOrder[fluentCellModelID][modelFacei]];
182 
183  found = false;
184 
185  // Loop through mesh faces until a match is found
186  forAll(localFaces, meshFacei)
187  {
188  if
189  (
190  !meshFaceUsed[meshFacei]
191  && localFaces[meshFacei].size() == curModelFace.size()
192  )
193  {
194  // A possible match. A mesh face will be rotated, so make a copy
195  labelList meshFaceLabels = localFaces[meshFacei];
196 
197  for
198  (
199  label rotation = 0;
200  rotation < meshFaceLabels.size();
201  rotation++
202  )
203  {
204  // try matching the face
205  label nMatchedLabels = 0;
206 
207  forAll(meshFaceLabels, pointi)
208  {
209  if
210  (
211  pointLabels[curModelFace[pointi]]
212  == meshFaceLabels[pointi]
213  )
214  {
215  nMatchedLabels++;
216  }
217  }
218 
219  if (nMatchedLabels >= 2)
220  {
221  // match!
222  found = true;
223  }
224 
225  if (found)
226  {
227  // match found. Insert mesh face
228  forAll(meshFaceLabels, pointi)
229  {
230  pointLabels[curModelFace[pointi]] =
231  meshFaceLabels[pointi];
232  }
233 
234  meshFaceUsed[meshFacei] = true;
235 
236  break;
237  }
238  else
239  {
240  // No match found. Rotate face
241  label firstLabel = meshFaceLabels[0];
242 
243  for (label i = 1; i < meshFaceLabels.size(); i++)
244  {
245  meshFaceLabels[i - 1] = meshFaceLabels[i];
246  }
247 
248  meshFaceLabels.last() = firstLabel;
249  }
250  }
251 
252  if (found) break;
253  }
254  }
255 
256  if (!found)
257  {
258  // A model face is not matched. Shape detection failed
260  << "Cannot find match for face "
261  << modelFacei
262  << ".\nModel: " << curModel.name() << " model face: "
263  << curModelFace << " Mesh faces: " << localFaces
264  << "Matched points: " << pointLabels
265  << abort(FatalError);
266  }
267  }
268 
269  return cellShape(curModel, pointLabels);
270 }
271 
272 
273 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
274 
275 } // End namespace Foam
276 
277 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
static const cellModel * lookup(const word &)
Look up a model by name and return a pointer to the model or NULL.
Definition: cellModeller.C:88
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
labelList pointLabels(nPoints,-1)
List< face > faceList
Definition: faceListFwd.H:43
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: UList.H:440
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void setSize(const label)
Reset size of List.
Definition: List.C:295
cellShape create3DCellShape(const label cellIndex, const labelList &faceLabels, const faceList &faces, const labelList &owner, const labelList &neighbour, const label fluentCellModelID)
Namespace for OpenFOAM.