FreeFOAM The Cross-Platform CFD Toolkit
collapseEdges.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 Application
25  collapseEdges
26 
27 Description
28  Collapse short edges and combines edges that are in line.
29 
30  - collapse short edges. Length of edges to collapse provided as argument.
31  - merge two edges if they are in line. Maximum angle provided as argument.
32  - remove unused points.
33 
34  Cannot remove cells. Can remove faces and points but does not check
35  for nonsense resulting topology.
36 
37  When collapsing an edge with one point on the boundary it will leave
38  the boundary point intact. When both points inside it chooses random. When
39  both points on boundary random again. Note: it should in fact use features
40  where if one point is on a feature it collapses to that one. Alas we don't
41  have features on a polyMesh.
42 
43 Usage
44 
45  - collapseEdges [OPTIONS] <edge length [m]> <merge angle [degrees]>
46 
47  @param <edge length [m]> \n
48  @todo Detailed description of argument.
49 
50  @param <merge angle [degrees]> \n
51  @todo Detailed description of argument.
52 
53  @param -overwrite \n
54  Overwrite existing data.
55 
56  @param -case <dir>\n
57  Case directory.
58 
59  @param -help \n
60  Display help message.
61 
62  @param -doc \n
63  Display Doxygen API documentation page for this application.
64 
65  @param -srcDoc \n
66  Display Doxygen source documentation page for this application.
67 
68 \*---------------------------------------------------------------------------*/
69 
70 #include <OpenFOAM/argList.H>
71 #include <OpenFOAM/Time.H>
75 #include <OpenFOAM/polyMesh.H>
76 #include <OpenFOAM/mapPolyMesh.H>
79 #include <OpenFOAM/SortableList.H>
80 
81 using namespace Foam;
82 
83 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
84 
85 // Get faceEdges in order of face points, i.e. faceEdges[0] is between
86 // f[0] and f[1]
87 labelList getSortedEdges
88 (
89  const edgeList& edges,
90  const labelList& f,
91  const labelList& edgeLabels
92 )
93 {
94  labelList faceEdges(edgeLabels.size(), -1);
95 
96  // Find starting pos in f for every edgeLabels
97  forAll(edgeLabels, i)
98  {
99  label edgeI = edgeLabels[i];
100 
101  const edge& e = edges[edgeI];
102 
103  label fp = findIndex(f, e[0]);
104  label fp1 = f.fcIndex(fp);
105 
106  if (f[fp1] == e[1])
107  {
108  // EdgeI between fp -> fp1
109  faceEdges[fp] = edgeI;
110  }
111  else
112  {
113  // EdgeI between fp-1 -> fp
114  faceEdges[f.rcIndex(fp)] = edgeI;
115  }
116  }
117 
118  return faceEdges;
119 }
120 
121 
122 // Merges edges which are in straight line. I.e. edge split by point.
123 label mergeEdges
124 (
125  const polyMesh& mesh,
126  const scalar maxCos,
127  edgeCollapser& collapser
128 )
129 {
130  const pointField& points = mesh.points();
131  const edgeList& edges = mesh.edges();
132  const labelListList& pointEdges = mesh.pointEdges();
133  const labelList& region = collapser.pointRegion();
134  const labelList& master = collapser.pointRegionMaster();
135 
136  label nCollapsed = 0;
137 
138  forAll(pointEdges, pointI)
139  {
140  const labelList& pEdges = pointEdges[pointI];
141 
142  if (pEdges.size() == 2)
143  {
144  const edge& leftE = edges[pEdges[0]];
145  const edge& rightE = edges[pEdges[1]];
146 
147  // Get the two vertices on both sides of the point
148  label leftV = leftE.otherVertex(pointI);
149  label rightV = rightE.otherVertex(pointI);
150 
151  // Collapse only if none of the points part of merge network
152  // or all of networks with different masters.
153  label midMaster = -1;
154  if (region[pointI] != -1)
155  {
156  midMaster = master[region[pointI]];
157  }
158 
159  label leftMaster = -2;
160  if (region[leftV] != -1)
161  {
162  leftMaster = master[region[leftV]];
163  }
164 
165  label rightMaster = -3;
166  if (region[rightV] != -1)
167  {
168  rightMaster = master[region[rightV]];
169  }
170 
171  if
172  (
173  midMaster != leftMaster
174  && midMaster != rightMaster
175  && leftMaster != rightMaster
176  )
177  {
178  // Check if the two edge are in line
179  vector leftVec = points[pointI] - points[leftV];
180  leftVec /= mag(leftVec) + VSMALL;
181 
182  vector rightVec = points[rightV] - points[pointI];
183  rightVec /= mag(rightVec) + VSMALL;
184 
185  if ((leftVec & rightVec) > maxCos)
186  {
187  // Collapse one (left) side of the edge. Make left vertex
188  // the master.
189  //if (collapser.unaffectedEdge(pEdges[0]))
190  {
191  collapser.collapseEdge(pEdges[0], leftV);
192  nCollapsed++;
193  }
194  }
195  }
196  }
197  }
198 
199  return nCollapsed;
200 }
201 
202 
203 // Return master point edge needs to be collapsed to (or -1)
204 label edgeMaster(const PackedBoolList& boundaryPoint, const edge& e)
205 {
206  label masterPoint = -1;
207 
208  // Collapse edge to boundary point.
209  if (boundaryPoint.get(e[0]))
210  {
211  if (boundaryPoint.get(e[1]))
212  {
213  // Both points on boundary. Choose one to collapse to.
214  // Note: should look at feature edges/points!
215  masterPoint = e[0];
216  }
217  else
218  {
219  masterPoint = e[0];
220  }
221  }
222  else
223  {
224  if (boundaryPoint.get(e[1]))
225  {
226  masterPoint = e[1];
227  }
228  else
229  {
230  // None on boundary. Choose arbitrary.
231  // Note: should look at geometry?
232  masterPoint = e[0];
233  }
234  }
235  return masterPoint;
236 }
237 
238 
239 label collapseSmallEdges
240 (
241  const polyMesh& mesh,
242  const PackedBoolList& boundaryPoint,
243  const scalar minLen,
244  edgeCollapser& collapser
245 )
246 {
247  const pointField& points = mesh.points();
248  const edgeList& edges = mesh.edges();
249 
250  // Collapse all edges that are too small. Choose intelligently which
251  // point to collapse edge to.
252 
253  label nCollapsed = 0;
254 
255  forAll(edges, edgeI)
256  {
257  const edge& e = edges[edgeI];
258 
259  if (e.mag(points) < minLen)
260  {
261  label master = edgeMaster(boundaryPoint, e);
262 
263  if (master != -1) // && collapser.unaffectedEdge(edgeI))
264  {
265  collapser.collapseEdge(edgeI, master);
266  nCollapsed++;
267  }
268  }
269  }
270  return nCollapsed;
271 }
272 
273 
274 // Faces which have edges just larger than collapse length but faces which
275 // are very small. This one tries to collapse them if it can be done with
276 // edge collapse. For faces where a face gets replace by two edges use
277 // collapseFaces
278 label collapseHighAspectFaces
279 (
280  const polyMesh& mesh,
281  const PackedBoolList& boundaryPoint,
282  const scalar areaFac,
283  const scalar edgeRatio,
284  edgeCollapser& collapser
285 )
286 {
287  const pointField& points = mesh.points();
288  const edgeList& edges = mesh.edges();
289  const faceList& faces = mesh.faces();
290  const labelListList& faceEdges = mesh.faceEdges();
291 
292  scalarField magArea(mag(mesh.faceAreas()));
293 
294  label maxIndex = findMax(magArea);
295 
296  scalar minArea = areaFac * magArea[maxIndex];
297 
298  Info<< "Max face area:" << magArea[maxIndex] << endl
299  << "Collapse area factor:" << areaFac << endl
300  << "Collapse area:" << minArea << endl;
301 
302  label nCollapsed = 0;
303 
304  forAll(faces, faceI)
305  {
306  if (magArea[faceI] < minArea)
307  {
308  const face& f = faces[faceI];
309 
310  // Get the edges in face point order
311  labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
312 
313  SortableList<scalar> lengths(fEdges.size());
314  forAll(fEdges, i)
315  {
316  lengths[i] = edges[fEdges[i]].mag(points);
317  }
318  lengths.sort();
319 
320 
321  label edgeI = -1;
322 
323  if (f.size() == 4)
324  {
325  // Compare second largest to smallest
326  if (lengths[2] > edgeRatio*lengths[0])
327  {
328  // Collapse smallest only. Triangle should be cleared
329  // next time around.
330  edgeI = fEdges[lengths.indices()[0]];
331  }
332  }
333  else if (f.size() == 3)
334  {
335  // Compare second largest to smallest
336  if (lengths[1] > edgeRatio*lengths[0])
337  {
338  edgeI = fEdges[lengths.indices()[0]];
339  }
340  }
341 
342 
343  if (edgeI != -1)
344  {
345  label master = edgeMaster(boundaryPoint, edges[edgeI]);
346 
347  if (master != -1)// && collapser.unaffectedEdge(edgeI))
348  {
349  collapser.collapseEdge(edgeI, master);
350  nCollapsed++;
351  }
352  }
353  }
354  }
355 
356  return nCollapsed;
357 }
358 
359 
360 void set(const labelList& elems, const bool val, boolList& status)
361 {
362  forAll(elems, i)
363  {
364  status[elems[i]] = val;
365  }
366 }
367 
368 
369 // Tries to simplify polygons to face of minSize (4=quad, 3=triangle)
370 label simplifyFaces
371 (
372  const polyMesh& mesh,
373  const PackedBoolList& boundaryPoint,
374  const label minSize,
375  const scalar lenGap,
376  edgeCollapser& collapser
377 )
378 {
379  const pointField& points = mesh.points();
380  const edgeList& edges = mesh.edges();
381  const faceList& faces = mesh.faces();
382  const cellList& cells = mesh.cells();
383  const labelListList& faceEdges = mesh.faceEdges();
384  const labelList& faceOwner = mesh.faceOwner();
385  const labelList& faceNeighbour = mesh.faceNeighbour();
386  const labelListList& pointCells = mesh.pointCells();
387  const labelListList& cellEdges = mesh.cellEdges();
388 
389  label nCollapsed = 0;
390 
391  boolList protectedEdge(mesh.nEdges(), false);
392 
393  forAll(faces, faceI)
394  {
395  const face& f = faces[faceI];
396 
397  if
398  (
399  f.size() > minSize
400  && cells[faceOwner[faceI]].size() >= 6
401  && (
402  mesh.isInternalFace(faceI)
403  && cells[faceNeighbour[faceI]].size() >= 6
404  )
405  )
406  {
407  // Get the edges in face point order
408  labelList fEdges(getSortedEdges(edges, f, faceEdges[faceI]));
409 
410  SortableList<scalar> lengths(fEdges.size());
411  forAll(fEdges, i)
412  {
413  lengths[i] = edges[fEdges[i]].mag(points);
414  }
415  lengths.sort();
416 
417 
418  // Now find a gap in length between consecutive elements greater
419  // than lenGap.
420 
421  label gapPos = -1;
422 
423  for (label i = f.size()-1-minSize; i >= 0; --i)
424  {
425  if (lengths[i+1] > lenGap*lengths[i])
426  {
427  gapPos = i;
428 
429  break;
430  }
431  }
432 
433  if (gapPos != -1)
434  {
435  //for (label i = gapPos; i >= 0; --i)
436  label i = 0; // Hack: collapse smallest edge only.
437  {
438  label edgeI = fEdges[lengths.indices()[i]];
439 
440  if (!protectedEdge[edgeI])
441  {
442  const edge& e = edges[edgeI];
443 
444  label master = edgeMaster(boundaryPoint, e);
445 
446  if (master != -1)
447  {
448  collapser.collapseEdge(edgeI, master);
449 
450  // Protect all other edges on all cells using edge
451  // points.
452 
453  const labelList& pCells0 = pointCells[e[0]];
454 
455  forAll(pCells0, i)
456  {
457  set(cellEdges[pCells0[i]], true, protectedEdge);
458  }
459  const labelList& pCells1 = pointCells[e[1]];
460 
461  forAll(pCells1, i)
462  {
463  set(cellEdges[pCells1[i]], true, protectedEdge);
464  }
465 
466  nCollapsed++;
467  }
468  }
469  }
470  }
471  }
472  }
473 
474  return nCollapsed;
475 }
476 
477 
478 // Main program:
479 
480 int main(int argc, char *argv[])
481 {
483  argList::validOptions.insert("overwrite", "");
484  argList::validArgs.append("edge length [m]");
485  argList::validArgs.append("merge angle (degrees)");
486 
487 # include <OpenFOAM/setRootCase.H>
488 # include <OpenFOAM/createTime.H>
489  runTime.functionObjects().off();
490 # include <OpenFOAM/createPolyMesh.H>
491  const word oldInstance = mesh.pointsInstance();
492 
493  scalar minLen(readScalar(IStringStream(args.additionalArgs()[0])()));
494  scalar angle(readScalar(IStringStream(args.additionalArgs()[1])()));
495  bool overwrite = args.optionFound("overwrite");
496 
497  scalar maxCos = Foam::cos(angle*mathematicalConstant::pi/180.0);
498 
499  Info<< "Merging:" << nl
500  << " edges with length less than " << minLen << " meters" << nl
501  << " edges split by a point with edges in line to within " << angle
502  << " degrees" << nl
503  << endl;
504 
505 
506  bool meshChanged = false;
507 
508  while (true)
509  {
510  const faceList& faces = mesh.faces();
511 
512  // Get all points on the boundary
513  PackedBoolList boundaryPoint(mesh.nPoints());
514 
515  label nIntFaces = mesh.nInternalFaces();
516  for (label faceI = nIntFaces; faceI < mesh.nFaces(); faceI++)
517  {
518  const face& f = faces[faceI];
519 
520  forAll(f, fp)
521  {
522  boundaryPoint.set(f[fp], 1);
523  }
524  }
525 
526  // Edge collapsing engine
527  edgeCollapser collapser(mesh);
528 
529 
530  // Collapse all edges that are too small.
531  label nCollapsed =
532  collapseSmallEdges
533  (
534  mesh,
535  boundaryPoint,
536  minLen,
537  collapser
538  );
539  Info<< "Collapsing " << nCollapsed << " small edges" << endl;
540 
541 
542  // Remove midpoints on straight edges.
543  if (nCollapsed == 0)
544  {
545  nCollapsed = mergeEdges(mesh, maxCos, collapser);
546  Info<< "Collapsing " << nCollapsed << " in line edges" << endl;
547  }
548 
549 
550  // Remove small sliver faces that can be collapsed to single edge
551  if (nCollapsed == 0)
552  {
553  nCollapsed =
554  collapseHighAspectFaces
555  (
556  mesh,
557  boundaryPoint,
558  1E-9, // factor of largest face area
559  5, // factor between smallest and largest edge on
560  // face
561  collapser
562  );
563  Info<< "Collapsing " << nCollapsed
564  << " small high aspect ratio faces" << endl;
565  }
566 
567  // Simplify faces to quads wherever possible
568  //if (nCollapsed == 0)
569  //{
570  // nCollapsed =
571  // simplifyFaces
572  // (
573  // mesh,
574  // boundaryPoint,
575  // 4, // minimum size of face
576  // 0.2, // gap in edge lengths on face
577  // collapser
578  // );
579  // Info<< "Collapsing " << nCollapsed << " polygonal faces" << endl;
580  //}
581 
582 
583  if (nCollapsed == 0)
584  {
585  break;
586  }
587 
588  polyTopoChange meshMod(mesh);
589 
590  // Insert mesh refinement into polyTopoChange.
591  collapser.setRefinement(meshMod);
592 
593  // Do all changes
594  Info<< "Morphing ..." << endl;
595 
596  autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
597 
598  collapser.updateMesh(morphMap());
599 
600  if (morphMap().hasMotionPoints())
601  {
602  mesh.movePoints(morphMap().preMotionPoints());
603  }
604 
605  meshChanged = true;
606  }
607 
608  if (meshChanged)
609  {
610  // Write resulting mesh
611  if (!overwrite)
612  {
613  runTime++;
614  }
615  else
616  {
617  mesh.setInstance(oldInstance);
618  }
619 
620  Info<< "Writing collapsed mesh to time " << runTime.timeName() << endl;
621 
622  mesh.write();
623  }
624 
625  Info << "End\n" << endl;
626 
627  return 0;
628 }
629 
630 
631 // ************************ vim: set sw=4 sts=4 et: ************************ //