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