FreeFOAM The Cross-Platform CFD Toolkit
edgeCollapser.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 \*---------------------------------------------------------------------------*/
25 
26 #include "edgeCollapser.H"
27 #include <OpenFOAM/polyMesh.H>
28 #include "polyTopoChange.H"
29 #include <OpenFOAM/ListOps.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 Foam::label Foam::edgeCollapser::findIndex
34 (
35  const labelList& elems,
36  const label start,
37  const label size,
38  const label val
39 )
40 {
41  for (label i = start; i < size; i++)
42  {
43  if (elems[i] == val)
44  {
45  return i;
46  }
47  }
48  return -1;
49 }
50 
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
54 // Changes region of connected set of points
55 Foam::label Foam::edgeCollapser::changePointRegion
56 (
57  const label pointI,
58  const label oldRegion,
59  const label newRegion
60 )
61 {
62  label nChanged = 0;
63 
64  if (pointRegion_[pointI] == oldRegion)
65  {
66  pointRegion_[pointI] = newRegion;
67  nChanged++;
68 
69  // Step to neighbouring point across edges with same region number
70 
71  const labelList& pEdges = mesh_.pointEdges()[pointI];
72 
73  forAll(pEdges, i)
74  {
75  label otherPointI = mesh_.edges()[pEdges[i]].otherVertex(pointI);
76 
77  nChanged += changePointRegion(otherPointI, oldRegion, newRegion);
78  }
79  }
80  return nChanged;
81 }
82 
83 
84 bool Foam::edgeCollapser::pointRemoved(const label pointI) const
85 {
86  label region = pointRegion_[pointI];
87 
88  if (region == -1 || pointRegionMaster_[region] == pointI)
89  {
90  return false;
91  }
92  else
93  {
94  return true;
95  }
96 }
97 
98 
99 void Foam::edgeCollapser::filterFace(const label faceI, face& f) const
100 {
101  label newFp = 0;
102 
103  forAll(f, fp)
104  {
105  label pointI = f[fp];
106 
107  label region = pointRegion_[pointI];
108 
109  if (region == -1)
110  {
111  f[newFp++] = pointI;
112  }
113  else
114  {
115  label master = pointRegionMaster_[region];
116 
117  if (findIndex(f, 0, newFp, master) == -1)
118  {
119  f[newFp++] = master;
120  }
121  }
122  }
123 
124  // Check for pinched face. Tries to correct
125  // - consecutive duplicate vertex. Removes duplicate vertex.
126  // - duplicate vertex with one other vertex in between (spike).
127  // Both of these should not really occur! and should be checked before
128  // collapsing edges.
129 
130  const label size = newFp;
131 
132  newFp = 2;
133 
134  for (label fp = 2; fp < size; fp++)
135  {
136  label fp1 = fp-1;
137  label fp2 = fp-2;
138 
139  label pointI = f[fp];
140 
141  // Search for previous occurrence.
142  label index = findIndex(f, 0, fp, pointI);
143 
144  if (index == fp1)
145  {
146  WarningIn
147  (
148  "Foam::edgeCollapser::filterFace(const label faceI, "
149  "face& f) const"
150  ) << "Removing consecutive duplicate vertex in face "
151  << f << endl;
152  // Don't store current pointI
153  }
154  else if (index == fp2)
155  {
156  WarningIn
157  (
158  "Foam::edgeCollapser::filterFace(const label faceI, "
159  "face& f) const"
160  ) << "Removing non-consecutive duplicate vertex in face "
161  << f << endl;
162  // Don't store current pointI and remove previous
163  newFp--;
164  }
165  else if (index != -1)
166  {
167  WarningIn
168  (
169  "Foam::edgeCollapser::filterFace(const label faceI, "
170  "face& f) const"
171  ) << "Pinched face " << f << endl;
172  f[newFp++] = pointI;
173  }
174  else
175  {
176  f[newFp++] = pointI;
177  }
178  }
179 
180  f.setSize(newFp);
181 }
182 
183 
184 // Debugging.
185 void Foam::edgeCollapser::printRegions() const
186 {
187  forAll(pointRegionMaster_, regionI)
188  {
189  label master = pointRegionMaster_[regionI];
190 
191  if (master != -1)
192  {
193  Info<< "Region:" << regionI << nl
194  << " master:" << master
195  << ' ' << mesh_.points()[master] << nl;
196 
197  forAll(pointRegion_, pointI)
198  {
199  if (pointRegion_[pointI] == regionI && pointI != master)
200  {
201  Info<< " slave:" << pointI
202  << ' ' << mesh_.points()[pointI] << nl;
203  }
204  }
205  }
206  }
207 }
208 
209 
210 
211 // Collapse list of edges
212 void Foam::edgeCollapser::collapseEdges(const labelList& edgeLabels)
213 {
214  const edgeList& edges = mesh_.edges();
215 
216  forAll(edgeLabels, i)
217  {
218  label edgeI = edgeLabels[i];
219  const edge& e = edges[edgeI];
220 
221  label region0 = pointRegion_[e[0]];
222  label region1 = pointRegion_[e[1]];
223 
224  if (region0 == -1)
225  {
226  if (region1 == -1)
227  {
228  // Both unaffected. Choose ad lib.
229  collapseEdge(edgeI, e[0]);
230  }
231  else
232  {
233  // Collapse to whatever e[1] collapses
234  collapseEdge(edgeI, e[1]);
235  }
236  }
237  else
238  {
239  if (region1 == -1)
240  {
241  // Collapse to whatever e[0] collapses
242  collapseEdge(edgeI, e[0]);
243  }
244  else
245  {
246  // Both collapsed.
247  if (pointRegionMaster_[region0] == e[0])
248  {
249  // e[0] is a master
250  collapseEdge(edgeI, e[0]);
251  }
252  else if (pointRegionMaster_[region1] == e[1])
253  {
254  // e[1] is a master
255  collapseEdge(edgeI, e[1]);
256  }
257  else
258  {
259  // Dont know
260  collapseEdge(edgeI, e[0]);
261  }
262  }
263  }
264  }
265 }
266 
267 
268 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
269 
270 // Construct from mesh
271 Foam::edgeCollapser::edgeCollapser(const polyMesh& mesh)
272 :
273  mesh_(mesh),
274  pointRegion_(mesh.nPoints(), -1),
275  pointRegionMaster_(mesh.nPoints() / 100),
276  freeRegions_()
277 {}
278 
279 
280 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
281 
282 bool Foam::edgeCollapser::unaffectedEdge(const label edgeI) const
283 {
284  const edge& e = mesh_.edges()[edgeI];
285 
286  return (pointRegion_[e[0]] == -1) && (pointRegion_[e[1]] == -1);
287 }
288 
289 
290 bool Foam::edgeCollapser::collapseEdge(const label edgeI, const label master)
291 {
292  const edge& e = mesh_.edges()[edgeI];
293 
294  label pointRegion0 = pointRegion_[e[0]];
295  label pointRegion1 = pointRegion_[e[1]];
296 
297  if (pointRegion0 == -1)
298  {
299  if (pointRegion1 == -1)
300  {
301  // Both endpoints not collapsed. Create new region.
302 
303  label freeRegion = -1;
304 
305  if (freeRegions_.size())
306  {
307  freeRegion = freeRegions_.removeHead();
308 
309  if (pointRegionMaster_[freeRegion] != -1)
310  {
312  ("edgeCollapser::collapseEdge(const label, const label)")
313  << "Problem : freeed region :" << freeRegion
314  << " has already master "
315  << pointRegionMaster_[freeRegion]
316  << abort(FatalError);
317  }
318  }
319  else
320  {
321  // If no region found create one. This is the only place where
322  // new regions are created.
323  freeRegion = pointRegionMaster_.size();
324  }
325 
326  pointRegion_[e[0]] = freeRegion;
327  pointRegion_[e[1]] = freeRegion;
328 
329  pointRegionMaster_(freeRegion) = master;
330  }
331  else
332  {
333  // e[1] is part of collapse network, e[0] not. Add e0 to e1 region.
334  pointRegion_[e[0]] = pointRegion1;
335 
336  pointRegionMaster_[pointRegion1] = master;
337  }
338  }
339  else
340  {
341  if (pointRegion1 == -1)
342  {
343  // e[0] is part of collapse network. Add e1 to e0 region
344  pointRegion_[e[1]] = pointRegion0;
345 
346  pointRegionMaster_[pointRegion0] = master;
347  }
348  else if (pointRegion0 != pointRegion1)
349  {
350  // Both part of collapse network. Merge the two regions.
351 
352  // Use the smaller region number for the whole network.
353  label minRegion = min(pointRegion0, pointRegion1);
354  label maxRegion = max(pointRegion0, pointRegion1);
355 
356  // Use minRegion as region for combined net, free maxRegion.
357  pointRegionMaster_[minRegion] = master;
358  pointRegionMaster_[maxRegion] = -1;
359  freeRegions_.insert(maxRegion);
360 
361  if (minRegion != pointRegion0)
362  {
363  changePointRegion(e[0], pointRegion0, minRegion);
364  }
365  if (minRegion != pointRegion1)
366  {
367  changePointRegion(e[1], pointRegion1, minRegion);
368  }
369  }
370  }
371 
372  return true;
373 }
374 
375 
377 {
378  const cellList& cells = mesh_.cells();
379  const labelList& faceOwner = mesh_.faceOwner();
380  const labelList& faceNeighbour = mesh_.faceNeighbour();
381  const labelListList& pointFaces = mesh_.pointFaces();
382  const labelListList& cellEdges = mesh_.cellEdges();
383 
384  // Print regions:
385  //printRegions()
386 
387  bool meshChanged = false;
388 
389 
390  // Current faces (is also collapseStatus: f.size() < 3)
391  faceList newFaces(mesh_.faces());
392 
393  // Current cellCollapse status
394  boolList cellRemoved(mesh_.nCells(), false);
395 
396 
397  do
398  {
399  // Update face collapse from edge collapses
400  forAll(newFaces, faceI)
401  {
402  filterFace(faceI, newFaces[faceI]);
403  }
404 
405  // Check if faces to be collapsed cause cells to become collapsed.
406  label nCellCollapsed = 0;
407 
408  forAll(cells, cellI)
409  {
410  if (!cellRemoved[cellI])
411  {
412  const cell& cFaces = cells[cellI];
413 
414  label nFaces = cFaces.size();
415 
416  forAll(cFaces, i)
417  {
418  label faceI = cFaces[i];
419 
420  if (newFaces[faceI].size() < 3)
421  {
422  --nFaces;
423 
424  if (nFaces < 4)
425  {
426  Info<< "Cell:" << cellI
427  << " uses faces:" << cFaces
428  << " of which too many are marked for removal:"
429  << endl
430  << " ";
431  forAll(cFaces, j)
432  {
433  if (newFaces[cFaces[j]].size() < 3)
434  {
435  Info<< ' '<< cFaces[j];
436  }
437  }
438  Info<< endl;
439 
440  cellRemoved[cellI] = true;
441 
442  // Collapse all edges of cell to nothing
443  collapseEdges(cellEdges[cellI]);
444 
445  nCellCollapsed++;
446 
447  break;
448  }
449  }
450  }
451  }
452  }
453 
454  if (nCellCollapsed == 0)
455  {
456  break;
457  }
458  } while(true);
459 
460 
461  // Keep track of faces that have been done already.
462  boolList doneFace(mesh_.nFaces(), false);
463 
464  {
465  // Mark points used.
466  boolList usedPoint(mesh_.nPoints(), false);
467 
468 
469  forAll(cellRemoved, cellI)
470  {
471  if (cellRemoved[cellI])
472  {
473  meshMod.removeCell(cellI, -1);
474  }
475  }
476 
477 
478  // Remove faces
479  forAll(newFaces, faceI)
480  {
481  const face& f = newFaces[faceI];
482 
483  if (f.size() < 3)
484  {
485  meshMod.removeFace(faceI, -1);
486  meshChanged = true;
487 
488  // Mark face as been done.
489  doneFace[faceI] = true;
490  }
491  else
492  {
493  // Kept face. Mark vertices
494  forAll(f, fp)
495  {
496  usedPoint[f[fp]] = true;
497  }
498  }
499  }
500 
501  // Remove unused vertices that have not been marked for removal already
502  forAll(usedPoint, pointI)
503  {
504  if (!usedPoint[pointI] && !pointRemoved(pointI))
505  {
506  meshMod.removePoint(pointI, -1);
507  meshChanged = true;
508  }
509  }
510  }
511 
512 
513 
514  // Remove points.
515  forAll(pointRegion_, pointI)
516  {
517  if (pointRemoved(pointI))
518  {
519  meshMod.removePoint(pointI, -1);
520  meshChanged = true;
521  }
522  }
523 
524 
525 
526  const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
527  const faceZoneMesh& faceZones = mesh_.faceZones();
528 
529 
530  // Renumber faces that use points
531  forAll(pointRegion_, pointI)
532  {
533  if (pointRemoved(pointI))
534  {
535  const labelList& changedFaces = pointFaces[pointI];
536 
537  forAll(changedFaces, changedFaceI)
538  {
539  label faceI = changedFaces[changedFaceI];
540 
541  if (!doneFace[faceI])
542  {
543  doneFace[faceI] = true;
544 
545  // Get current zone info
546  label zoneID = faceZones.whichZone(faceI);
547 
548  bool zoneFlip = false;
549 
550  if (zoneID >= 0)
551  {
552  const faceZone& fZone = faceZones[zoneID];
553 
554  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
555  }
556 
557  // Get current connectivity
558  label own = faceOwner[faceI];
559  label nei = -1;
560  label patchID = -1;
561 
562  if (mesh_.isInternalFace(faceI))
563  {
564  nei = faceNeighbour[faceI];
565  }
566  else
567  {
568  patchID = boundaryMesh.whichPatch(faceI);
569  }
570 
571  meshMod.modifyFace
572  (
573  newFaces[faceI], // face
574  faceI, // faceI to change
575  own, // owner
576  nei, // neighbour
577  false, // flipFaceFlux
578  patchID, // patch
579  zoneID,
580  zoneFlip
581  );
582  meshChanged = true;
583  }
584  }
585  }
586  }
587 
588  return meshChanged;
589 }
590 
591 
593 {
594  pointRegion_.setSize(mesh_.nPoints());
595  pointRegion_ = -1;
596  // Reset count, do not remove underlying storage
597  pointRegionMaster_.clear();
598  freeRegions_.clear();
599 }
600 
601 
602 // ************************ vim: set sw=4 sts=4 et: ************************ //