faceCollapser.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-2012 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 "faceCollapser.H"
27 #include "polyMesh.H"
28 #include "polyTopoChange.H"
29 #include "polyModifyPoint.H"
30 #include "polyModifyFace.H"
31 #include "polyRemoveFace.H"
32 #include "SortableList.H"
33 #include "meshTools.H"
34 #include "OFstream.H"
35 
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 // Insert labelList into labelHashSet. Optional excluded element.
40 void Foam::faceCollapser::insert
41 (
42  const labelList& elems,
43  const label excludeElem,
44  labelHashSet& set
45 )
46 {
47  forAll(elems, i)
48  {
49  if (elems[i] != excludeElem)
50  {
51  set.insert(elems[i]);
52  }
53  }
54 }
55 
56 
57 // Find edge amongst candidate edges. FatalError if none.
58 Foam::label Foam::faceCollapser::findEdge
59 (
60  const edgeList& edges,
61  const labelList& edgeLabels,
62  const label v0,
63  const label v1
64 )
65 {
66  forAll(edgeLabels, i)
67  {
68  label edgeI = edgeLabels[i];
69 
70  const edge& e = edges[edgeI];
71 
72  if
73  (
74  (e[0] == v0 && e[1] == v1)
75  || (e[0] == v1 && e[1] == v0)
76  )
77  {
78  return edgeI;
79  }
80  }
81 
82  FatalErrorIn("findEdge") << "Cannot find edge between vertices " << v0
83  << " and " << v1 << " in edge labels " << edgeLabels
84  << abort(FatalError);
85 
86  return -1;
87 }
88 
89 
90 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
91 
92 // Replace vertices in face
93 void Foam::faceCollapser::filterFace
94 (
95  const Map<labelList>& splitEdges,
96  const label faceI,
97  polyTopoChange& meshMod
98 ) const
99 {
100  const face& f = mesh_.faces()[faceI];
101  const labelList& fEdges = mesh_.faceEdges()[faceI];
102 
103  // Space for replaced vertices and split edges.
104  DynamicList<label> newFace(10 * f.size());
105 
106  forAll(f, fp)
107  {
108  label v0 = f[fp];
109 
110  newFace.append(v0);
111 
112  // Look ahead one to get edge.
113  label fp1 = f.fcIndex(fp);
114 
115  label v1 = f[fp1];
116 
117  // Get split on edge if any.
118  label edgeI = findEdge(mesh_.edges(), fEdges, v0, v1);
119 
121  splitEdges.find(edgeI);
122 
123  if (edgeFnd != splitEdges.end())
124  {
125  // edgeI has been split (by introducing new vertices).
126  // Insert new vertices in face in correct order
127  // (splitEdges was constructed to be from edge.start() to end())
128 
129  const labelList& extraVerts = edgeFnd();
130 
131  if (v0 == mesh_.edges()[edgeI].start())
132  {
133  forAll(extraVerts, i)
134  {
135  newFace.append(extraVerts[i]);
136  }
137  }
138  else
139  {
140  forAllReverse(extraVerts, i)
141  {
142  newFace.append(extraVerts[i]);
143  }
144  }
145  }
146  }
147  face newF(newFace.shrink());
148 
149  //Pout<< "Modifying face:" << faceI << " from " << f << " to " << newFace
150  // << endl;
151 
152  if (newF != f)
153  {
154  label nei = -1;
155 
156  label patchI = -1;
157 
158  if (mesh_.isInternalFace(faceI))
159  {
160  nei = mesh_.faceNeighbour()[faceI];
161  }
162  else
163  {
164  patchI = mesh_.boundaryMesh().whichPatch(faceI);
165  }
166 
167  // Get current zone info
168  label zoneID = mesh_.faceZones().whichZone(faceI);
169 
170  bool zoneFlip = false;
171 
172  if (zoneID >= 0)
173  {
174  const faceZone& fZone = mesh_.faceZones()[zoneID];
175 
176  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
177  }
178 
179  meshMod.setAction
180  (
181  polyModifyFace
182  (
183  newF, // modified face
184  faceI, // label of face being modified
185  mesh_.faceOwner()[faceI], // owner
186  nei, // neighbour
187  false, // face flip
188  patchI, // patch for face
189  false, // remove from zone
190  zoneID, // zone for face
191  zoneFlip // face flip in zone
192  )
193  );
194  }
195 }
196 
197 
198 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
199 
200 // Construct from mesh
201 Foam::faceCollapser::faceCollapser(const polyMesh& mesh)
202 :
203  mesh_(mesh)
204 {}
205 
206 
207 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
208 
210 (
211  const labelList& faceLabels,
212  const labelList& fpStart,
213  const labelList& fpEnd,
214  polyTopoChange& meshMod
215 ) const
216 {
217  const pointField& points = mesh_.points();
218  const edgeList& edges = mesh_.edges();
219  const faceList& faces = mesh_.faces();
220  const labelListList& edgeFaces = mesh_.edgeFaces();
221 
222 
223  // From split edge to newly introduced point(s). Can be more than one per
224  // edge!
225  Map<labelList> splitEdges(faceLabels.size());
226 
227  // Mark faces in any way affect by modifying points/edges. Used later on
228  // to prevent having to redo all faces.
229  labelHashSet affectedFaces(4*faceLabels.size());
230 
231 
232  //
233  // Add/remove vertices and construct mapping
234  //
235 
236  forAll(faceLabels, i)
237  {
238  const label faceI = faceLabels[i];
239 
240  const face& f = faces[faceI];
241 
242  const label fpA = fpStart[i];
243  const label fpB = fpEnd[i];
244 
245  const point& pA = points[f[fpA]];
246  const point& pB = points[f[fpB]];
247 
248  Pout<< "Face:" << f << " collapsed to fp:" << fpA << ' ' << fpB
249  << " with points:" << pA << ' ' << pB
250  << endl;
251 
252  // Create line from fpA to fpB
253  linePointRef lineAB(pA, pB);
254 
255  // Get projections of all vertices onto line.
256 
257  // Distance(squared) to pA for every point on face.
258  SortableList<scalar> dist(f.size());
259 
260  dist[fpA] = 0;
261  dist[fpB] = magSqr(pB - pA);
262 
263  // Step from fpA to fpB
264  // ~~~~~~~~~~~~~~~~~~~~
265  // (by incrementing)
266 
267  label fpMin1 = fpA;
268  label fp = f.fcIndex(fpMin1);
269 
270  while (fp != fpB)
271  {
272  // See where fp sorts. Make sure it is above fpMin1!
273  pointHit near = lineAB.nearestDist(points[f[fp]]);
274 
275  scalar w = magSqr(near.rawPoint() - pA);
276 
277  if (w <= dist[fpMin1])
278  {
279  // Offset.
280  w = dist[fpMin1] + 1e-6*(dist[fpB] - dist[fpA]);
281 
282  point newPoint
283  (
284  pA + Foam::sqrt(w / (dist[fpB] - dist[fpA]))*(pB - pA)
285  );
286 
287  Pout<< "Adapting position of vertex " << f[fp] << " on face "
288  << f << " from " << near.rawPoint() << " to " << newPoint
289  << endl;
290 
291  near.setPoint(newPoint);
292  }
293 
294  // Responsability of caller to make sure polyModifyPoint is only
295  // called once per point. (so max only one collapse face per
296  // edge)
297  meshMod.setAction
298  (
300  (
301  f[fp],
302  near.rawPoint(),
303  false,
304  -1,
305  true
306  )
307  );
308 
309  dist[fp] = w;
310 
311  // Step to next vertex.
312  fpMin1 = fp;
313  fp = f.fcIndex(fpMin1);
314  }
315 
316  // Step from fpA to fpB
317  // ~~~~~~~~~~~~~~~~~~~~
318  // (by decrementing)
319 
320  fpMin1 = fpA;
321  fp = f.rcIndex(fpMin1);
322 
323  while (fp != fpB)
324  {
325  // See where fp sorts. Make sure it is below fpMin1!
326  pointHit near = lineAB.nearestDist(points[f[fp]]);
327 
328  scalar w = magSqr(near.rawPoint() - pA);
329 
330  if (w <= dist[fpMin1])
331  {
332  // Offset.
333  w = dist[fpMin1] + 1e-6*(dist[fpB] - dist[fpA]);
334 
335  point newPoint
336  (
337  pA + Foam::sqrt(w / (dist[fpB] - dist[fpA]))*(pB - pA)
338  );
339 
340  Pout<< "Adapting position of vertex " << f[fp] << " on face "
341  << f << " from " << near.rawPoint() << " to " << newPoint
342  << endl;
343 
344  near.setPoint(newPoint);
345  }
346 
347  // Responsability of caller to make sure polyModifyPoint is only
348  // called once per point. (so max only one collapse face per
349  // edge)
350  meshMod.setAction
351  (
353  (
354  f[fp],
355  near.rawPoint(),
356  false,
357  -1,
358  true
359  )
360  );
361 
362  dist[fp] = w;
363 
364  // Step to previous vertex.
365  fpMin1 = fp;
366  fp = f.rcIndex(fpMin1);
367  }
368 
369  dist.sort();
370 
371  // Check that fpB sorts latest.
372  if (dist.indices()[dist.size()-1] != fpB)
373  {
374  OFstream str("conflictingFace.obj");
375  meshTools::writeOBJ(str, faceList(1, f), points);
376 
377  FatalErrorIn("faceCollapser::setRefinement")
378  << "Trying to collapse face:" << faceI << " vertices:" << f
379  << " to edges between vertices " << f[fpA] << " and "
380  << f[fpB] << " but " << f[fpB] << " does not seem to be the"
381  << " vertex furthest away from " << f[fpA] << endl
382  << "Dumped conflicting face to obj file conflictingFace.obj"
383  << abort(FatalError);
384  }
385 
386 
387  // From fp to index in sort:
388  Pout<< "Face:" << f << " fpA:" << fpA << " fpB:" << fpB << nl;
389 
390  labelList sortedFp(f.size());
391  forAll(dist.indices(), i)
392  {
393  label fp = dist.indices()[i];
394 
395  Pout<< " fp:" << fp << " distance:" << dist[i] << nl;
396 
397  sortedFp[fp] = i;
398  }
399 
400  const labelList& fEdges = mesh_.faceEdges()[faceI];
401 
402  // Now look up all edges in the face and see if they get extra
403  // vertices inserted and build an edge-to-intersected-points table.
404 
405  // Order of inserted points is in edge order (from e.start to
406  // e.end)
407 
408  forAll(f, fp)
409  {
410  label fp1 = f.fcIndex(fp);
411 
412  // Get index in sorted list
413  label sorted0 = sortedFp[fp];
414  label sorted1 = sortedFp[fp1];
415 
416  // Get indices in increasing order
417  DynamicList<label> edgePoints(f.size());
418 
419  // Whether edgePoints are from fp to fp1
420  bool fpToFp1;
421 
422  if (sorted0 < sorted1)
423  {
424  fpToFp1 = true;
425 
426  for (label j = sorted0+1; j < sorted1; j++)
427  {
428  edgePoints.append(f[dist.indices()[j]]);
429  }
430  }
431  else
432  {
433  fpToFp1 = false;
434 
435  for (label j = sorted1+1; j < sorted0; j++)
436  {
437  edgePoints.append(f[dist.indices()[j]]);
438  }
439  }
440 
441  if (edgePoints.size())
442  {
443  edgePoints.shrink();
444 
445  label edgeI = findEdge(edges, fEdges, f[fp], f[fp1]);
446 
447  const edge& e = edges[edgeI];
448 
449  if (fpToFp1 == (f[fp] == e.start()))
450  {
451  splitEdges.insert(edgeI, edgePoints);
452  }
453  else
454  {
455  reverse(edgePoints);
456  splitEdges.insert(edgeI, edgePoints);
457  }
458 
459  // Mark all faces affected
460  insert(edgeFaces[edgeI], faceI, affectedFaces);
461  }
462  }
463  }
464 
465  forAllConstIter(Map<labelList>, splitEdges, iter)
466  {
467  Pout<< "Split edge:" << iter.key()
468  << " verts:" << mesh_.edges()[iter.key()]
469  << " in:" << nl;
470  const labelList& edgePoints = iter();
471 
472  forAll(edgePoints, i)
473  {
474  Pout<< " " << edgePoints[i] << nl;
475  }
476  }
477 
478 
479  //
480  // Remove faces.
481  //
482 
483  forAll(faceLabels, i)
484  {
485  const label faceI = faceLabels[i];
486 
487  meshMod.setAction(polyRemoveFace(faceI));
488 
489  // Update list of faces we still have to modify
490  affectedFaces.erase(faceI);
491  }
492 
493 
494  //
495  // Modify faces affected (but not removed)
496  //
497 
498  forAllConstIter(labelHashSet, affectedFaces, iter)
499  {
500  filterFace(splitEdges, iter.key(), meshMod);
501  }
502 }
503 
504 
505 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
Output to file stream.
Definition: OFstream.H:81
PointHit< Point > nearestDist(const Point &p) const
Return nearest distance to line from a given point.
Definition: lineI.H:95
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
const pointField & points
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
const labelList & indices() const
Return the list of sorted indices. Updated every sort.
Definition: SortableList.H:90
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:232
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
void setPoint(const Point &p)
Definition: PointHit.H:181
void setRefinement(const labelList &faceLabels, const labelList &fpA, const labelList &fpB, polyTopoChange &) const
Collapse faces along endpoints. Play commands into.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49
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
Class describing modification of a point.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
List< face > faceList
Definition: faceListFwd.H:43
dynamicFvMesh & mesh
void sort()
(stable) sort the list (if changed after construction time)
Definition: SortableList.C:95
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
const labelListList & faceEdges() const
Class containing data for face removal.
A line primitive.
Definition: line.H:56
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
A list that is sorted upon construction or when explicitly requested with the sort() method...
Definition: List.H:65
#define forAll(list, i)
Definition: UList.H:421
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
List< edge > edgeList
Definition: edgeList.H:38
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
#define forAllReverse(list, i)
Definition: UList.H:424
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
error FatalError
const labelListList & edgeFaces() const
List< label > labelList
A List of labels.
Definition: labelList.H:56
void reverse(UList< T > &, const label n)
Definition: UListI.H:322
Direct mesh changes based on v1.3 polyTopoChange syntax.
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
label start() const
Return start vertex label.
Definition: edgeI.H:81
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1060
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:139
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1079
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:209
HashTable< T, label, Hash< label > >::const_iterator const_iterator
Definition: Map.H:59
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.