extendedEdgeMeshTemplates.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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 \*---------------------------------------------------------------------------*/
25 
26 #include "extendedEdgeMesh.H"
27 #include "ListListOps.H"
28 #include "unitConversion.H"
29 #include "PackedBoolList.H"
30 #include "PatchTools.H"
31 #include "searchableBox.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 template<class Patch>
37 (
38  const Patch& surf,
39  const labelList& featureEdges,
40  const labelList& regionFeatureEdges,// subset of featureEdges: inter-region
41  const labelList& featurePoints
42 )
43 {
44  const pointField& sFeatLocalPts(surf.localPoints());
45  const edgeList& sFeatEds(surf.edges());
46  const labelListList edgeFaces = PatchTools::sortedEdgeFaces(surf);
47  const vectorField& faceNormals = surf.faceNormals();
48  const labelListList pointEdges = PatchTools::sortedPointEdges(surf);
49 
50  // Extract and reorder the data from surfaceFeatures
51 
52  // References to the surfaceFeatures data
53 
54  // Filling the extendedEdgeMesh with the raw geometrical data.
55 
56  label nFeatEds = featureEdges.size();
57  label nFeatPts = featurePoints.size();
58 
59  DynamicList<point> tmpPts;
60  edgeList eds(nFeatEds);
61  DynamicList<vector> norms;
62  vectorField edgeDirections(nFeatEds);
63  labelListList edgeNormals(nFeatEds);
64  labelListList normalDirections(nFeatEds);
65  DynamicList<label> regionEdges;
66 
67  // Keep track of the ordered feature point feature edges
68  labelListList featurePointFeatureEdges(nFeatPts);
69  forAll(featurePointFeatureEdges, pI)
70  {
71  featurePointFeatureEdges[pI] = pointEdges[featurePoints[pI]];
72  }
73 
74  // Mapping between old and new indices, there is entry in the map for each
75  // of surf.localPoints, -1 means that this point hasn't been used (yet),
76  // >= 0 corresponds to the index
77  labelList pointMap(sFeatLocalPts.size(), -1);
78 
79  // Mapping between surface edge index and its feature edge index. -1 if it
80  // is not a feature edge
81  labelList edgeMap(sFeatEds.size(), -1);
82 
83  // Noting when the normal of a face has been used so not to duplicate
84  labelList faceMap(surf.size(), -1);
85 
86  // Collecting the status of edge for subsequent sorting
87  List<edgeStatus> edStatus(nFeatEds, NONE);
88 
89  forAll(featurePoints, i)
90  {
91  label sFPI = featurePoints[i];
92 
93  tmpPts.append(sFeatLocalPts[sFPI]);
94 
95  pointMap[sFPI] = tmpPts.size() - 1;
96  }
97 
98  // All feature points have been added
99  nonFeatureStart_ = tmpPts.size();
100 
101  PackedBoolList isRegionFeatureEdge(regionFeatureEdges);
102 
103  forAll(featureEdges, i)
104  {
105  label sFEI = featureEdges[i];
106 
107  edgeMap[sFEI] = i;
108 
109  const edge& fE = sFeatEds[sFEI];
110 
111  edgeDirections[i] = fE.vec(sFeatLocalPts);
112 
113  // Check to see if the points have been already used
114  if (pointMap[fE.start()] == -1)
115  {
116  tmpPts.append(sFeatLocalPts[fE.start()]);
117 
118  pointMap[fE.start()] = tmpPts.size() - 1;
119  }
120 
121  eds[i].start() = pointMap[fE.start()];
122 
123  if (pointMap[fE.end()] == -1)
124  {
125  tmpPts.append(sFeatLocalPts[fE.end()]);
126 
127  pointMap[fE.end()] = tmpPts.size() - 1;
128  }
129 
130  eds[i].end() = pointMap[fE.end()];
131 
132  // Pick up the faces adjacent to the feature edge
133  const labelList& eFaces = edgeFaces[sFEI];
134 
135  edgeNormals[i].setSize(eFaces.size());
136  normalDirections[i].setSize(eFaces.size());
137 
138  forAll(eFaces, j)
139  {
140  label eFI = eFaces[j];
141 
142  // Check to see if the points have been already used
143  if (faceMap[eFI] == -1)
144  {
145  norms.append(faceNormals[eFI]);
146 
147  faceMap[eFI] = norms.size() - 1;
148  }
149 
150  edgeNormals[i][j] = faceMap[eFI];
151 
152  const vector cross = (faceNormals[eFI] ^ edgeDirections[i]);
153  const vector fC0tofE0 =
154  surf[eFI].centre(surf.points())
155  - sFeatLocalPts[fE.start()];
156 
157  normalDirections[i][j] =
158  (
159  (
160  (cross/(mag(cross) + vSmall))
161  & (fC0tofE0/(mag(fC0tofE0)+ vSmall))
162  )
163  > 0.0
164  ? 1
165  : -1
166  );
167  }
168 
169  vector fC0tofC1(Zero);
170 
171  if (eFaces.size() == 2)
172  {
173  fC0tofC1 =
174  surf[eFaces[1]].centre(surf.points())
175  - surf[eFaces[0]].centre(surf.points());
176  }
177 
178  edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
179 
180  if (isRegionFeatureEdge[i])
181  {
182  regionEdges.append(i);
183  }
184  }
185 
186  // Populate feature point feature edges
187  DynamicList<label> newFeatureEdges;
188 
189  forAll(featurePointFeatureEdges, pI)
190  {
191  const labelList& fpfe = featurePointFeatureEdges[pI];
192 
193  newFeatureEdges.setCapacity(fpfe.size());
194 
195  forAll(fpfe, eI)
196  {
197  const label oldEdgeIndex = fpfe[eI];
198 
199  const label newFeatureEdgeIndex = edgeMap[oldEdgeIndex];
200 
201  if (newFeatureEdgeIndex != -1)
202  {
203  newFeatureEdges.append(newFeatureEdgeIndex);
204  }
205  }
206 
207  featurePointFeatureEdges[pI].transfer(newFeatureEdges);
208  }
209 
210  // Reorder the edges by classification
211  List<DynamicList<label>> allEds(nEdgeTypes);
212 
213  DynamicList<label>& externalEds(allEds[0]);
214  DynamicList<label>& internalEds(allEds[1]);
215  DynamicList<label>& flatEds(allEds[2]);
216  DynamicList<label>& openEds(allEds[3]);
217  DynamicList<label>& multipleEds(allEds[4]);
218 
219  forAll(eds, i)
220  {
221  edgeStatus eStat = edStatus[i];
222 
223  if (eStat == EXTERNAL)
224  {
225  externalEds.append(i);
226  }
227  else if (eStat == INTERNAL)
228  {
229  internalEds.append(i);
230  }
231  else if (eStat == FLAT)
232  {
233  flatEds.append(i);
234  }
235  else if (eStat == OPEN)
236  {
237  openEds.append(i);
238  }
239  else if (eStat == MULTIPLE)
240  {
241  multipleEds.append(i);
242  }
243  else if (eStat == NONE)
244  {
246  << nl << "classifyEdge returned NONE on edge "
247  << eds[i]
248  << ". There is a problem with definition of this edge."
249  << nl << abort(FatalError);
250  }
251  }
252 
253  internalStart_ = externalEds.size();
254  flatStart_ = internalStart_ + internalEds.size();
255  openStart_ = flatStart_ + flatEds.size();
256  multipleStart_ = openStart_ + openEds.size();
257 
258  labelList edMap
259  (
260  ListListOps::combine<labelList>
261  (
262  allEds,
264  )
265  );
266 
267  edMap = invert(edMap.size(), edMap);
268 
269  inplaceReorder(edMap, eds);
270  inplaceReorder(edMap, edStatus);
271  inplaceReorder(edMap, edgeDirections);
272  inplaceReorder(edMap, edgeNormals);
273  inplaceReorder(edMap, normalDirections);
274  inplaceRenumber(edMap, regionEdges);
275 
276  forAll(featurePointFeatureEdges, pI)
277  {
278  inplaceRenumber(edMap, featurePointFeatureEdges[pI]);
279  }
280 
281  pointField pts(tmpPts);
282 
283  // Initialise the edgeMesh
284  edgeMesh::operator=(edgeMesh(pts, eds));
285 
286  // Initialise sorted edge related data
287  edgeDirections_ = edgeDirections/(mag(edgeDirections) + vSmall);
288  edgeNormals_ = edgeNormals;
289  normalDirections_ = normalDirections;
290  regionEdges_ = regionEdges;
291 
292  // Normals are all now found and indirectly addressed, can also be stored
293  normals_ = vectorField(norms);
294 
295 
296  // Reorder the feature points by classification
297 
298  List<DynamicList<label>> allPts(3);
299 
300  DynamicList<label>& convexPts(allPts[0]);
301  DynamicList<label>& concavePts(allPts[1]);
302  DynamicList<label>& mixedPts(allPts[2]);
303 
304  for (label i = 0; i < nonFeatureStart_; i++)
305  {
306  pointStatus ptStatus = classifyFeaturePoint(i);
307 
308  if (ptStatus == CONVEX)
309  {
310  convexPts.append(i);
311  }
312  else if (ptStatus == CONCAVE)
313  {
314  concavePts.append(i);
315  }
316  else if (ptStatus == MIXED)
317  {
318  mixedPts.append(i);
319  }
320  else if (ptStatus == NONFEATURE)
321  {
323  << nl << "classifyFeaturePoint returned NONFEATURE on point at "
324  << points()[i]
325  << ". There is a problem with definition of this feature point."
326  << nl << abort(FatalError);
327  }
328  }
329 
330  concaveStart_ = convexPts.size();
331  mixedStart_ = concaveStart_ + concavePts.size();
332 
333  labelList ftPtMap
334  (
335  ListListOps::combine<labelList>
336  (
337  allPts,
339  )
340  );
341 
342  ftPtMap = invert(ftPtMap.size(), ftPtMap);
343 
344  // Creating the ptMap from the ftPtMap with identity values up to the size
345  // of pts to create an oldToNew map for inplaceReorder
346 
347  labelList ptMap(identity(pts.size()));
348 
349  forAll(ftPtMap, i)
350  {
351  ptMap[i] = ftPtMap[i];
352  }
353 
354  inplaceReorder(ptMap, pts);
355  inplaceReorder(ptMap, featurePointFeatureEdges);
356 
357  forAll(eds, i)
358  {
359  inplaceRenumber(ptMap, eds[i]);
360  }
361 
362  // Reinitialise the edgeMesh with sorted feature points and
363  // renumbered edges
364  reset(move(pts), move(eds));
365 
366  // Generate the featurePointNormals
367 
368  labelListList featurePointNormals(nonFeatureStart_);
369 
370  for (label i = 0; i < nonFeatureStart_; i++)
371  {
372  DynamicList<label> tmpFtPtNorms;
373 
374  const labelList& ptEds = edgeMesh::pointEdges()[i];
375 
376  forAll(ptEds, j)
377  {
378  const labelList& ptEdNorms(edgeNormals[ptEds[j]]);
379 
380  forAll(ptEdNorms, k)
381  {
382  if (findIndex(tmpFtPtNorms, ptEdNorms[k]) == -1)
383  {
384  bool addNormal = true;
385 
386  // Check that the normal direction is unique at this feature
387  forAll(tmpFtPtNorms, q)
388  {
389  if
390  (
391  (normals_[ptEdNorms[k]] & normals_[tmpFtPtNorms[q]])
392  > cosNormalAngleTol_
393  )
394  {
395  // Parallel to an existing normal, do not add
396  addNormal = false;
397 
398  break;
399  }
400  }
401 
402  if (addNormal)
403  {
404  tmpFtPtNorms.append(ptEdNorms[k]);
405  }
406  }
407  }
408  }
409 
410  featurePointNormals[i] = tmpFtPtNorms;
411  }
412 
413  featurePointNormals_ = featurePointNormals;
414  featurePointEdges_ = featurePointFeatureEdges;
415 }
416 
417 
418 // ************************************************************************* //
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
void sortPointsAndEdges(const Patch &, const labelList &featureEdges, const labelList &regionFeatureEdges, const labelList &feaurePoints)
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
volVectorField vectorField(fieldObject, mesh)
void cross(FieldField< Field1, typename crossProduct< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
label k
Boltzmann constant.
iterator end()
Return an iterator to end traversing the UList.
Definition: UListI.H:224
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt >> &) const
Return *this (used for point which is a typedef to Vector<scalar>.
Definition: VectorI.H:116
labelList invert(const label len, const labelUList &)
Invert one-to-one map. Unmapped elements will be -1.
Definition: ListOps.C:37
const pointField & points
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
void setCapacity(const label)
Alter the size of the underlying storage.
Definition: DynamicListI.H:130
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
vector vec(const pointField &) const
Return the vector (end - start)
Definition: edgeI.H:175
static const char nl
Definition: Ostream.H:260
Points connected by edges.
Definition: edgeMesh.H:69
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void setSize(const label)
Reset size of List.
Definition: List.C:281
A bit-packed bool list.
label end() const
Return end vertex label.
Definition: edgeI.H:92
dimensioned< scalar > mag(const dimensioned< Type > &)
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
label start() const
Return start vertex label.
Definition: edgeI.H:81