polyMeshFilter.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) 2012-2016 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 "polyMeshFilter.H"
27 #include "polyMesh.H"
28 #include "fvMesh.H"
29 #include "unitConversion.H"
30 #include "edgeCollapser.H"
31 #include "syncTools.H"
32 #include "polyTopoChange.H"
33 #include "globalIndex.H"
34 #include "PackedBoolList.H"
35 #include "pointSet.H"
36 #include "faceSet.H"
37 #include "cellSet.H"
38 
39 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
40 
41 namespace Foam
42 {
43 defineTypeNameAndDebug(polyMeshFilter, 0);
44 }
45 
46 
47 void Foam::polyMeshFilter::updateSets(const mapPolyMesh& map)
48 {
49  updateSets<pointSet>(map);
50  updateSets<faceSet>(map);
51  updateSets<cellSet>(map);
52 }
53 
54 
55 void Foam::polyMeshFilter::copySets
56 (
57  const polyMesh& oldMesh,
58  const polyMesh& newMesh
59 )
60 {
61  copySets<pointSet>(oldMesh, newMesh);
62  copySets<faceSet>(oldMesh, newMesh);
63  copySets<cellSet>(oldMesh, newMesh);
64 }
65 
66 
68 {
69  polyTopoChange originalMeshToNewMesh(mesh);
70 
71  autoPtr<fvMesh> meshCopy;
72  autoPtr<mapPolyMesh> mapPtr = originalMeshToNewMesh.makeMesh
73  (
74  meshCopy,
75  IOobject
76  (
77  mesh.name(),
78  mesh.polyMesh::instance(),
79  mesh.time(),
82  false
83  ),
84  mesh,
85  true // parallel sync
86  );
87 
88  const mapPolyMesh& map = mapPtr();
89 
90  // Update fields
91  meshCopy().updateMesh(map);
92  if (map.hasMotionPoints())
93  {
94  meshCopy().movePoints(map.preMotionPoints());
95  }
96 
97  copySets(mesh, meshCopy());
98 
99  return meshCopy;
100 }
101 
102 
103 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
104 
105 Foam::label Foam::polyMeshFilter::filterFacesLoop(const label nOriginalBadFaces)
106 {
107  label nBadFaces = labelMax;
108  label nOuterIterations = 0;
109 
110  // Maintain the number of times a point has been part of a bad face
111  labelList pointErrorCount(mesh_.nPoints(), 0);
112 
113  PackedBoolList newErrorPoint(mesh_.nPoints());
115  (
116  mesh_,
117  meshQualityCoeffDict(),
118  newErrorPoint
119  );
120 
121  bool newBadFaces = true;
122 
123  // Main loop
124  // ~~~~~~~~~
125  // It tries and do some collapses, checks the resulting mesh and
126  // 'freezes' some edges (by marking in minEdgeLen) and tries again.
127  // This will iterate ultimately to the situation where every edge is
128  // frozen and nothing gets collapsed.
129  while
130  (
131  nOuterIterations < maxIterations()
132  //&& nBadFaces > nOriginalBadFaces
133  && newBadFaces
134  )
135  {
136  Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
137  << endl;
138 
139  printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
140  printScalarFieldStats("Face Filter Factor", faceFilterFactor_);
141 
142  // Reset the new mesh to the old mesh
143  newMeshPtr_ = copyMesh(mesh_);
144  fvMesh& newMesh = newMeshPtr_();
145 
146  scalarField newMeshFaceFilterFactor = faceFilterFactor_;
147  pointPriority_.reset(new labelList(originalPointPriority_));
148 
149  labelList origToCurrentPointMap(identity(newMesh.nPoints()));
150 
151  {
152  label nInnerIterations = 0;
153  label nPrevLocalCollapse = labelMax;
154 
155  Info<< incrIndent;
156 
157  while (true)
158  {
159  Info<< nl << indent << "Inner iteration = "
160  << nInnerIterations++ << nl << incrIndent << endl;
161 
162  label nLocalCollapse = filterFaces
163  (
164  newMesh,
165  newMeshFaceFilterFactor,
166  origToCurrentPointMap
167  );
168  Info<< decrIndent;
169 
170  if
171  (
172  nLocalCollapse >= nPrevLocalCollapse
173  || nLocalCollapse == 0
174  )
175  {
176  Info<< decrIndent;
177  break;
178  }
179  else
180  {
181  nPrevLocalCollapse = nLocalCollapse;
182  }
183  }
184  }
185 
186  scalarField newMeshMinEdgeLen = minEdgeLen_;
187 
188  {
189  label nInnerIterations = 0;
190  label nPrevLocalCollapse = labelMax;
191 
192  Info<< incrIndent;
193 
194  while (true)
195  {
196  Info<< nl << indent << "Inner iteration = "
197  << nInnerIterations++ << nl << incrIndent << endl;
198 
199  label nLocalCollapse = filterEdges
200  (
201  newMesh,
202  newMeshMinEdgeLen,
203  origToCurrentPointMap
204  );
205  Info<< decrIndent;
206 
207  if
208  (
209  nLocalCollapse >= nPrevLocalCollapse
210  || nLocalCollapse == 0
211  )
212  {
213  Info<< decrIndent;
214  break;
215  }
216  else
217  {
218  nPrevLocalCollapse = nLocalCollapse;
219  }
220  }
221  }
222 
223 
224  // Mesh check
225  // ~~~~~~~~~~~~~~~~~~
226  // Do not allow collapses in regions of error.
227  // Updates minEdgeLen, nRelaxedEdges
228 
229  if (controlMeshQuality())
230  {
231  PackedBoolList isErrorPoint(newMesh.nPoints());
233  (
234  newMesh,
235  meshQualityCoeffDict(),
236  isErrorPoint
237  );
238 
239  Info<< nl << " Number of bad faces : " << nBadFaces << nl
240  << " Number of marked points : "
241  << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
242  << endl;
243 
244  updatePointErrorCount
245  (
246  isErrorPoint,
247  origToCurrentPointMap,
248  pointErrorCount
249  );
250 
251  checkMeshEdgesAndRelaxEdges
252  (
253  newMesh,
254  origToCurrentPointMap,
255  isErrorPoint,
256  pointErrorCount
257  );
258 
259  checkMeshFacesAndRelaxEdges
260  (
261  newMesh,
262  origToCurrentPointMap,
263  isErrorPoint,
264  pointErrorCount
265  );
266 
267  newBadFaces = false;
268  forAll(mesh_.points(), pI)
269  {
270  if
271  (
272  origToCurrentPointMap[pI] >= 0
273  && isErrorPoint[origToCurrentPointMap[pI]]
274  )
275  {
276  if (!newErrorPoint[pI])
277  {
278  newBadFaces = true;
279  break;
280  }
281  }
282  }
283 
284  reduce(newBadFaces, orOp<bool>());
285  }
286  else
287  {
288  return -1;
289  }
290  }
291 
292  return nBadFaces;
293 }
294 
295 
296 Foam::label Foam::polyMeshFilter::filterFaces
297 (
298  polyMesh& newMesh,
299  scalarField& newMeshFaceFilterFactor,
300  labelList& origToCurrentPointMap
301 )
302 {
303  // Per edge collapse status
305 
306  Map<point> collapsePointToLocation(newMesh.nPoints());
307 
308  edgeCollapser collapser(newMesh, collapseFacesCoeffDict());
309 
310  {
311  // Collapse faces
312  labelPair nCollapsedPtEdge = collapser.markSmallSliverFaces
313  (
314  newMeshFaceFilterFactor,
315  pointPriority_(),
316  collapseEdge,
317  collapsePointToLocation
318  );
319 
320  label nCollapsed = 0;
321  forAll(nCollapsedPtEdge, collapseTypeI)
322  {
323  nCollapsed += nCollapsedPtEdge[collapseTypeI];
324  }
325 
326  reduce(nCollapsed, sumOp<label>());
327 
328  label nToPoint = returnReduce(nCollapsedPtEdge.first(), sumOp<label>());
329  label nToEdge = returnReduce(nCollapsedPtEdge.second(), sumOp<label>());
330 
331  Info<< indent
332  << "Collapsing " << nCollapsed << " faces "
333  << "(to point = " << nToPoint << ", to edge = " << nToEdge << ")"
334  << endl;
335 
336  if (nCollapsed == 0)
337  {
338  return 0;
339  }
340  }
341 
342  // Merge edge collapses into consistent collapse-network.
343  // Make sure no cells get collapsed.
344  List<pointEdgeCollapse> allPointInfo;
345  const globalIndex globalPoints(newMesh.nPoints());
346 
347  collapser.consistentCollapse
348  (
349  globalPoints,
350  pointPriority_(),
351  collapsePointToLocation,
352  collapseEdge,
353  allPointInfo
354  );
355 
356  label nLocalCollapse = collapseEdge.count();
357 
358  reduce(nLocalCollapse, sumOp<label>());
359  Info<< nl << indent << "Collapsing " << nLocalCollapse
360  << " edges after synchronisation and PointEdgeWave" << endl;
361 
362  if (nLocalCollapse == 0)
363  {
364  return 0;
365  }
366 
367  {
368  // Apply collapses to current mesh
369  polyTopoChange newMeshMod(newMesh);
370 
371  // Insert mesh refinement into polyTopoChange.
372  collapser.setRefinement(allPointInfo, newMeshMod);
373 
374  Info<< indent << "Apply changes to the current mesh" << endl;
375 
376  // Apply changes to current mesh
377  autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
378  (
379  newMesh,
380  false
381  );
382  const mapPolyMesh& newMap = newMapPtr();
383 
384  // Update fields
385  newMesh.updateMesh(newMap);
386  if (newMap.hasMotionPoints())
387  {
388  newMesh.movePoints(newMap.preMotionPoints());
389  }
390  updateSets(newMap);
391 
392  updatePointPriorities(newMesh, newMap.pointMap());
393 
394  mapOldMeshFaceFieldToNewMesh
395  (
396  newMesh,
397  newMap.faceMap(),
398  newMeshFaceFilterFactor
399  );
400 
401  updateOldToNewPointMap
402  (
403  newMap.reversePointMap(),
404  origToCurrentPointMap
405  );
406  }
407 
408  return nLocalCollapse;
409 }
410 
411 
412 Foam::label Foam::polyMeshFilter::filterEdges
413 (
414  polyMesh& newMesh,
415  scalarField& newMeshMinEdgeLen,
416  labelList& origToCurrentPointMap
417 )
418 {
419  // Per edge collapse status
421 
422  Map<point> collapsePointToLocation(newMesh.nPoints());
423 
424  edgeCollapser collapser(newMesh, collapseFacesCoeffDict());
425 
426  // Work out which edges to collapse
427  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
428  // This is by looking at minEdgeLen (to avoid frozen edges)
429  // and marking in collapseEdge.
430  label nSmallCollapsed = collapser.markSmallEdges
431  (
432  newMeshMinEdgeLen,
433  pointPriority_(),
434  collapseEdge,
435  collapsePointToLocation
436  );
437 
438  reduce(nSmallCollapsed, sumOp<label>());
439  Info<< indent << "Collapsing " << nSmallCollapsed
440  << " small edges" << endl;
441 
442  // Merge inline edges
443  label nMerged = collapser.markMergeEdges
444  (
445  maxCos(),
446  pointPriority_(),
447  collapseEdge,
448  collapsePointToLocation
449  );
450 
451  reduce(nMerged, sumOp<label>());
452  Info<< indent << "Collapsing " << nMerged << " in line edges"
453  << endl;
454 
455  if (nMerged + nSmallCollapsed == 0)
456  {
457  return 0;
458  }
459 
460  // Merge edge collapses into consistent collapse-network.
461  // Make sure no cells get collapsed.
462  List<pointEdgeCollapse> allPointInfo;
463  const globalIndex globalPoints(newMesh.nPoints());
464 
465  collapser.consistentCollapse
466  (
467  globalPoints,
468  pointPriority_(),
469  collapsePointToLocation,
470  collapseEdge,
471  allPointInfo
472  );
473 
474  label nLocalCollapse = collapseEdge.count();
475 
476  reduce(nLocalCollapse, sumOp<label>());
477  Info<< nl << indent << "Collapsing " << nLocalCollapse
478  << " edges after synchronisation and PointEdgeWave" << endl;
479 
480  if (nLocalCollapse == 0)
481  {
482  return 0;
483  }
484 
485  // Apply collapses to current mesh
486  polyTopoChange newMeshMod(newMesh);
487 
488  // Insert mesh refinement into polyTopoChange.
489  collapser.setRefinement(allPointInfo, newMeshMod);
490 
491  Info<< indent << "Apply changes to the current mesh" << endl;
492 
493  // Apply changes to current mesh
494  autoPtr<mapPolyMesh> newMapPtr = newMeshMod.changeMesh
495  (
496  newMesh,
497  false
498  );
499  const mapPolyMesh& newMap = newMapPtr();
500 
501  // Update fields
502  newMesh.updateMesh(newMap);
503  if (newMap.hasMotionPoints())
504  {
505  newMesh.movePoints(newMap.preMotionPoints());
506  }
507  updateSets(newMap);
508 
509  // Synchronise the factors
510  mapOldMeshEdgeFieldToNewMesh
511  (
512  newMesh,
513  newMap.pointMap(),
514  newMeshMinEdgeLen
515  );
516 
517  updateOldToNewPointMap
518  (
519  newMap.reversePointMap(),
520  origToCurrentPointMap
521  );
522 
523  updatePointPriorities(newMesh, newMap.pointMap());
524 
525  return nLocalCollapse;
526 }
527 
528 
529 void Foam::polyMeshFilter::updatePointErrorCount
530 (
531  const PackedBoolList& isErrorPoint,
532  const labelList& oldToNewMesh,
533  labelList& pointErrorCount
534 ) const
535 {
536  forAll(mesh_.points(), pI)
537  {
538  if (isErrorPoint[oldToNewMesh[pI]])
539  {
540  pointErrorCount[pI]++;
541  }
542  }
543 }
544 
545 
546 void Foam::polyMeshFilter::checkMeshEdgesAndRelaxEdges
547 (
548  const polyMesh& newMesh,
549  const labelList& oldToNewMesh,
550  const PackedBoolList& isErrorPoint,
551  const labelList& pointErrorCount
552 )
553 {
554  const edgeList& edges = mesh_.edges();
555 
556  forAll(edges, edgeI)
557  {
558  const edge& e = edges[edgeI];
559  label newStart = oldToNewMesh[e[0]];
560  label newEnd = oldToNewMesh[e[1]];
561 
562  if
563  (
564  pointErrorCount[e[0]] >= maxPointErrorCount()
565  || pointErrorCount[e[1]] >= maxPointErrorCount()
566  )
567  {
568  minEdgeLen_[edgeI] = -1;
569  }
570 
571  if
572  (
573  (newStart >= 0 && isErrorPoint[newStart])
574  || (newEnd >= 0 && isErrorPoint[newEnd])
575  )
576  {
577  minEdgeLen_[edgeI] *= edgeReductionFactor();
578  }
579  }
580 
581  syncTools::syncEdgeList(mesh_, minEdgeLen_, minEqOp<scalar>(), scalar(0.0));
582 
583  for (label smoothIter = 0; smoothIter < maxSmoothIters(); ++smoothIter)
584  {
585  // Smooth minEdgeLen
586  forAll(mesh_.edges(), edgeI)
587  {
588  const edge& e = mesh_.edges()[edgeI];
589 
590  scalar sumMinEdgeLen = 0;
591  label nEdges = 0;
592 
593  forAll(e, pointi)
594  {
595  const labelList& pEdges = mesh_.pointEdges()[e[pointi]];
596 
597  forAll(pEdges, pEdgeI)
598  {
599  const label pEdge = pEdges[pEdgeI];
600  sumMinEdgeLen += minEdgeLen_[pEdge];
601  nEdges++;
602  }
603  }
604 
605  minEdgeLen_[edgeI] = min
606  (
607  minEdgeLen_[edgeI],
608  sumMinEdgeLen/nEdges
609  );
610  }
611 
613  (
614  mesh_,
615  minEdgeLen_,
616  minEqOp<scalar>(),
617  scalar(0.0)
618  );
619  }
620 }
621 
622 
623 void Foam::polyMeshFilter::checkMeshFacesAndRelaxEdges
624 (
625  const polyMesh& newMesh,
626  const labelList& oldToNewMesh,
627  const PackedBoolList& isErrorPoint,
628  const labelList& pointErrorCount
629 )
630 {
631  const faceList& faces = mesh_.faces();
632 
633  forAll(faces, facei)
634  {
635  const face& f = faces[facei];
636 
637  forAll(f, fpI)
638  {
639  const label ptIndex = oldToNewMesh[f[fpI]];
640 
641  if (pointErrorCount[f[fpI]] >= maxPointErrorCount())
642  {
643  faceFilterFactor_[facei] = -1;
644  }
645 
646  if (isErrorPoint[ptIndex])
647  {
648  faceFilterFactor_[facei] *= faceReductionFactor();
649 
650  break;
651  }
652  }
653  }
654 
655  syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());
656 
657  for (label smoothIter = 0; smoothIter < maxSmoothIters(); ++smoothIter)
658  {
659  // Smooth faceFilterFactor
660  forAll(faces, facei)
661  {
662  const labelList& fEdges = mesh_.faceEdges()[facei];
663 
664  scalar sumFaceFilterFactors = 0;
665  label nFaces = 0;
666 
667  // This is important: Only smooth around faces that share an
668  // edge with a bad face
669  bool skipFace = true;
670 
671  forAll(fEdges, fEdgeI)
672  {
673  const labelList& eFaces = mesh_.edgeFaces()[fEdges[fEdgeI]];
674 
675  forAll(eFaces, eFacei)
676  {
677  const label eFace = eFaces[eFacei];
678 
679  const face& f = faces[eFace];
680 
681  forAll(f, fpI)
682  {
683  const label ptIndex = oldToNewMesh[f[fpI]];
684 
685  if (isErrorPoint[ptIndex])
686  {
687  skipFace = false;
688  break;
689  }
690  }
691 
692  if (eFace != facei)
693  {
694  sumFaceFilterFactors += faceFilterFactor_[eFace];
695  nFaces++;
696  }
697  }
698  }
699 
700  if (skipFace)
701  {
702  continue;
703  }
704 
705  faceFilterFactor_[facei] = min
706  (
707  faceFilterFactor_[facei],
708  sumFaceFilterFactors/nFaces
709  );
710  }
711 
712  // Face filter factor needs to be synchronised!
713  syncTools::syncFaceList(mesh_, faceFilterFactor_, minEqOp<scalar>());
714  }
715 }
716 
717 
718 void Foam::polyMeshFilter::updatePointPriorities
719 (
720  const polyMesh& newMesh,
721  const labelList& pointMap
722 )
723 {
724  labelList newPointPriority(newMesh.nPoints(), labelMin);
725  const labelList& currPointPriority = pointPriority_();
726 
727  forAll(newPointPriority, ptI)
728  {
729  const label newPointToOldPoint = pointMap[ptI];
730  const label origPointPriority = currPointPriority[newPointToOldPoint];
731 
732  newPointPriority[ptI] = max(origPointPriority, newPointPriority[ptI]);
733  }
734 
736  (
737  newMesh,
738  newPointPriority,
739  maxEqOp<label>(),
740  labelMin
741  );
742 
743  pointPriority_.reset(new labelList(newPointPriority));
744 }
745 
746 
747 void Foam::polyMeshFilter::printScalarFieldStats
748 (
749  const string desc,
750  const scalarField& fld
751 ) const
752 {
753  scalar sum = 0;
754  scalar validElements = 0;
755  scalar min = GREAT;
756  scalar max = -GREAT;
757 
758  forAll(fld, i)
759  {
760  const scalar fldElement = fld[i];
761 
762  if (fldElement >= 0)
763  {
764  sum += fldElement;
765 
766  if (fldElement < min)
767  {
768  min = fldElement;
769  }
770 
771  if (fldElement > max)
772  {
773  max = fldElement;
774  }
775 
776  validElements++;
777  }
778  }
779 
780  reduce(sum, sumOp<scalar>());
781  reduce(min, minOp<scalar>());
782  reduce(max, maxOp<scalar>());
783  reduce(validElements, sumOp<label>());
784  const label totFieldSize = returnReduce(fld.size(), sumOp<label>());
785 
786  Info<< incrIndent << indent << desc
787  << ": min = " << min
788  << " av = " << sum/(validElements + SMALL)
789  << " max = " << max << nl
790  << indent
791  << " " << validElements << " / " << totFieldSize << " elements used"
792  << decrIndent << endl;
793 }
794 
795 
796 void Foam::polyMeshFilter::mapOldMeshEdgeFieldToNewMesh
797 (
798  const polyMesh& newMesh,
799  const labelList& pointMap,
800  scalarField& newMeshMinEdgeLen
801 ) const
802 {
803  scalarField tmp(newMesh.nEdges());
804 
805  const edgeList& newEdges = newMesh.edges();
806 
807  forAll(newEdges, newEdgeI)
808  {
809  const edge& newEdge = newEdges[newEdgeI];
810  const label pStart = newEdge.start();
811  const label pEnd = newEdge.end();
812 
813  tmp[newEdgeI] = min
814  (
815  newMeshMinEdgeLen[pointMap[pStart]],
816  newMeshMinEdgeLen[pointMap[pEnd]]
817  );
818  }
819 
820  newMeshMinEdgeLen.transfer(tmp);
821 
823  (
824  newMesh,
825  newMeshMinEdgeLen,
826  maxEqOp<scalar>(),
827  scalar(0.0)
828  );
829 }
830 
831 
832 void Foam::polyMeshFilter::mapOldMeshFaceFieldToNewMesh
833 (
834  const polyMesh& newMesh,
835  const labelList& faceMap,
836  scalarField& newMeshFaceFilterFactor
837 ) const
838 {
839  scalarField tmp(newMesh.nFaces());
840 
841  forAll(faceMap, newFacei)
842  {
843  const label oldFacei = faceMap[newFacei];
844 
845  tmp[newFacei] = newMeshFaceFilterFactor[oldFacei];
846  }
847 
848  newMeshFaceFilterFactor.transfer(tmp);
849 
851  (
852  newMesh,
853  newMeshFaceFilterFactor,
855  );
856 }
857 
858 
859 void Foam::polyMeshFilter::updateOldToNewPointMap
860 (
861  const labelList& currToNew,
862  labelList& origToCurrentPointMap
863 ) const
864 {
865  forAll(origToCurrentPointMap, origPointi)
866  {
867  label oldPointi = origToCurrentPointMap[origPointi];
868 
869  if (oldPointi != -1)
870  {
871  label newPointi = currToNew[oldPointi];
872 
873  if (newPointi >= 0)
874  {
875  origToCurrentPointMap[origPointi] = newPointi;
876  }
877  else if (newPointi == -1)
878  {
879  origToCurrentPointMap[origPointi] = -1;
880  }
881  else
882  {
883  origToCurrentPointMap[origPointi] = -newPointi-2;
884  }
885  }
886  }
887 }
888 
889 
890 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
891 
892 Foam::polyMeshFilter::polyMeshFilter(const fvMesh& mesh)
893 :
895  (
897  (
898  IOobject
899  (
900  "collapseDict",
901  mesh.time().system(),
902  mesh.time(),
903  IOobject::MUST_READ,
904  IOobject::NO_WRITE
905  )
906  )
907  ),
908  mesh_(mesh),
909  newMeshPtr_(),
910  originalPointPriority_(mesh.nPoints(), labelMin),
911  pointPriority_(),
912  minEdgeLen_(),
913  faceFilterFactor_()
914 {
916 }
917 
918 
919 Foam::polyMeshFilter::polyMeshFilter
920 (
921  const fvMesh& mesh,
922  const labelList& pointPriority
923 )
924 :
926  (
928  (
929  IOobject
930  (
931  "collapseDict",
932  mesh.time().system(),
933  mesh.time(),
936  )
937  )
938  ),
939  mesh_(mesh),
940  newMeshPtr_(),
941  originalPointPriority_(pointPriority),
942  pointPriority_(),
943  minEdgeLen_(),
944  faceFilterFactor_()
945 {
947 }
948 
949 Foam::polyMeshFilter::polyMeshFilter
950 (
951  const fvMesh& mesh,
952  const labelList& pointPriority,
953  const dictionary& dict
954 )
955 :
957  mesh_(mesh),
958  newMeshPtr_(),
959  originalPointPriority_(pointPriority),
960  pointPriority_(),
961  minEdgeLen_(),
962  faceFilterFactor_()
963 {
965 }
966 
967 
968 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
969 
971 {}
972 
973 
974 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
975 
977 {
978  minEdgeLen_.resize(mesh_.nEdges(), minLen());
979  faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());
980 
981  return filterFacesLoop(nOriginalBadFaces);
982 }
983 
984 
986 {
987  minEdgeLen_.resize(mesh_.nEdges(), minLen());
988  faceFilterFactor_.resize(mesh_.nFaces(), initialFaceLengthFactor());
989 
990  forAll(faceFilterFactor_, fI)
991  {
992  if (!fSet.found(fI))
993  {
994  faceFilterFactor_[fI] = -1;
995  }
996  }
997 
998  return filterFacesLoop(0);
999 }
1000 
1001 
1002 Foam::label Foam::polyMeshFilter::filterEdges
1004  const label nOriginalBadFaces
1005 )
1006 {
1007  label nBadFaces = labelMax/2;
1008  label nPreviousBadFaces = labelMax;
1009  label nOuterIterations = 0;
1010 
1011  minEdgeLen_.resize(mesh_.nEdges(), minLen());
1012  faceFilterFactor_.resize(0);
1013 
1014  labelList pointErrorCount(mesh_.nPoints(), 0);
1015 
1016  // Main loop
1017  // ~~~~~~~~~
1018  // It tries and do some collapses, checks the resulting mesh and
1019  // 'freezes' some edges (by marking in minEdgeLen) and tries again.
1020  // This will iterate ultimately to the situation where every edge is
1021  // frozen and nothing gets collapsed.
1022  while
1023  (
1024  nOuterIterations < maxIterations()
1025  && nBadFaces > nOriginalBadFaces
1026  && nBadFaces < nPreviousBadFaces
1027  )
1028  {
1029  Info<< nl << "Outer Iteration = " << nOuterIterations++ << nl
1030  << endl;
1031 
1032  printScalarFieldStats("Edge Filter Factor", minEdgeLen_);
1033 
1034  nPreviousBadFaces = nBadFaces;
1035 
1036  // Reset the new mesh to the old mesh
1037  newMeshPtr_ = copyMesh(mesh_);
1038  fvMesh& newMesh = newMeshPtr_();
1039 
1040  scalarField newMeshMinEdgeLen = minEdgeLen_;
1041  pointPriority_.reset(new labelList(originalPointPriority_));
1042 
1043  labelList origToCurrentPointMap(identity(newMesh.nPoints()));
1044 
1045  Info<< incrIndent;
1046 
1047  label nInnerIterations = 0;
1048  label nPrevLocalCollapse = labelMax;
1049 
1050  while (true)
1051  {
1052  Info<< nl
1053  << indent << "Inner iteration = " << nInnerIterations++ << nl
1054  << incrIndent << endl;
1055 
1056  label nLocalCollapse = filterEdges
1057  (
1058  newMesh,
1059  newMeshMinEdgeLen,
1060  origToCurrentPointMap
1061  );
1062 
1063  Info<< decrIndent;
1064 
1065  if
1066  (
1067  nLocalCollapse >= nPrevLocalCollapse
1068  || nLocalCollapse == 0
1069  )
1070  {
1071  Info<< decrIndent;
1072  break;
1073  }
1074  else
1075  {
1076  nPrevLocalCollapse = nLocalCollapse;
1077  }
1078  }
1079 
1080  // Mesh check
1081  // ~~~~~~~~~~~~~~~~~~
1082  // Do not allow collapses in regions of error.
1083  // Updates minEdgeLen, nRelaxedEdges
1084 
1085  if (controlMeshQuality())
1086  {
1087  PackedBoolList isErrorPoint(newMesh.nPoints());
1089  (
1090  newMesh,
1092  isErrorPoint
1093  );
1094 
1095  Info<< nl << " Number of bad faces : " << nBadFaces << nl
1096  << " Number of marked points : "
1097  << returnReduce(isErrorPoint.count(), sumOp<unsigned int>())
1098  << endl;
1099 
1100  updatePointErrorCount
1101  (
1102  isErrorPoint,
1103  origToCurrentPointMap,
1104  pointErrorCount
1105  );
1106 
1107  checkMeshEdgesAndRelaxEdges
1108  (
1109  newMesh,
1110  origToCurrentPointMap,
1111  isErrorPoint,
1112  pointErrorCount
1113  );
1114  }
1115  else
1116  {
1117  return -1;
1118  }
1119  }
1120 
1121  return nBadFaces;
1122 }
1123 
1124 
1126 {
1127  return newMeshPtr_;
1128 }
1129 
1130 
1133 {
1134  return pointPriority_;
1135 }
1136 
1137 
1138 // ************************************************************************* //
unsigned int count() const
Count number of bits set, O(log(n))
Definition: PackedList.C:55
label end() const
Return end vertex label.
Definition: edgeI.H:92
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: polyMesh.C:1048
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
label filter(const label nOriginalBadFaces)
Filter edges and faces.
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
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:223
A list of face labels.
Definition: faceSet.H:48
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const double e
Elementary charge.
Definition: doubleFloat.H:78
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
autoPtr< mapPolyMesh > makeMesh(autoPtr< fvMesh > &newMesh, const IOobject &io, const polyMesh &mesh, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Create new mesh with old mesh patches.
Unit conversion functions.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
const autoPtr< labelList > & pointPriority() const
Return the new pointPriority list.
const Type & second() const
Return second.
Definition: Pair.H:99
polyMeshFilterSettings(const dictionary &dict)
Construct from dictionary.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
bool hasMotionPoints() const
Has valid preMotionPoints?
Definition: mapPolyMesh.H:613
Does polyTopoChanges to remove edges. Can remove faces due to edge collapse but can not remove cells ...
Definition: edgeCollapser.H:66
const labelList & pointMap() const
Old point map.
Definition: mapPolyMesh.H:390
label collapseEdge(triSurface &surf, const scalar minLen)
Keep collapsing all edges < minLen.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
const Type & first() const
Return first.
Definition: Pair.H:87
const Switch & controlMeshQuality() const
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:50
const autoPtr< fvMesh > & filteredMesh() const
Return reference to the filtered mesh. Does not check if the.
static autoPtr< fvMesh > copyMesh(const fvMesh &mesh)
Return a copy of an fvMesh.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
void writeSettings(Ostream &os) const
Write the settings to a stream.
const dictionary & meshQualityCoeffDict() const
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
virtual void updateMesh(const mapPolyMesh &mpm)
Update the mesh corresponding to given map.
label nPoints
const scalar & initialFaceLengthFactor() const
static const label labelMax
Definition: label.H:62
List< label > labelList
A List of labels.
Definition: labelList.H:56
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
label start() const
Return start vertex label.
Definition: edgeI.H:81
static const char nl
Definition: Ostream.H:262
static label checkMeshQuality(const polyMesh &mesh, const dictionary &meshQualityDict, PackedBoolList &isErrorPoint)
Check mesh and mark points on faces in error.
Definition: edgeCollapser.C:82
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:237
labelList f(nPoints)
~polyMeshFilter()
Destructor.
label nEdges() const
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
label nFaces() const
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:381
Calculates points shared by more than two processor patches or cyclic patches.
Definition: globalPoints.H:100
A bit-packed bool list.
label newPointi
Definition: readKivaGrid.H:501
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
Direct mesh changes based on v1.3 polyTopoChange syntax.
const label & maxIterations() const
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
label nPoints() const
messageStream Info
const labelList & faceMap() const
Old face map.
Definition: mapPolyMesh.H:404
const pointField & preMotionPoints() const
Pre-motion point positions.
Definition: mapPolyMesh.H:607
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const word & system() const
Return system name.
Definition: TimePaths.H:114
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
Class to store the settings for the polyMeshFilter class.
A class for managing temporary objects.
Definition: PtrList.H:54
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:230
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:365
const word & name() const
Return reference to name.
Definition: fvMesh.H:257
static const label labelMin
Definition: label.H:61
Namespace for OpenFOAM.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
int system(const std::string &command)
Execute the specified command.
Definition: POSIX.C:1020
const labelList & reversePointMap() const
Reverse point map.
Definition: mapPolyMesh.H:463
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49