FreeFOAM The Cross-Platform CFD Toolkit
removePoints.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 
27 #include "removePoints.H"
29 #include <OpenFOAM/polyMesh.H>
30 #include "polyTopoChange.H"
34 #include <OpenFOAM/syncTools.H>
35 #include <meshTools/wallPoint.H> // only to use 'greatPoint'
36 #include <meshTools/faceSet.H>
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42 
43 defineTypeNameAndDebug(removePoints, 0);
44 
45 }
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
49 // Change the vertices of the face whilst keeping everything else the same.
50 void Foam::removePoints::modifyFace
51 (
52  const label faceI,
53  const face& newFace,
54  polyTopoChange& meshMod
55 ) const
56 {
57  // Get other face data.
58  label patchI = -1;
59  label owner = mesh_.faceOwner()[faceI];
60  label neighbour = -1;
61 
62  if (mesh_.isInternalFace(faceI))
63  {
64  neighbour = mesh_.faceNeighbour()[faceI];
65  }
66  else
67  {
68  patchI = mesh_.boundaryMesh().whichPatch(faceI);
69  }
70 
71  label zoneID = mesh_.faceZones().whichZone(faceI);
72 
73  bool zoneFlip = false;
74 
75  if (zoneID >= 0)
76  {
77  const faceZone& fZone = mesh_.faceZones()[zoneID];
78 
79  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
80  }
81 
82  meshMod.setAction
83  (
84  polyModifyFace
85  (
86  newFace, // modified face
87  faceI, // label of face being modified
88  owner, // owner
89  neighbour, // neighbour
90  false, // face flip
91  patchI, // patch for face
92  false, // remove from zone
93  zoneID, // zone for face
94  zoneFlip // face flip in zone
95  )
96  );
97 }
98 
99 
100 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
101 
102 // Construct from mesh
103 Foam::removePoints::removePoints
104 (
105  const polyMesh& mesh,
106  const bool undoable
107 )
108 :
109  mesh_(mesh),
110  undoable_(undoable)
111 {}
112 
113 
114 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
115 
117 (
118  const scalar minCos,
119  boolList& pointCanBeDeleted
120 ) const
121 {
122  // Containers to store two edges per point:
123  // -1 : not filled
124  // >= 0 : edge label
125  // -2 : more than two edges for point
126  labelList edge0(mesh_.nPoints(), -1);
127  labelList edge1(mesh_.nPoints(), -1);
128 
129  const edgeList& edges = mesh_.edges();
130 
131  forAll(edges, edgeI)
132  {
133  const edge& e = edges[edgeI];
134 
135  forAll(e, eI)
136  {
137  label pointI = e[eI];
138 
139  if (edge0[pointI] == -2)
140  {
141  // Already too many edges
142  }
143  else if (edge0[pointI] == -1)
144  {
145  // Store first edge using point
146  edge0[pointI] = edgeI;
147  }
148  else
149  {
150  // Already one edge using point. Check second container.
151 
152  if (edge1[pointI] == -1)
153  {
154  // Store second edge using point
155  edge1[pointI] = edgeI;
156  }
157  else
158  {
159  // Third edge using point. Mark.
160  edge0[pointI] = -2;
161  edge1[pointI] = -2;
162  }
163  }
164  }
165  }
166 
167 
168  // Check the ones used by only 2 edges that these are sufficiently in line.
169  const pointField& points = mesh_.points();
170 
171  pointCanBeDeleted.setSize(mesh_.nPoints());
172  pointCanBeDeleted = false;
173  label nDeleted = 0;
174 
175  forAll(edge0, pointI)
176  {
177  if (edge0[pointI] >= 0 && edge1[pointI] >= 0)
178  {
179  // Point used by two edges exactly
180 
181  const edge& e0 = edges[edge0[pointI]];
182  const edge& e1 = edges[edge1[pointI]];
183 
184  label common = e0.commonVertex(e1);
185  label vLeft = e0.otherVertex(common);
186  label vRight = e1.otherVertex(common);
187 
188  vector e0Vec = points[common] - points[vLeft];
189  e0Vec /= mag(e0Vec) + VSMALL;
190 
191  vector e1Vec = points[vRight] - points[common];
192  e1Vec /= mag(e1Vec) + VSMALL;
193 
194  if ((e0Vec & e1Vec) > minCos)
195  {
196  pointCanBeDeleted[pointI] = true;
197  nDeleted++;
198  }
199  }
200  else if (edge0[pointI] == -1)
201  {
202  // point not used at all
203  pointCanBeDeleted[pointI] = true;
204  nDeleted++;
205  }
206  }
207  edge0.clear();
208  edge1.clear();
209 
210 
211  // Protect any points on faces that would collapse down to nothing
212  // No particular intelligence so might protect too many points
213  forAll(mesh_.faces(), faceI)
214  {
215  const face& f = mesh_.faces()[faceI];
216 
217  label nCollapse = 0;
218  forAll(f, fp)
219  {
220  if (pointCanBeDeleted[f[fp]])
221  {
222  nCollapse++;
223  }
224  }
225 
226  if ((f.size() - nCollapse) < 3)
227  {
228  // Just unmark enough points
229  forAll(f, fp)
230  {
231  if (pointCanBeDeleted[f[fp]])
232  {
233  pointCanBeDeleted[f[fp]] = false;
234  --nCollapse;
235  if (nCollapse == 0)
236  {
237  break;
238  }
239  }
240  }
241  }
242  }
243 
244 
245  // Point can be deleted only if all processors want to delete it
247  (
248  mesh_,
249  pointCanBeDeleted,
250  andEqOp<bool>(),
251  true, // null value
252  false // no separation
253  );
254 
255  return returnReduce(nDeleted, sumOp<label>());
256 }
257 
258 
260 (
261  const boolList& pointCanBeDeleted,
262  polyTopoChange& meshMod
263 )
264 {
265  // Count deleted points
266  label nDeleted = 0;
267  forAll(pointCanBeDeleted, pointI)
268  {
269  if (pointCanBeDeleted[pointI])
270  {
271  nDeleted++;
272  }
273  }
274 
275  // Faces (in mesh face labels) affected by points removed. Will hopefully
276  // be only a few.
277  labelHashSet facesAffected(4*nDeleted);
278 
279 
280  // Undo: from global mesh point to index in savedPoint_
281  Map<label> pointToSaved;
282 
283  // Size undo storage
284  if (undoable_)
285  {
286  savedPoints_.setSize(nDeleted);
287  pointToSaved.resize(2*nDeleted);
288  }
289 
290 
291  // Remove points
292  // ~~~~~~~~~~~~~
293 
294  nDeleted = 0;
295 
296  forAll(pointCanBeDeleted, pointI)
297  {
298  if (pointCanBeDeleted[pointI])
299  {
300  if (undoable_)
301  {
302  pointToSaved.insert(pointI, nDeleted);
303  savedPoints_[nDeleted++] = mesh_.points()[pointI];
304  }
305  meshMod.setAction(polyRemovePoint(pointI));
306 
307  // Store faces affected
308  const labelList& pFaces = mesh_.pointFaces()[pointI];
309 
310  forAll(pFaces, i)
311  {
312  facesAffected.insert(pFaces[i]);
313  }
314  }
315  }
316 
317 
318 
319  // Update faces
320  // ~~~~~~~~~~~~
321 
322 
323  if (undoable_)
324  {
325  savedFaceLabels_.setSize(facesAffected.size());
326  savedFaces_.setSize(facesAffected.size());
327  }
328  label nSaved = 0;
329 
330  forAllConstIter(labelHashSet, facesAffected, iter)
331  {
332  label faceI = iter.key();
333 
334  const face& f = mesh_.faces()[faceI];
335 
336  face newFace(f.size());
337 
338  label newI = 0;
339 
340  forAll(f, fp)
341  {
342  label pointI = f[fp];
343 
344  if (!pointCanBeDeleted[pointI])
345  {
346  newFace[newI++] = pointI;
347  }
348  }
349  newFace.setSize(newI);
350 
351  // Actually change the face to the new vertices
352  modifyFace(faceI, newFace, meshMod);
353 
354  // Save the face. Negative indices are into savedPoints_
355  if (undoable_)
356  {
357  savedFaceLabels_[nSaved] = faceI;
358 
359  face& savedFace = savedFaces_[nSaved++];
360  savedFace.setSize(f.size());
361 
362  forAll(f, fp)
363  {
364  label pointI = f[fp];
365 
366  if (pointCanBeDeleted[pointI])
367  {
368  savedFace[fp] = -pointToSaved[pointI]-1;
369  }
370  else
371  {
372  savedFace[fp] = pointI;
373  }
374  }
375  }
376  }
377 
378  if (undoable_)
379  {
380  // DEBUG: Compare the stored faces with the current ones.
381  if (debug)
382  {
383  forAll(savedFaceLabels_, saveI)
384  {
385  // Points from the mesh
386  List<point> meshPoints
387  (
389  (
390  mesh_.points(),
391  mesh_.faces()[savedFaceLabels_[saveI]] // mesh face
392  )
393  );
394 
395  // Points from the saved face
396  List<point> keptPoints
397  (
399  (
400  mesh_.points(),
401  savedPoints_,
402  savedFaces_[saveI] // saved face
403  )
404  );
405 
406  if (meshPoints != keptPoints)
407  {
408  FatalErrorIn("setRefinement")
409  << "faceI:" << savedFaceLabels_[saveI] << nl
410  << "meshPoints:" << meshPoints << nl
411  << "keptPoints:" << keptPoints << nl
412  << abort(FatalError);
413  }
414  }
415  }
416  }
417 }
418 
419 
421 {
422  if (undoable_)
423  {
424  forAll(savedFaceLabels_, localI)
425  {
426  if (savedFaceLabels_[localI] >= 0)
427  {
428  label newFaceI = map.reverseFaceMap()[savedFaceLabels_[localI]];
429 
430  if (newFaceI == -1)
431  {
433  (
434  "removePoints::updateMesh(const mapPolyMesh&)"
435  ) << "Old face " << savedFaceLabels_[localI]
436  << " seems to have dissapeared."
437  << abort(FatalError);
438  }
439  savedFaceLabels_[localI] = newFaceI;
440  }
441  }
442 
443 
444  // Renumber mesh vertices (indices >=0). Leave saved vertices
445  // (<0) intact.
446  forAll(savedFaces_, i)
447  {
448  face& f = savedFaces_[i];
449 
450  forAll(f, fp)
451  {
452  label pointI = f[fp];
453 
454  if (pointI >= 0)
455  {
456  f[fp] = map.reversePointMap()[pointI];
457 
458  if (f[fp] == -1)
459  {
461  (
462  "removePoints::updateMesh(const mapPolyMesh&)"
463  ) << "Old point " << pointI
464  << " seems to have dissapeared."
465  << abort(FatalError);
466  }
467  }
468  }
469  }
470 
471 
472  // DEBUG: Compare the stored faces with the current ones.
473  if (debug)
474  {
475  forAll(savedFaceLabels_, saveI)
476  {
477  if (savedFaceLabels_[saveI] >= 0)
478  {
479  const face& f = mesh_.faces()[savedFaceLabels_[saveI]];
480 
481  // Get kept points of saved faces.
482  const face& savedFace = savedFaces_[saveI];
483 
484  face keptFace(savedFace.size());
485  label keptFp = 0;
486 
487  forAll(savedFace, fp)
488  {
489  label pointI = savedFace[fp];
490 
491  if (pointI >= 0)
492  {
493  keptFace[keptFp++] = pointI;
494  }
495  }
496  keptFace.setSize(keptFp);
497 
498  // Compare as faces (since f might have rotated and
499  // face::operator== takes care of this)
500  if (keptFace != f)
501  {
502  FatalErrorIn("setRefinement")
503  << "faceI:" << savedFaceLabels_[saveI] << nl
504  << "face:" << f << nl
505  << "keptFace:" << keptFace << nl
506  << "saved points:"
508  (
509  mesh_.points(),
510  savedPoints_,
511  savedFace
512  )() << nl
513  << abort(FatalError);
514  }
515  }
516  }
517  }
518  }
519 }
520 
521 
522 // Given list of faces to undo picks up the local indices of the faces
523 // to restore. Additionally it also picks up all the faces that use
524 // any of the deleted points.
526 (
527  const labelList& undoFaces,
528  labelList& localFaces,
529  labelList& localPoints
530 ) const
531 {
532  if (!undoable_)
533  {
535  (
536  "removePoints::getUnrefimentSet(const labelList&"
537  ", labelList&, labelList&) const"
538  ) << "removePoints not constructed with"
539  << " unrefinement capability."
540  << abort(FatalError);
541  }
542 
543  if (debug)
544  {
545  // Check if synced.
546  faceSet undoFacesSet(mesh_, "undoFacesSet", undoFaces);
547  label sz = undoFacesSet.size();
548 
549  undoFacesSet.sync(mesh_);
550  if (sz != undoFacesSet.size())
551  {
553  (
554  "removePoints::getUnrefimentSet(const labelList&"
555  ", labelList&, labelList&) const"
556  ) << "undoFaces not synchronised across coupled faces." << endl
557  << "Size before sync:" << sz
558  << " after sync:" << undoFacesSet.size()
559  << abort(FatalError);
560  }
561  }
562 
563 
564  // Problem: input undoFaces are synced. But e.g.
565  // two faces, A (uncoupled) and B(coupled) share a deleted point. A gets
566  // passed in to be restored. Now picking up the deleted point and the
567  // original faces using it picks up face B. But not its coupled neighbour!
568  // Problem is that we cannot easily synchronise the deleted points
569  // (e.g. syncPointList) since they're not in the mesh anymore - only the
570  // faces are. So we have to collect the points-to-restore as indices
571  // in the faces (which is information we can synchronise)
572 
573 
574 
575  // Mark points-to-restore
576  labelHashSet localPointsSet(undoFaces.size());
577 
578  {
579  // Create set of faces to be restored
580  labelHashSet undoFacesSet(undoFaces);
581 
582  forAll(savedFaceLabels_, saveI)
583  {
584  if (savedFaceLabels_[saveI] < 0)
585  {
587  (
588  "removePoints::getUnrefimentSet(const labelList&"
589  ", labelList&, labelList&) const"
590  ) << "Illegal face label " << savedFaceLabels_[saveI]
591  << " at index " << saveI
592  << abort(FatalError);
593  }
594 
595  if (undoFacesSet.found(savedFaceLabels_[saveI]))
596  {
597  const face& savedFace = savedFaces_[saveI];
598 
599  forAll(savedFace, fp)
600  {
601  if (savedFace[fp] < 0)
602  {
603  label savedPointI = -savedFace[fp]-1;
604 
605  if (savedPoints_[savedPointI] == wallPoint::greatPoint)
606  {
608  (
609  "removePoints::getUnrefimentSet"
610  "(const labelList&, labelList&, labelList&)"
611  " const"
612  ) << "Trying to restore point " << savedPointI
613  << " from mesh face " << savedFaceLabels_[saveI]
614  << " saved face:" << savedFace
615  << " which has already been undone."
616  << abort(FatalError);
617  }
618 
619  localPointsSet.insert(savedPointI);
620  }
621  }
622  }
623  }
624 
625 
626  // Per boundary face, per index in face whether the point needs
627  // restoring. Note that this is over all saved faces, not just over
628  // the ones in undoFaces.
629 
630  boolListList faceVertexRestore(mesh_.nFaces()-mesh_.nInternalFaces());
631 
632  // Populate with my local points-to-restore.
633  forAll(savedFaces_, saveI)
634  {
635  label bFaceI = savedFaceLabels_[saveI] - mesh_.nInternalFaces();
636 
637  if (bFaceI >= 0)
638  {
639  const face& savedFace = savedFaces_[saveI];
640 
641  boolList& fRestore = faceVertexRestore[bFaceI];
642 
643  fRestore.setSize(savedFace.size());
644  fRestore = false;
645 
646  forAll(savedFace, fp)
647  {
648  if (savedFace[fp] < 0)
649  {
650  label savedPointI = -savedFace[fp]-1;
651 
652  if (localPointsSet.found(savedPointI))
653  {
654  fRestore[fp] = true;
655  }
656  }
657  }
658  }
659  }
660 
662  (
663  mesh_,
664  faceVertexRestore,
665  faceEqOp<bool, orEqOp>(), // special operator to handle faces
666  false // no separation
667  );
668 
669  // So now if any of the points-to-restore is used by any coupled face
670  // anywhere the corresponding index in faceVertexRestore will be set.
671 
672  // Now combine the localPointSet and the (sychronised)
673  // boundary-points-to-restore.
674 
675  forAll(savedFaces_, saveI)
676  {
677  label bFaceI = savedFaceLabels_[saveI] - mesh_.nInternalFaces();
678 
679  if (bFaceI >= 0)
680  {
681  const boolList& fRestore = faceVertexRestore[bFaceI];
682 
683  const face& savedFace = savedFaces_[saveI];
684 
685  forAll(fRestore, fp)
686  {
687  // Does neighbour require point restored?
688  if (fRestore[fp])
689  {
690  if (savedFace[fp] >= 0)
691  {
693  (
694  "removePoints::getUnrefimentSet"
695  "(const labelList&, labelList&, labelList&)"
696  " const"
697  ) << "Problem: on coupled face:"
698  << savedFaceLabels_[saveI]
699  << " fc:"
700  << mesh_.faceCentres()[savedFaceLabels_[saveI]]
701  << endl
702  << " my neighbour tries to restore the vertex"
703  << " at index " << fp
704  << " whereas my saved face:" << savedFace
705  << " does not indicate a deleted vertex"
706  << " at that position."
707  << abort(FatalError);
708  }
709 
710  label savedPointI = -savedFace[fp]-1;
711 
712  localPointsSet.insert(savedPointI);
713  }
714  }
715  }
716  }
717  }
718 
719  localPoints = localPointsSet.toc();
720 
721 
722  // Collect all saved faces using any localPointsSet
723  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
724 
725  labelHashSet localFacesSet(2*undoFaces.size());
726 
727  forAll(savedFaces_, saveI)
728  {
729  const face& savedFace = savedFaces_[saveI];
730 
731  forAll(savedFace, fp)
732  {
733  if (savedFace[fp] < 0)
734  {
735  label savedPointI = -savedFace[fp]-1;
736 
737  if (localPointsSet.found(savedPointI))
738  {
739  localFacesSet.insert(saveI);
740  }
741  }
742  }
743  }
744  localFaces = localFacesSet.toc();
745 
746 
747  // Note that at this point the localFaces to restore can contain points
748  // that are not going to be restored! The localFaces though will
749  // be guaranteed to be all the ones affected by the restoration of the
750  // localPoints.
751 }
752 
753 
755 (
756  const labelList& localFaces,
757  const labelList& localPoints,
758  polyTopoChange& meshMod
759 )
760 {
761  if (!undoable_)
762  {
764  (
765  "removePoints::setUnrefinement(const labelList&"
766  ", labelList&, polyTopoChange&)"
767  ) << "removePoints not constructed with"
768  << " unrefinement capability."
769  << abort(FatalError);
770  }
771 
772 
773  // Per savedPoint -1 or the restored point label
774  labelList addedPoints(savedPoints_.size(), -1);
775 
776  forAll(localPoints, i)
777  {
778  label localI = localPoints[i];
779 
780  if (savedPoints_[localI] == wallPoint::greatPoint)
781  {
783  (
784  "removePoints::setUnrefinement(const labelList&"
785  ", labelList&, polyTopoChange&)"
786  ) << "Saved point " << localI << " already restored!"
787  << abort(FatalError);
788  }
789 
790  addedPoints[localI] = meshMod.setAction
791  (
793  (
794  savedPoints_[localI], // point
795  -1, // master point
796  -1, // zone for point
797  true // supports a cell
798  )
799  );
800 
801  // Mark the restored points so they are not restored again.
802  savedPoints_[localI] = wallPoint::greatPoint;
803  }
804 
805  forAll(localFaces, i)
806  {
807  label saveI = localFaces[i];
808 
809  // Modify indices into saved points (so < 0) to point to the
810  // added points.
811  face& savedFace = savedFaces_[saveI];
812 
813  face newFace(savedFace.size());
814  label newFp = 0;
815 
816  bool hasSavedPoints = false;
817 
818  forAll(savedFace, fp)
819  {
820  if (savedFace[fp] < 0)
821  {
822  label addedPointI = addedPoints[-savedFace[fp]-1];
823 
824  if (addedPointI != -1)
825  {
826  savedFace[fp] = addedPointI;
827  newFace[newFp++] = addedPointI;
828  }
829  else
830  {
831  hasSavedPoints = true;
832  }
833  }
834  else
835  {
836  newFace[newFp++] = savedFace[fp];
837  }
838  }
839  newFace.setSize(newFp);
840 
841  modifyFace(savedFaceLabels_[saveI], newFace, meshMod);
842 
843  if (!hasSavedPoints)
844  {
845  // Face fully restored. Mark for compaction later on
846  savedFaceLabels_[saveI] = -1;
847  savedFaces_[saveI].clear();
848  }
849  }
850 
851 
852  // Compact face labels
853  label newSaveI = 0;
854 
855  forAll(savedFaceLabels_, saveI)
856  {
857  if (savedFaceLabels_[saveI] != -1)
858  {
859  if (newSaveI != saveI)
860  {
861  savedFaceLabels_[newSaveI] = savedFaceLabels_[saveI];
862  savedFaces_[newSaveI].transfer(savedFaces_[saveI]);
863  }
864  newSaveI++;
865  }
866  }
867 
868  savedFaceLabels_.setSize(newSaveI);
869  savedFaces_.setSize(newSaveI);
870 
871 
872  // Check that all faces have been restored that use any restored points
873  if (debug)
874  {
875  forAll(savedFaceLabels_, saveI)
876  {
877  const face& savedFace = savedFaces_[saveI];
878 
879  forAll(savedFace, fp)
880  {
881  if (savedFace[fp] < 0)
882  {
883  label addedPointI = addedPoints[-savedFace[fp]-1];
884 
885  if (addedPointI != -1)
886  {
887  FatalErrorIn("setUnrefinement")
888  << "Face:" << savedFaceLabels_[saveI]
889  << " savedVerts:" << savedFace
890  << " uses restored point:" << -savedFace[fp]-1
891  << " with new pointlabel:" << addedPointI
892  << abort(FatalError);
893  }
894  }
895  }
896  }
897  }
898 }
899 
900 
901 // ************************ vim: set sw=4 sts=4 et: ************************ //