FreeFOAM The Cross-Platform CFD Toolkit
autoLayerDriver.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  All to do with adding cell layers
26 
27 \*----------------------------------------------------------------------------*/
28 
29 #include "autoLayerDriver.H"
30 #include <finiteVolume/fvMesh.H>
31 #include <OpenFOAM/Time.H>
34 #include <OpenFOAM/pointFields.H>
37 #include <meshTools/pointSet.H>
38 #include <meshTools/faceSet.H>
39 #include <meshTools/cellSet.H>
41 #include <OpenFOAM/mapPolyMesh.H>
44 #include <OpenFOAM/OFstream.H>
47 #include <OpenFOAM/IOmanip.H>
48 
49 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
50 
51 namespace Foam
52 {
53 
54 defineTypeNameAndDebug(autoLayerDriver, 0);
55 
56 } // End namespace Foam
57 
58 
59 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
60 
61 Foam::label Foam::autoLayerDriver::mergePatchFacesUndo
62 (
63  const scalar minCos,
64  const scalar concaveCos,
65  const dictionary& motionDict
66 )
67 {
68  fvMesh& mesh = meshRefiner_.mesh();
69 
70  // Patch face merging engine
71  combineFaces faceCombiner(mesh, true);
72 
73  // Pick up all candidate cells on boundary
74  labelHashSet boundaryCells(mesh.nFaces()-mesh.nInternalFaces());
75 
76  {
77  labelList patchIDs(meshRefiner_.meshedPatches());
78 
79  const polyBoundaryMesh& patches = mesh.boundaryMesh();
80 
81  forAll(patchIDs, i)
82  {
83  label patchI = patchIDs[i];
84 
85  const polyPatch& patch = patches[patchI];
86 
87  if (!patch.coupled())
88  {
89  forAll(patch, i)
90  {
91  boundaryCells.insert(mesh.faceOwner()[patch.start()+i]);
92  }
93  }
94  }
95  }
96 
97  // Get all sets of faces that can be merged
98  labelListList allFaceSets
99  (
100  faceCombiner.getMergeSets
101  (
102  minCos,
103  concaveCos,
104  boundaryCells
105  )
106  );
107 
108  label nFaceSets = returnReduce(allFaceSets.size(), sumOp<label>());
109 
110  Info<< "Merging " << nFaceSets << " sets of faces." << nl << endl;
111 
112  if (nFaceSets > 0)
113  {
114  if (debug)
115  {
116  faceSet allSets(mesh, "allFaceSets", allFaceSets.size());
117  forAll(allFaceSets, setI)
118  {
119  forAll(allFaceSets[setI], i)
120  {
121  allSets.insert(allFaceSets[setI][i]);
122  }
123  }
124  Pout<< "Writing all faces to be merged to set "
125  << allSets.objectPath() << endl;
126  allSets.write();
127  }
128 
129 
130  // Topology changes container
131  polyTopoChange meshMod(mesh);
132 
133  // Merge all faces of a set into the first face of the set.
134  faceCombiner.setRefinement(allFaceSets, meshMod);
135 
136  // Experimental: store data for all the points that have been deleted
137  meshRefiner_.storeData
138  (
139  faceCombiner.savedPointLabels(), // points to store
140  labelList(0), // faces to store
141  labelList(0) // cells to store
142  );
143 
144  // Change the mesh (no inflation)
145  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
146 
147  // Update fields
148  mesh.updateMesh(map);
149 
150  // Move mesh (since morphing does not do this)
151  if (map().hasMotionPoints())
152  {
153  mesh.movePoints(map().preMotionPoints());
154  }
155  else
156  {
157  // Delete mesh volumes.
158  mesh.clearOut();
159  }
160 
161  if (meshRefiner_.overwrite())
162  {
163  mesh.setInstance(meshRefiner_.oldInstance());
164  }
165 
166  faceCombiner.updateMesh(map);
167 
168  meshRefiner_.updateMesh(map, labelList(0));
169 
170 
171  for (label iteration = 0; iteration < 100; iteration++)
172  {
173  Info<< nl
174  << "Undo iteration " << iteration << nl
175  << "----------------" << endl;
176 
177 
178  // Check mesh for errors
179  // ~~~~~~~~~~~~~~~~~~~~~
180 
181  faceSet errorFaces
182  (
183  mesh,
184  "errorFaces",
185  mesh.nFaces()-mesh.nInternalFaces()
186  );
187  bool hasErrors = motionSmoother::checkMesh
188  (
189  false, // report
190  mesh,
191  motionDict,
192  errorFaces
193  );
194  //if (checkEdgeConnectivity)
195  //{
196  // Info<< "Checking edge-face connectivity (duplicate faces"
197  // << " or non-consecutive shared vertices)" << endl;
198  //
199  // label nOldSize = errorFaces.size();
200  //
201  // hasErrors =
202  // mesh.checkFaceFaces
203  // (
204  // false,
205  // &errorFaces
206  // )
207  // || hasErrors;
208  //
209  // Info<< "Detected additional "
210  // << returnReduce(errorFaces.size()-nOldSize, sumOp<label>())
211  // << " faces with illegal face-face connectivity"
212  // << endl;
213  //}
214 
215  if (!hasErrors)
216  {
217  break;
218  }
219 
220 
221  if (debug)
222  {
223  Pout<< "Writing all faces in error to faceSet "
224  << errorFaces.objectPath() << nl << endl;
225  errorFaces.write();
226  }
227 
228 
229  // Check any master cells for using any of the error faces
230  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
231 
232  DynamicList<label> mastersToRestore(allFaceSets.size());
233 
234  forAll(allFaceSets, setI)
235  {
236  label masterFaceI = faceCombiner.masterFace()[setI];
237 
238  if (masterFaceI != -1)
239  {
240  label masterCellII = mesh.faceOwner()[masterFaceI];
241 
242  const cell& cFaces = mesh.cells()[masterCellII];
243 
244  forAll(cFaces, i)
245  {
246  if (errorFaces.found(cFaces[i]))
247  {
248  mastersToRestore.append(masterFaceI);
249  break;
250  }
251  }
252  }
253  }
254  mastersToRestore.shrink();
255 
256  label nRestore = returnReduce
257  (
258  mastersToRestore.size(),
259  sumOp<label>()
260  );
261 
262  Info<< "Masters that need to be restored:"
263  << nRestore << endl;
264 
265  if (debug)
266  {
267  faceSet restoreSet(mesh, "mastersToRestore", mastersToRestore);
268  Pout<< "Writing all " << mastersToRestore.size()
269  << " masterfaces to be restored to set "
270  << restoreSet.objectPath() << endl;
271  restoreSet.write();
272  }
273 
274 
275  if (nRestore == 0)
276  {
277  break;
278  }
279 
280 
281  // Undo
282  // ~~~~
283 
284  // Topology changes container
285  polyTopoChange meshMod(mesh);
286 
287  // Merge all faces of a set into the first face of the set.
288  // Experimental:mark all points/faces/cells that have been restored.
289  Map<label> restoredPoints(0);
290  Map<label> restoredFaces(0);
291  Map<label> restoredCells(0);
292 
293  faceCombiner.setUnrefinement
294  (
295  mastersToRestore,
296  meshMod,
297  restoredPoints,
298  restoredFaces,
299  restoredCells
300  );
301 
302  // Change the mesh (no inflation)
303  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
304 
305  // Update fields
306  mesh.updateMesh(map);
307 
308  // Move mesh (since morphing does not do this)
309  if (map().hasMotionPoints())
310  {
311  mesh.movePoints(map().preMotionPoints());
312  }
313  else
314  {
315  // Delete mesh volumes.
316  mesh.clearOut();
317  }
318 
319  if (meshRefiner_.overwrite())
320  {
321  mesh.setInstance(meshRefiner_.oldInstance());
322  }
323 
324  faceCombiner.updateMesh(map);
325 
326  // Renumber restore maps
327  inplaceMapKey(map().reversePointMap(), restoredPoints);
328  inplaceMapKey(map().reverseFaceMap(), restoredFaces);
329  inplaceMapKey(map().reverseCellMap(), restoredCells);
330 
331  // Experimental:restore all points/face/cells in maps
332  meshRefiner_.updateMesh
333  (
334  map,
335  labelList(0), // changedFaces
336  restoredPoints,
337  restoredFaces,
338  restoredCells
339  );
340 
341  Info<< endl;
342  }
343 
344  if (debug)
345  {
346  Pout<< "Writing merged-faces mesh to time "
347  << meshRefiner_.timeName() << nl << endl;
348  mesh.write();
349  }
350  }
351  else
352  {
353  Info<< "No faces merged ..." << endl;
354  }
355 
356  return nFaceSets;
357 }
358 
359 
360 // Remove points. pointCanBeDeleted is parallel synchronised.
361 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRemovePoints
362 (
363  removePoints& pointRemover,
364  const boolList& pointCanBeDeleted
365 )
366 {
367  fvMesh& mesh = meshRefiner_.mesh();
368 
369  // Topology changes container
370  polyTopoChange meshMod(mesh);
371 
372  pointRemover.setRefinement(pointCanBeDeleted, meshMod);
373 
374  // Change the mesh (no inflation)
375  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
376 
377  // Update fields
378  mesh.updateMesh(map);
379 
380  // Move mesh (since morphing does not do this)
381  if (map().hasMotionPoints())
382  {
383  mesh.movePoints(map().preMotionPoints());
384  }
385  else
386  {
387  // Delete mesh volumes.
388  mesh.clearOut();
389  }
390 
391  if (meshRefiner_.overwrite())
392  {
393  mesh.setInstance(meshRefiner_.oldInstance());
394  }
395 
396  pointRemover.updateMesh(map);
397  meshRefiner_.updateMesh(map, labelList(0));
398 
399  return map;
400 }
401 
402 
403 // Restore faces (which contain removed points)
404 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoLayerDriver::doRestorePoints
405 (
406  removePoints& pointRemover,
407  const labelList& facesToRestore
408 )
409 {
410  fvMesh& mesh = meshRefiner_.mesh();
411 
412  // Topology changes container
413  polyTopoChange meshMod(mesh);
414 
415  // Determine sets of points and faces to restore
416  labelList localFaces, localPoints;
417  pointRemover.getUnrefimentSet
418  (
419  facesToRestore,
420  localFaces,
421  localPoints
422  );
423 
424  // Undo the changes on the faces that are in error.
425  pointRemover.setUnrefinement
426  (
427  localFaces,
428  localPoints,
429  meshMod
430  );
431 
432  // Change the mesh (no inflation)
433  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
434 
435  // Update fields
436  mesh.updateMesh(map);
437 
438  // Move mesh (since morphing does not do this)
439  if (map().hasMotionPoints())
440  {
441  mesh.movePoints(map().preMotionPoints());
442  }
443  else
444  {
445  // Delete mesh volumes.
446  mesh.clearOut();
447  }
448 
449  if (meshRefiner_.overwrite())
450  {
451  mesh.setInstance(meshRefiner_.oldInstance());
452  }
453 
454  pointRemover.updateMesh(map);
455  meshRefiner_.updateMesh(map, labelList(0));
456 
457  return map;
458 }
459 
460 
461 // Collect all faces that are both in candidateFaces and in the set.
462 // If coupled face also collects the coupled face.
463 Foam::labelList Foam::autoLayerDriver::collectFaces
464 (
465  const labelList& candidateFaces,
466  const labelHashSet& set
467 ) const
468 {
469  const fvMesh& mesh = meshRefiner_.mesh();
470 
471  // Has face been selected?
472  boolList selected(mesh.nFaces(), false);
473 
474  forAll(candidateFaces, i)
475  {
476  label faceI = candidateFaces[i];
477 
478  if (set.found(faceI))
479  {
480  selected[faceI] = true;
481  }
482  }
484  (
485  mesh,
486  selected,
487  orEqOp<bool>(), // combine operator
488  false // separation
489  );
490 
491  labelList selectedFaces(findIndices(selected, true));
492 
493  return selectedFaces;
494 }
495 
496 
497 // Pick up faces of cells of faces in set.
498 Foam::labelList Foam::autoLayerDriver::growFaceCellFace
499 (
500  const labelHashSet& set
501 ) const
502 {
503  const fvMesh& mesh = meshRefiner_.mesh();
504 
505  boolList selected(mesh.nFaces(), false);
506 
507  forAllConstIter(faceSet, set, iter)
508  {
509  label faceI = iter.key();
510 
511  label own = mesh.faceOwner()[faceI];
512 
513  const cell& ownFaces = mesh.cells()[own];
514  forAll(ownFaces, ownFaceI)
515  {
516  selected[ownFaces[ownFaceI]] = true;
517  }
518 
519  if (mesh.isInternalFace(faceI))
520  {
521  label nbr = mesh.faceNeighbour()[faceI];
522 
523  const cell& nbrFaces = mesh.cells()[nbr];
524  forAll(nbrFaces, nbrFaceI)
525  {
526  selected[nbrFaces[nbrFaceI]] = true;
527  }
528  }
529  }
531  (
532  mesh,
533  selected,
534  orEqOp<bool>(), // combine operator
535  false // separation
536  );
537  return findIndices(selected, true);
538 }
539 
540 
541 // Remove points not used by any face or points used by only two faces where
542 // the edges are in line
543 Foam::label Foam::autoLayerDriver::mergeEdgesUndo
544 (
545  const scalar minCos,
546  const dictionary& motionDict
547 )
548 {
549  fvMesh& mesh = meshRefiner_.mesh();
550 
551  Info<< nl
552  << "Merging all points on surface that" << nl
553  << "- are used by only two boundary faces and" << nl
554  << "- make an angle with a cosine of more than " << minCos
555  << "." << nl << endl;
556 
557  // Point removal analysis engine with undo
558  removePoints pointRemover(mesh, true);
559 
560  // Count usage of points
561  boolList pointCanBeDeleted;
562  label nRemove = pointRemover.countPointUsage(minCos, pointCanBeDeleted);
563 
564  if (nRemove > 0)
565  {
566  Info<< "Removing " << nRemove
567  << " straight edge points ..." << nl << endl;
568 
569  // Remove points
570  // ~~~~~~~~~~~~~
571 
572  doRemovePoints(pointRemover, pointCanBeDeleted);
573 
574 
575  for (label iteration = 0; iteration < 100; iteration++)
576  {
577  Info<< nl
578  << "Undo iteration " << iteration << nl
579  << "----------------" << endl;
580 
581 
582  // Check mesh for errors
583  // ~~~~~~~~~~~~~~~~~~~~~
584 
585  faceSet errorFaces
586  (
587  mesh,
588  "errorFaces",
589  mesh.nFaces()-mesh.nInternalFaces()
590  );
591  bool hasErrors = motionSmoother::checkMesh
592  (
593  false, // report
594  mesh,
595  motionDict,
596  errorFaces
597  );
598  //if (checkEdgeConnectivity)
599  //{
600  // Info<< "Checking edge-face connectivity (duplicate faces"
601  // << " or non-consecutive shared vertices)" << endl;
602  //
603  // label nOldSize = errorFaces.size();
604  //
605  // hasErrors =
606  // mesh.checkFaceFaces
607  // (
608  // false,
609  // &errorFaces
610  // )
611  // || hasErrors;
612  //
613  // Info<< "Detected additional "
614  // << returnReduce(errorFaces.size()-nOldSize,sumOp<label>())
615  // << " faces with illegal face-face connectivity"
616  // << endl;
617  //}
618 
619  if (!hasErrors)
620  {
621  break;
622  }
623 
624  if (debug)
625  {
626  Pout<< "**Writing all faces in error to faceSet "
627  << errorFaces.objectPath() << nl << endl;
628  errorFaces.write();
629  }
630 
631  labelList masterErrorFaces
632  (
633  collectFaces
634  (
635  pointRemover.savedFaceLabels(),
636  errorFaces
637  )
638  );
639 
640  label n = returnReduce(masterErrorFaces.size(), sumOp<label>());
641 
642  Info<< "Detected " << n
643  << " error faces on boundaries that have been merged."
644  << " These will be restored to their original faces." << nl
645  << endl;
646 
647  if (n == 0)
648  {
649  if (hasErrors)
650  {
651  Info<< "Detected "
652  << returnReduce(errorFaces.size(), sumOp<label>())
653  << " error faces in mesh."
654  << " Restoring neighbours of faces in error." << nl
655  << endl;
656 
657  labelList expandedErrorFaces
658  (
659  growFaceCellFace
660  (
661  errorFaces
662  )
663  );
664 
665  doRestorePoints(pointRemover, expandedErrorFaces);
666  }
667 
668  break;
669  }
670 
671  doRestorePoints(pointRemover, masterErrorFaces);
672  }
673 
674  if (debug)
675  {
676  Pout<< "Writing merged-edges mesh to time "
677  << meshRefiner_.timeName() << nl << endl;
678  mesh.write();
679  }
680  }
681  else
682  {
683  Info<< "No straight edges simplified and no points removed ..." << endl;
684  }
685 
686  return nRemove;
687 }
688 
689 
690 // For debugging: Dump displacement to .obj files
691 void Foam::autoLayerDriver::dumpDisplacement
692 (
693  const fileName& prefix,
694  const indirectPrimitivePatch& pp,
695  const vectorField& patchDisp,
696  const List<extrudeMode>& extrudeStatus
697 )
698 {
699  OFstream dispStr(prefix + "_disp.obj");
700  Info<< "Writing all displacements to " << dispStr.name() << nl << endl;
701 
702  label vertI = 0;
703 
704  forAll(patchDisp, patchPointI)
705  {
706  const point& pt = pp.localPoints()[patchPointI];
707 
708  meshTools::writeOBJ(dispStr, pt); vertI++;
709  meshTools::writeOBJ(dispStr, pt + patchDisp[patchPointI]); vertI++;
710 
711  dispStr << "l " << vertI-1 << ' ' << vertI << nl;
712  }
713 
714 
715  OFstream illStr(prefix + "_illegal.obj");
716  Info<< "Writing invalid displacements to " << illStr.name() << nl << endl;
717 
718  vertI = 0;
719 
720  forAll(patchDisp, patchPointI)
721  {
722  if (extrudeStatus[patchPointI] != EXTRUDE)
723  {
724  const point& pt = pp.localPoints()[patchPointI];
725 
726  meshTools::writeOBJ(illStr, pt); vertI++;
727  meshTools::writeOBJ(illStr, pt + patchDisp[patchPointI]); vertI++;
728 
729  illStr << "l " << vertI-1 << ' ' << vertI << nl;
730  }
731  }
732 }
733 
734 
735 // Check that primitivePatch is not multiply connected. Collect non-manifold
736 // points in pointSet.
737 void Foam::autoLayerDriver::checkManifold
738 (
739  const indirectPrimitivePatch& fp,
740  pointSet& nonManifoldPoints
741 )
742 {
743  // Check for non-manifold points (surface pinched at point)
744  fp.checkPointManifold(false, &nonManifoldPoints);
745 
746  // Check for edge-faces (surface pinched at edge)
747  const labelListList& edgeFaces = fp.edgeFaces();
748 
749  forAll(edgeFaces, edgeI)
750  {
751  const labelList& eFaces = edgeFaces[edgeI];
752 
753  if (eFaces.size() > 2)
754  {
755  const edge& e = fp.edges()[edgeI];
756 
757  nonManifoldPoints.insert(fp.meshPoints()[e[0]]);
758  nonManifoldPoints.insert(fp.meshPoints()[e[1]]);
759  }
760  }
761 }
762 
763 
764 void Foam::autoLayerDriver::checkMeshManifold() const
765 {
766  const fvMesh& mesh = meshRefiner_.mesh();
767 
768  Info<< nl << "Checking mesh manifoldness ..." << endl;
769 
770  // Get all outside faces
771  labelList outsideFaces(mesh.nFaces() - mesh.nInternalFaces());
772 
773  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
774  {
775  outsideFaces[faceI - mesh.nInternalFaces()] = faceI;
776  }
777 
778  pointSet nonManifoldPoints
779  (
780  mesh,
781  "nonManifoldPoints",
782  mesh.nPoints() / 100
783  );
784 
785  // Build primitivePatch out of faces and check it for problems.
786  checkManifold
787  (
789  (
790  IndirectList<face>(mesh.faces(), outsideFaces),
791  mesh.points()
792  ),
793  nonManifoldPoints
794  );
795 
796  label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
797 
798  if (nNonManif > 0)
799  {
800  Info<< "Outside of mesh is multiply connected across edges or"
801  << " points." << nl
802  << "This is not a fatal error but might cause some unexpected"
803  << " behaviour." << nl
804  << "Writing " << nNonManif
805  << " points where this happens to pointSet "
806  << nonManifoldPoints.name() << endl;
807 
808  nonManifoldPoints.write();
809  }
810  Info<< endl;
811 }
812 
813 
814 
815 // Unset extrusion on point. Returns true if anything unset.
816 bool Foam::autoLayerDriver::unmarkExtrusion
817 (
818  const label patchPointI,
819  pointField& patchDisp,
820  labelList& patchNLayers,
821  List<extrudeMode>& extrudeStatus
822 )
823 {
824  if (extrudeStatus[patchPointI] == EXTRUDE)
825  {
826  extrudeStatus[patchPointI] = NOEXTRUDE;
827  patchNLayers[patchPointI] = 0;
828  patchDisp[patchPointI] = vector::zero;
829  return true;
830  }
831  else if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
832  {
833  extrudeStatus[patchPointI] = NOEXTRUDE;
834  patchNLayers[patchPointI] = 0;
835  patchDisp[patchPointI] = vector::zero;
836  return true;
837  }
838  else
839  {
840  return false;
841  }
842 }
843 
844 
845 // Unset extrusion on face. Returns true if anything unset.
846 bool Foam::autoLayerDriver::unmarkExtrusion
847 (
848  const face& localFace,
849  pointField& patchDisp,
850  labelList& patchNLayers,
851  List<extrudeMode>& extrudeStatus
852 )
853 {
854  bool unextruded = false;
855 
856  forAll(localFace, fp)
857  {
858  if
859  (
860  unmarkExtrusion
861  (
862  localFace[fp],
863  patchDisp,
864  patchNLayers,
865  extrudeStatus
866  )
867  )
868  {
869  unextruded = true;
870  }
871  }
872  return unextruded;
873 }
874 
875 
876 // No extrusion at non-manifold points.
877 void Foam::autoLayerDriver::handleNonManifolds
878 (
879  const indirectPrimitivePatch& pp,
880  const labelList& meshEdges,
881  pointField& patchDisp,
882  labelList& patchNLayers,
883  List<extrudeMode>& extrudeStatus
884 ) const
885 {
886  const fvMesh& mesh = meshRefiner_.mesh();
887 
888  Info<< nl << "Handling non-manifold points ..." << endl;
889 
890  // Detect non-manifold points
891  Info<< nl << "Checking patch manifoldness ..." << endl;
892 
893  pointSet nonManifoldPoints(mesh, "nonManifoldPoints", pp.nPoints());
894 
895  // 1. Local check
896  checkManifold(pp, nonManifoldPoints);
897 
898  label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
899 
900  Info<< "Outside of local patch is multiply connected across edges or"
901  << " points at " << nNonManif << " points." << endl;
902 
903 
904  const labelList& meshPoints = pp.meshPoints();
905 
906 
907  // 2. Parallel check
908  // For now disable extrusion at any shared edges.
909  const labelHashSet sharedEdgeSet(mesh.globalData().sharedEdgeLabels());
910 
911  forAll(pp.edges(), edgeI)
912  {
913  if (sharedEdgeSet.found(meshEdges[edgeI]))
914  {
915  const edge& e = mesh.edges()[meshEdges[edgeI]];
916 
917  Pout<< "Disabling extrusion at edge "
918  << mesh.points()[e[0]] << mesh.points()[e[1]]
919  << " since it is non-manifold across coupled patches."
920  << endl;
921 
922  nonManifoldPoints.insert(e[0]);
923  nonManifoldPoints.insert(e[1]);
924  }
925  }
926 
927  // 3b. extrusion can produce multiple faces between 2 cells
928  // across processor boundary
929  // This occurs when a coupled face shares more than 1 edge with a
930  // non-coupled boundary face.
931  // This is now correctly handled by addPatchCellLayer in that it
932  // extrudes a single face from the stringed up edges.
933 
934 
935  nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
936 
937  if (nNonManif > 0)
938  {
939  Info<< "Outside of patches is multiply connected across edges or"
940  << " points at " << nNonManif << " points." << nl
941  << "Writing " << nNonManif
942  << " points where this happens to pointSet "
943  << nonManifoldPoints.name() << nl
944  << "and setting layer thickness to zero on these points."
945  << endl;
946 
947  nonManifoldPoints.write();
948 
949  forAll(meshPoints, patchPointI)
950  {
951  if (nonManifoldPoints.found(meshPoints[patchPointI]))
952  {
953  unmarkExtrusion
954  (
955  patchPointI,
956  patchDisp,
957  patchNLayers,
958  extrudeStatus
959  );
960  }
961  }
962  }
963 
964  Info<< "Set displacement to zero for all " << nNonManif
965  << " non-manifold points" << endl;
966 }
967 
968 
969 // Parallel feature edge detection. Assumes non-manifold edges already handled.
970 void Foam::autoLayerDriver::handleFeatureAngle
971 (
972  const indirectPrimitivePatch& pp,
973  const labelList& meshEdges,
974  const scalar minCos,
975  pointField& patchDisp,
976  labelList& patchNLayers,
977  List<extrudeMode>& extrudeStatus
978 ) const
979 {
980  const fvMesh& mesh = meshRefiner_.mesh();
981 
982  Info<< nl << "Handling feature edges ..." << endl;
983 
984  if (minCos < 1-SMALL)
985  {
986  // Normal component of normals of connected faces.
987  vectorField edgeNormal(mesh.nEdges(), wallPoint::greatPoint);
988 
989  const labelListList& edgeFaces = pp.edgeFaces();
990 
991  forAll(edgeFaces, edgeI)
992  {
993  const labelList& eFaces = pp.edgeFaces()[edgeI];
994 
995  label meshEdgeI = meshEdges[edgeI];
996 
997  forAll(eFaces, i)
998  {
999  nomalsCombine()
1000  (
1001  edgeNormal[meshEdgeI],
1002  pp.faceNormals()[eFaces[i]]
1003  );
1004  }
1005  }
1006 
1008  (
1009  mesh,
1010  edgeNormal,
1011  nomalsCombine(),
1012  wallPoint::greatPoint, // null value
1013  false // no separation
1014  );
1015 
1016  label vertI = 0;
1017  autoPtr<OFstream> str;
1018  if (debug)
1019  {
1020  str.reset(new OFstream(mesh.time().path()/"featureEdges.obj"));
1021  Info<< "Writing feature edges to " << str().name() << endl;
1022  }
1023 
1024  label nFeats = 0;
1025 
1026  // Now on coupled edges the edgeNormal will have been truncated and
1027  // only be still be the old value where two faces have the same normal
1028  forAll(edgeFaces, edgeI)
1029  {
1030  const labelList& eFaces = pp.edgeFaces()[edgeI];
1031 
1032  label meshEdgeI = meshEdges[edgeI];
1033 
1034  const vector& n = edgeNormal[meshEdgeI];
1035 
1036  if (n != wallPoint::greatPoint)
1037  {
1038  scalar cos = n & pp.faceNormals()[eFaces[0]];
1039 
1040  if (cos < minCos)
1041  {
1042  const edge& e = pp.edges()[edgeI];
1043 
1044  unmarkExtrusion
1045  (
1046  e[0],
1047  patchDisp,
1048  patchNLayers,
1049  extrudeStatus
1050  );
1051  unmarkExtrusion
1052  (
1053  e[1],
1054  patchDisp,
1055  patchNLayers,
1056  extrudeStatus
1057  );
1058 
1059  nFeats++;
1060 
1061  if (str.valid())
1062  {
1063  meshTools::writeOBJ(str(), pp.localPoints()[e[0]]);
1064  vertI++;
1065  meshTools::writeOBJ(str(), pp.localPoints()[e[1]]);
1066  vertI++;
1067  str()<< "l " << vertI-1 << ' ' << vertI << nl;
1068  }
1069  }
1070  }
1071  }
1072 
1073  Info<< "Set displacement to zero for points on "
1074  << returnReduce(nFeats, sumOp<label>())
1075  << " feature edges" << endl;
1076  }
1077 }
1078 
1079 
1080 // No extrusion on cells with warped faces. Calculates the thickness of the
1081 // layer and compares it to the space the warped face takes up. Disables
1082 // extrusion if layer thickness is more than faceRatio of the thickness of
1083 // the face.
1084 void Foam::autoLayerDriver::handleWarpedFaces
1085 (
1086  const indirectPrimitivePatch& pp,
1087  const scalar faceRatio,
1088  const scalar edge0Len,
1089  const labelList& cellLevel,
1090  pointField& patchDisp,
1091  labelList& patchNLayers,
1092  List<extrudeMode>& extrudeStatus
1093 ) const
1094 {
1095  const fvMesh& mesh = meshRefiner_.mesh();
1096 
1097  Info<< nl << "Handling cells with warped patch faces ..." << nl;
1098 
1099  const pointField& points = mesh.points();
1100 
1101  label nWarpedFaces = 0;
1102 
1103  forAll(pp, i)
1104  {
1105  const face& f = pp[i];
1106 
1107  if (f.size() > 3)
1108  {
1109  label faceI = pp.addressing()[i];
1110 
1111  label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
1112  scalar edgeLen = edge0Len/(1<<ownLevel);
1113 
1114  // Normal distance to face centre plane
1115  const point& fc = mesh.faceCentres()[faceI];
1116  const vector& fn = pp.faceNormals()[i];
1117 
1118  scalarField vProj(f.size());
1119 
1120  forAll(f, fp)
1121  {
1122  vector n = points[f[fp]] - fc;
1123  vProj[fp] = (n & fn);
1124  }
1125 
1126  // Get normal 'span' of face
1127  scalar minVal = min(vProj);
1128  scalar maxVal = max(vProj);
1129 
1130  if ((maxVal - minVal) > faceRatio * edgeLen)
1131  {
1132  if
1133  (
1134  unmarkExtrusion
1135  (
1136  pp.localFaces()[i],
1137  patchDisp,
1138  patchNLayers,
1139  extrudeStatus
1140  )
1141  )
1142  {
1143  nWarpedFaces++;
1144  }
1145  }
1146  }
1147  }
1148 
1149  Info<< "Set displacement to zero on "
1150  << returnReduce(nWarpedFaces, sumOp<label>())
1151  << " warped faces since layer would be > " << faceRatio
1152  << " of the size of the bounding box." << endl;
1153 }
1154 
1155 
1156 
1159 //void Foam::autoLayerDriver::handleMultiplePatchFaces
1160 //(
1161 // const indirectPrimitivePatch& pp,
1162 // pointField& patchDisp,
1163 // labelList& patchNLayers,
1164 // List<extrudeMode>& extrudeStatus
1165 //) const
1166 //{
1167 // const fvMesh& mesh = meshRefiner_.mesh();
1168 //
1169 // Info<< nl << "Handling cells with multiple patch faces ..." << nl;
1170 //
1171 // const labelListList& pointFaces = pp.pointFaces();
1172 //
1173 // // Cells that should not get an extrusion layer
1174 // cellSet multiPatchCells(mesh, "multiPatchCells", pp.size());
1175 //
1176 // // Detect points that use multiple faces on same cell.
1177 // forAll(pointFaces, patchPointI)
1178 // {
1179 // const labelList& pFaces = pointFaces[patchPointI];
1180 //
1181 // labelHashSet pointCells(pFaces.size());
1182 //
1183 // forAll(pFaces, i)
1184 // {
1185 // label cellI = mesh.faceOwner()[pp.addressing()[pFaces[i]]];
1186 //
1187 // if (!pointCells.insert(cellI))
1188 // {
1189 // // Second or more occurrence of cell so cell has two or more
1190 // // pp faces connected to this point.
1191 // multiPatchCells.insert(cellI);
1192 // }
1193 // }
1194 // }
1195 //
1196 // label nMultiPatchCells = returnReduce
1197 // (
1198 // multiPatchCells.size(),
1199 // sumOp<label>()
1200 // );
1201 //
1202 // Info<< "Detected " << nMultiPatchCells
1203 // << " cells with multiple (connected) patch faces." << endl;
1204 //
1205 // label nChanged = 0;
1206 //
1207 // if (nMultiPatchCells > 0)
1208 // {
1209 // Info<< "Writing " << nMultiPatchCells
1210 // << " cells with multiple (connected) patch faces to cellSet "
1211 // << multiPatchCells.objectPath() << endl;
1212 // multiPatchCells.write();
1213 //
1214 //
1215 // // Go through all points and remove extrusion on any cell in
1216 // // multiPatchCells
1217 // // (has to be done in separate loop since having one point on
1218 // // multipatches has to reset extrusion on all points of cell)
1219 //
1220 // forAll(pointFaces, patchPointI)
1221 // {
1222 // if (extrudeStatus[patchPointI] != NOEXTRUDE)
1223 // {
1224 // const labelList& pFaces = pointFaces[patchPointI];
1225 //
1226 // forAll(pFaces, i)
1227 // {
1228 // label cellI =
1229 // mesh.faceOwner()[pp.addressing()[pFaces[i]]];
1230 //
1231 // if (multiPatchCells.found(cellI))
1232 // {
1233 // if
1234 // (
1235 // unmarkExtrusion
1236 // (
1237 // patchPointI,
1238 // patchDisp,
1239 // patchNLayers,
1240 // extrudeStatus
1241 // )
1242 // )
1243 // {
1244 // nChanged++;
1245 // }
1246 // }
1247 // }
1248 // }
1249 // }
1250 //
1251 // reduce(nChanged, sumOp<label>());
1252 // }
1253 //
1254 // Info<< "Prevented extrusion on " << nChanged
1255 // << " points due to multiple patch faces." << nl << endl;
1256 //}
1257 
1258 
1259 // No extrusion on faces with differing number of layers for points
1260 void Foam::autoLayerDriver::setNumLayers
1261 (
1262  const labelList& patchToNLayers,
1263  const labelList& patchIDs,
1264  const indirectPrimitivePatch& pp,
1265  pointField& patchDisp,
1266  labelList& patchNLayers,
1267  List<extrudeMode>& extrudeStatus
1268 ) const
1269 {
1270  const fvMesh& mesh = meshRefiner_.mesh();
1271 
1272  Info<< nl << "Handling points with inconsistent layer specification ..."
1273  << endl;
1274 
1275  // Get for every point (really only nessecary on patch external points)
1276  // the max and min of any patch faces using it.
1277  labelList maxLayers(patchNLayers.size(), labelMin);
1278  labelList minLayers(patchNLayers.size(), labelMax);
1279 
1280  forAll(patchIDs, i)
1281  {
1282  label patchI = patchIDs[i];
1283 
1284  const labelList& meshPoints = mesh.boundaryMesh()[patchI].meshPoints();
1285 
1286  label wantedLayers = patchToNLayers[patchI];
1287 
1288  forAll(meshPoints, patchPointI)
1289  {
1290  label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
1291 
1292  maxLayers[ppPointI] = max(wantedLayers, maxLayers[ppPointI]);
1293  minLayers[ppPointI] = min(wantedLayers, minLayers[ppPointI]);
1294  }
1295  }
1296 
1298  (
1299  mesh,
1300  pp.meshPoints(),
1301  maxLayers,
1302  maxEqOp<label>(),
1303  labelMin, // null value
1304  false // no separation
1305  );
1307  (
1308  mesh,
1309  pp.meshPoints(),
1310  minLayers,
1311  minEqOp<label>(),
1312  labelMax, // null value
1313  false // no separation
1314  );
1315 
1316  // Unmark any point with different min and max
1317  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1318 
1319  //label nConflicts = 0;
1320 
1321  forAll(maxLayers, i)
1322  {
1323  if (maxLayers[i] == labelMin || minLayers[i] == labelMax)
1324  {
1325  FatalErrorIn("setNumLayers(..)")
1326  << "Patchpoint:" << i << " coord:" << pp.localPoints()[i]
1327  << " maxLayers:" << maxLayers
1328  << " minLayers:" << minLayers
1329  << abort(FatalError);
1330  }
1331  else if (maxLayers[i] == minLayers[i])
1332  {
1333  // Ok setting.
1334  patchNLayers[i] = maxLayers[i];
1335  }
1336  else
1337  {
1338  // Inconsistent num layers between patch faces using point
1339  //if
1340  //(
1341  // unmarkExtrusion
1342  // (
1343  // i,
1344  // patchDisp,
1345  // patchNLayers,
1346  // extrudeStatus
1347  // )
1348  //)
1349  //{
1350  // nConflicts++;
1351  //}
1352  patchNLayers[i] = maxLayers[i];
1353  }
1354  }
1355 
1356  //reduce(nConflicts, sumOp<label>());
1357  //
1358  //Info<< "Set displacement to zero for " << nConflicts
1359  // << " points due to points being on multiple regions"
1360  // << " with inconsistent nLayers specification." << endl;
1361 }
1362 
1363 
1364 // Grow no-extrusion layer
1365 void Foam::autoLayerDriver::growNoExtrusion
1366 (
1367  const indirectPrimitivePatch& pp,
1368  pointField& patchDisp,
1369  labelList& patchNLayers,
1370  List<extrudeMode>& extrudeStatus
1371 )
1372 {
1373  Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
1374 
1375  List<extrudeMode> grownExtrudeStatus(extrudeStatus);
1376 
1377  const faceList& localFaces = pp.localFaces();
1378 
1379  label nGrown = 0;
1380 
1381  forAll(localFaces, faceI)
1382  {
1383  const face& f = localFaces[faceI];
1384 
1385  bool hasSqueeze = false;
1386  forAll(f, fp)
1387  {
1388  if (extrudeStatus[f[fp]] == NOEXTRUDE)
1389  {
1390  hasSqueeze = true;
1391  break;
1392  }
1393  }
1394 
1395  if (hasSqueeze)
1396  {
1397  // Squeeze all points of face
1398  forAll(f, fp)
1399  {
1400  if
1401  (
1402  extrudeStatus[f[fp]] == NOEXTRUDE
1403  && grownExtrudeStatus[f[fp]] != NOEXTRUDE
1404  )
1405  {
1406  grownExtrudeStatus[f[fp]] = NOEXTRUDE;
1407  nGrown++;
1408  }
1409  }
1410  }
1411  }
1412 
1413  extrudeStatus.transfer(grownExtrudeStatus);
1414 
1415  forAll(extrudeStatus, patchPointI)
1416  {
1417  if (extrudeStatus[patchPointI] == NOEXTRUDE)
1418  {
1419  patchDisp[patchPointI] = vector::zero;
1420  patchNLayers[patchPointI] = 0;
1421  }
1422  }
1423 
1424  reduce(nGrown, sumOp<label>());
1425 
1426  Info<< "Set displacement to zero for an additional " << nGrown
1427  << " points." << endl;
1428 }
1429 
1430 
1431 void Foam::autoLayerDriver::calculateLayerThickness
1432 (
1433  const indirectPrimitivePatch& pp,
1434  const labelList& patchIDs,
1435  const scalarField& patchExpansionRatio,
1436 
1437  const bool relativeSizes,
1438  const scalarField& patchFinalLayerThickness,
1439  const scalarField& patchMinThickness,
1440 
1441  const labelList& cellLevel,
1442  const labelList& patchNLayers,
1443  const scalar edge0Len,
1444 
1445  scalarField& thickness,
1446  scalarField& minThickness,
1447  scalarField& expansionRatio
1448 ) const
1449 {
1450  const fvMesh& mesh = meshRefiner_.mesh();
1451  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1452 
1453 
1454  // Rework patch-wise layer parameters into minimum per point
1455  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1456 
1457  // Reuse input fields
1458  expansionRatio.setSize(pp.nPoints());
1459  expansionRatio = GREAT;
1460  thickness.setSize(pp.nPoints());
1461  thickness = GREAT;
1462  minThickness.setSize(pp.nPoints());
1463  minThickness = GREAT;
1464 
1465  forAll(patchIDs, i)
1466  {
1467  label patchI = patchIDs[i];
1468 
1469  const labelList& meshPoints = patches[patchI].meshPoints();
1470 
1471  forAll(meshPoints, patchPointI)
1472  {
1473  label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
1474 
1475  expansionRatio[ppPointI] = min
1476  (
1477  expansionRatio[ppPointI],
1478  patchExpansionRatio[patchI]
1479  );
1480  thickness[ppPointI] = min
1481  (
1482  thickness[ppPointI],
1483  patchFinalLayerThickness[patchI]
1484  );
1485  minThickness[ppPointI] = min
1486  (
1487  minThickness[ppPointI],
1488  patchMinThickness[patchI]
1489  );
1490  }
1491  }
1492 
1494  (
1495  mesh,
1496  pp.meshPoints(),
1497  expansionRatio,
1498  minEqOp<scalar>(),
1499  GREAT, // null value
1500  false // no separation
1501  );
1503  (
1504  mesh,
1505  pp.meshPoints(),
1506  thickness,
1507  minEqOp<scalar>(),
1508  GREAT, // null value
1509  false // no separation
1510  );
1512  (
1513  mesh,
1514  pp.meshPoints(),
1515  minThickness,
1516  minEqOp<scalar>(),
1517  GREAT, // null value
1518  false // no separation
1519  );
1520 
1521 
1522  // Now the thicknesses are set according to the minimum of connected
1523  // patches.
1524 
1525 
1526  // Rework relative thickness into absolute
1527  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1528  // by multiplying with the internal cell size.
1529 
1530  if (relativeSizes)
1531  {
1532  if (min(patchMinThickness) < 0 || max(patchMinThickness) > 2)
1533  {
1534  FatalErrorIn("calculateLayerThickness(..)")
1535  << "Thickness should be factor of local undistorted cell size."
1536  << " Valid values are [0..2]." << nl
1537  << " minThickness:" << patchMinThickness
1538  << exit(FatalError);
1539  }
1540 
1541 
1542  // Determine per point the max cell level of connected cells
1543  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1544 
1545  labelList maxPointLevel(pp.nPoints(), labelMin);
1546 
1547  forAll(pp, i)
1548  {
1549  label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
1550 
1551  const face& f = pp.localFaces()[i];
1552 
1553  forAll(f, fp)
1554  {
1555  maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
1556  }
1557  }
1558 
1560  (
1561  mesh,
1562  pp.meshPoints(),
1563  maxPointLevel,
1564  maxEqOp<label>(),
1565  labelMin, // null value
1566  false // no separation
1567  );
1568 
1569 
1570  forAll(maxPointLevel, pointI)
1571  {
1572  // Find undistorted edge size for this level.
1573  scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
1574  thickness[pointI] *= edgeLen;
1575  minThickness[pointI] *= edgeLen;
1576  }
1577  }
1578 
1579 
1580 
1581  // Rework thickness (of final layer) into overall thickness of all layers
1582  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1583 
1584  forAll(thickness, pointI)
1585  {
1586  // Calculate layer thickness based on expansion ratio
1587  // and final layer height
1588  if (expansionRatio[pointI] == 1)
1589  {
1590  thickness[pointI] *= patchNLayers[pointI];
1591  }
1592  else
1593  {
1594 
1595  scalar invExpansion = 1.0 / expansionRatio[pointI];
1596  label nLay = patchNLayers[pointI];
1597  thickness[pointI] *=
1598  (1.0 - pow(invExpansion, nLay))
1599  / (1.0 - invExpansion);
1600  }
1601  }
1602 
1603 
1604  //Info<< "calculateLayerThickness : min:" << gMin(thickness)
1605  // << " max:" << gMax(thickness) << endl;
1606 }
1607 
1608 
1609 // Synchronize displacement among coupled patches.
1610 void Foam::autoLayerDriver::syncPatchDisplacement
1611 (
1612  const motionSmoother& meshMover,
1613  const scalarField& minThickness,
1614  pointField& patchDisp,
1615  labelList& patchNLayers,
1616  List<extrudeMode>& extrudeStatus
1617 ) const
1618 {
1619  const fvMesh& mesh = meshRefiner_.mesh();
1620  const labelList& meshPoints = meshMover.patch().meshPoints();
1621 
1622  label nChangedTotal = 0;
1623 
1624  while (true)
1625  {
1626  label nChanged = 0;
1627 
1628  // Sync displacement (by taking min)
1630  (
1631  mesh,
1632  meshPoints,
1633  patchDisp,
1634  minEqOp<vector>(),
1635  wallPoint::greatPoint, // null value
1636  false // no separation
1637  );
1638 
1639  // Unmark if displacement too small
1640  forAll(patchDisp, i)
1641  {
1642  if (mag(patchDisp[i]) < minThickness[i])
1643  {
1644  if
1645  (
1646  unmarkExtrusion
1647  (
1648  i,
1649  patchDisp,
1650  patchNLayers,
1651  extrudeStatus
1652  )
1653  )
1654  {
1655  nChanged++;
1656  }
1657  }
1658  }
1659 
1660  labelList syncPatchNLayers(patchNLayers);
1661 
1663  (
1664  mesh,
1665  meshPoints,
1666  syncPatchNLayers,
1667  minEqOp<label>(),
1668  labelMax, // null value
1669  false // no separation
1670  );
1671 
1672  // Reset if differs
1673  forAll(syncPatchNLayers, i)
1674  {
1675  if (syncPatchNLayers[i] != patchNLayers[i])
1676  {
1677  if
1678  (
1679  unmarkExtrusion
1680  (
1681  i,
1682  patchDisp,
1683  patchNLayers,
1684  extrudeStatus
1685  )
1686  )
1687  {
1688  nChanged++;
1689  }
1690  }
1691  }
1692 
1694  (
1695  mesh,
1696  meshPoints,
1697  syncPatchNLayers,
1698  maxEqOp<label>(),
1699  labelMin, // null value
1700  false // no separation
1701  );
1702 
1703  // Reset if differs
1704  forAll(syncPatchNLayers, i)
1705  {
1706  if (syncPatchNLayers[i] != patchNLayers[i])
1707  {
1708  if
1709  (
1710  unmarkExtrusion
1711  (
1712  i,
1713  patchDisp,
1714  patchNLayers,
1715  extrudeStatus
1716  )
1717  )
1718  {
1719  nChanged++;
1720  }
1721  }
1722  }
1723  nChangedTotal += nChanged;
1724 
1725  if (!returnReduce(nChanged, sumOp<label>()))
1726  {
1727  break;
1728  }
1729  }
1730 
1731  Info<< "Prevented extrusion on "
1732  << returnReduce(nChangedTotal, sumOp<label>())
1733  << " coupled patch points during syncPatchDisplacement." << endl;
1734 }
1735 
1736 
1737 // Calculate displacement vector for all patch points. Uses pointNormal.
1738 // Checks that displaced patch point would be visible from all centres
1739 // of the faces using it.
1740 // extrudeStatus is both input and output and gives the status of each
1741 // patch point.
1742 void Foam::autoLayerDriver::getPatchDisplacement
1743 (
1744  const motionSmoother& meshMover,
1745  const scalarField& thickness,
1746  const scalarField& minThickness,
1747  pointField& patchDisp,
1748  labelList& patchNLayers,
1749  List<extrudeMode>& extrudeStatus
1750 ) const
1751 {
1752  Info<< nl << "Determining displacement for added points"
1753  << " according to pointNormal ..." << endl;
1754 
1755  const fvMesh& mesh = meshRefiner_.mesh();
1756  const indirectPrimitivePatch& pp = meshMover.patch();
1757  const vectorField& faceNormals = pp.faceNormals();
1758  const labelListList& pointFaces = pp.pointFaces();
1759  const pointField& localPoints = pp.localPoints();
1760  const labelList& meshPoints = pp.meshPoints();
1761 
1762  // Determine pointNormal
1763  // ~~~~~~~~~~~~~~~~~~~~~
1764 
1765  pointField pointNormals(pp.nPoints(), vector::zero);
1766  {
1767  labelList nPointFaces(pp.nPoints(), 0);
1768 
1769  forAll(faceNormals, faceI)
1770  {
1771  const face& f = pp.localFaces()[faceI];
1772 
1773  forAll(f, fp)
1774  {
1775  pointNormals[f[fp]] += faceNormals[faceI];
1776  nPointFaces[f[fp]] ++;
1777  }
1778  }
1779 
1781  (
1782  mesh,
1783  meshPoints,
1784  pointNormals,
1785  plusEqOp<vector>(),
1786  vector::zero, // null value
1787  false // no separation
1788  );
1789 
1791  (
1792  mesh,
1793  meshPoints,
1794  nPointFaces,
1795  plusEqOp<label>(),
1796  0, // null value
1797  false // no separation
1798  );
1799 
1800  forAll(pointNormals, i)
1801  {
1802  pointNormals[i] /= nPointFaces[i];
1803  }
1804  }
1805 
1806 
1807  // Determine local length scale on patch
1808  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1809 
1810  // Start off from same thickness everywhere (except where no extrusion)
1811  patchDisp = thickness*pointNormals;
1812 
1813  // Check if no extrude possible.
1814  forAll(pointNormals, patchPointI)
1815  {
1816  label meshPointI = pp.meshPoints()[patchPointI];
1817 
1818  if (extrudeStatus[patchPointI] == NOEXTRUDE)
1819  {
1820  // Do not use unmarkExtrusion; forcibly set to zero extrusion.
1821  patchNLayers[patchPointI] = 0;
1822  patchDisp[patchPointI] = vector::zero;
1823  }
1824  else
1825  {
1826  // Get normal
1827  const vector& n = pointNormals[patchPointI];
1828 
1829  if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointI]))
1830  {
1831  Pout<< "No valid normal for point " << meshPointI
1832  << ' ' << pp.points()[meshPointI]
1833  << "; setting displacement to " << patchDisp[patchPointI]
1834  << endl;
1835 
1836  extrudeStatus[patchPointI] = EXTRUDEREMOVE;
1837  }
1838  }
1839  }
1840 
1841  // At illegal points make displacement average of new neighbour positions
1842  forAll(extrudeStatus, patchPointI)
1843  {
1844  if (extrudeStatus[patchPointI] == EXTRUDEREMOVE)
1845  {
1846  point avg(vector::zero);
1847  label nPoints = 0;
1848 
1849  const labelList& pEdges = pp.pointEdges()[patchPointI];
1850 
1851  forAll(pEdges, i)
1852  {
1853  label edgeI = pEdges[i];
1854 
1855  label otherPointI = pp.edges()[edgeI].otherVertex(patchPointI);
1856 
1857  if (extrudeStatus[otherPointI] != NOEXTRUDE)
1858  {
1859  avg += localPoints[otherPointI] + patchDisp[otherPointI];
1860  nPoints++;
1861  }
1862  }
1863 
1864  if (nPoints > 0)
1865  {
1866  Pout<< "Displacement at illegal point "
1867  << localPoints[patchPointI]
1868  << " set to " << (avg / nPoints - localPoints[patchPointI])
1869  << endl;
1870 
1871  patchDisp[patchPointI] =
1872  avg / nPoints
1873  - localPoints[patchPointI];
1874  }
1875  }
1876  }
1877 
1878  // Make sure displacement is equal on both sides of coupled patches.
1879  syncPatchDisplacement
1880  (
1881  meshMover,
1882  minThickness,
1883  patchDisp,
1884  patchNLayers,
1885  extrudeStatus
1886  );
1887 
1888  Info<< endl;
1889 }
1890 
1891 
1892 // Truncates displacement
1893 // - for all patchFaces in the faceset displacement gets set to zero
1894 // - all displacement < minThickness gets set to zero
1895 Foam::label Foam::autoLayerDriver::truncateDisplacement
1896 (
1897  const motionSmoother& meshMover,
1898  const scalarField& minThickness,
1899  const faceSet& illegalPatchFaces,
1900  pointField& patchDisp,
1901  labelList& patchNLayers,
1902  List<extrudeMode>& extrudeStatus
1903 ) const
1904 {
1905  const polyMesh& mesh = meshMover.mesh();
1906  const indirectPrimitivePatch& pp = meshMover.patch();
1907 
1908  label nChanged = 0;
1909 
1910  const Map<label>& meshPointMap = pp.meshPointMap();
1911 
1912  forAllConstIter(faceSet, illegalPatchFaces, iter)
1913  {
1914  label faceI = iter.key();
1915 
1916  if (mesh.isInternalFace(faceI))
1917  {
1918  FatalErrorIn("truncateDisplacement(..)")
1919  << "Faceset " << illegalPatchFaces.name()
1920  << " contains internal face " << faceI << nl
1921  << "It should only contain patch faces" << abort(FatalError);
1922  }
1923 
1924  const face& f = mesh.faces()[faceI];
1925 
1926 
1927  forAll(f, fp)
1928  {
1929  if (meshPointMap.found(f[fp]))
1930  {
1931  label patchPointI = meshPointMap[f[fp]];
1932 
1933  if (extrudeStatus[patchPointI] != NOEXTRUDE)
1934  {
1935  unmarkExtrusion
1936  (
1937  patchPointI,
1938  patchDisp,
1939  patchNLayers,
1940  extrudeStatus
1941  );
1942  nChanged++;
1943  }
1944  }
1945  }
1946  }
1947 
1948  forAll(patchDisp, patchPointI)
1949  {
1950  if (mag(patchDisp[patchPointI]) < minThickness[patchPointI])
1951  {
1952  if
1953  (
1954  unmarkExtrusion
1955  (
1956  patchPointI,
1957  patchDisp,
1958  patchNLayers,
1959  extrudeStatus
1960  )
1961  )
1962  {
1963  nChanged++;
1964  }
1965  }
1966  else if (extrudeStatus[patchPointI] == NOEXTRUDE)
1967  {
1968  // Make sure displacement is 0. Should already be so but ...
1969  patchDisp[patchPointI] = vector::zero;
1970  patchNLayers[patchPointI] = 0;
1971  }
1972  }
1973 
1974 
1975  const faceList& localFaces = pp.localFaces();
1976 
1977  while (true)
1978  {
1979  syncPatchDisplacement
1980  (
1981  meshMover,
1982  minThickness,
1983  patchDisp,
1984  patchNLayers,
1985  extrudeStatus
1986  );
1987 
1988  // Make sure that a face doesn't have two non-consecutive areas
1989  // not extruded (e.g. quad where vertex 0 and 2 are not extruded
1990  // but 1 and 3 are) since this gives topological errors.
1991 
1992  label nPinched = 0;
1993 
1994  forAll(localFaces, i)
1995  {
1996  const face& localF = localFaces[i];
1997 
1998  // Count number of transitions from unsnapped to snapped.
1999  label nTrans = 0;
2000 
2001  extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
2002 
2003  forAll(localF, fp)
2004  {
2005  extrudeMode fpMode = extrudeStatus[localF[fp]];
2006 
2007  if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
2008  {
2009  nTrans++;
2010  }
2011  prevMode = fpMode;
2012  }
2013 
2014  if (nTrans > 1)
2015  {
2016  // Multiple pinches. Reset whole face as unextruded.
2017  if
2018  (
2019  unmarkExtrusion
2020  (
2021  localF,
2022  patchDisp,
2023  patchNLayers,
2024  extrudeStatus
2025  )
2026  )
2027  {
2028  nPinched++;
2029  nChanged++;
2030  }
2031  }
2032  }
2033 
2034  reduce(nPinched, sumOp<label>());
2035 
2036  Info<< "truncateDisplacement : Unextruded " << nPinched
2037  << " faces due to non-consecutive vertices being extruded." << endl;
2038 
2039 
2040  // Make sure that a face has consistent number of layers for all
2041  // its vertices.
2042 
2043  label nDiffering = 0;
2044 
2045  //forAll(localFaces, i)
2046  //{
2047  // const face& localF = localFaces[i];
2048  //
2049  // label numLayers = -1;
2050  //
2051  // forAll(localF, fp)
2052  // {
2053  // if (patchNLayers[localF[fp]] > 0)
2054  // {
2055  // if (numLayers == -1)
2056  // {
2057  // numLayers = patchNLayers[localF[fp]];
2058  // }
2059  // else if (numLayers != patchNLayers[localF[fp]])
2060  // {
2061  // // Differing number of layers
2062  // if
2063  // (
2064  // unmarkExtrusion
2065  // (
2066  // localF,
2067  // patchDisp,
2068  // patchNLayers,
2069  // extrudeStatus
2070  // )
2071  // )
2072  // {
2073  // nDiffering++;
2074  // nChanged++;
2075  // }
2076  // break;
2077  // }
2078  // }
2079  // }
2080  //}
2081  //
2082  //reduce(nDiffering, sumOp<label>());
2083  //
2084  //Info<< "truncateDisplacement : Unextruded " << nDiffering
2085  // << " faces due to having differing number of layers." << endl;
2086 
2087  if (nPinched+nDiffering == 0)
2088  {
2089  break;
2090  }
2091  }
2092 
2093  return nChanged;
2094 }
2095 
2096 
2097 // Setup layer information (at points and faces) to modify mesh topology in
2098 // regions where layer mesh terminates.
2099 void Foam::autoLayerDriver::setupLayerInfoTruncation
2100 (
2101  const motionSmoother& meshMover,
2102  const labelList& patchNLayers,
2103  const List<extrudeMode>& extrudeStatus,
2104  const label nBufferCellsNoExtrude,
2105  labelList& nPatchPointLayers,
2106  labelList& nPatchFaceLayers
2107 ) const
2108 {
2109  Info<< nl << "Setting up information for layer truncation ..." << endl;
2110 
2111  const indirectPrimitivePatch& pp = meshMover.patch();
2112  const polyMesh& mesh = meshMover.mesh();
2113 
2114  if (nBufferCellsNoExtrude < 0)
2115  {
2116  Info<< nl << "Performing no layer truncation."
2117  << " nBufferCellsNoExtrude set to less than 0 ..." << endl;
2118 
2119  // Face layers if any point get extruded
2120  forAll(pp.localFaces(), patchFaceI)
2121  {
2122  const face& f = pp.localFaces()[patchFaceI];
2123 
2124  forAll(f, fp)
2125  {
2126  if (patchNLayers[f[fp]] > 0)
2127  {
2128  nPatchFaceLayers[patchFaceI] = patchNLayers[f[fp]];
2129  break;
2130  }
2131  }
2132  }
2133  nPatchPointLayers = patchNLayers;
2134  }
2135  else
2136  {
2137  // Determine max point layers per face.
2138  labelList maxLevel(pp.size(), 0);
2139 
2140  forAll(pp.localFaces(), patchFaceI)
2141  {
2142  const face& f = pp.localFaces()[patchFaceI];
2143 
2144  // find patch faces where layer terminates (i.e contains extrude
2145  // and noextrude points).
2146 
2147  bool noExtrude = false;
2148  label mLevel = 0;
2149 
2150  forAll(f, fp)
2151  {
2152  if (extrudeStatus[f[fp]] == NOEXTRUDE)
2153  {
2154  noExtrude = true;
2155  }
2156  mLevel = max(mLevel, patchNLayers[f[fp]]);
2157  }
2158 
2159  if (mLevel > 0)
2160  {
2161  // So one of the points is extruded. Check if all are extruded
2162  // or is a mix.
2163 
2164  if (noExtrude)
2165  {
2166  nPatchFaceLayers[patchFaceI] = 1;
2167  maxLevel[patchFaceI] = mLevel;
2168  }
2169  else
2170  {
2171  maxLevel[patchFaceI] = mLevel;
2172  }
2173  }
2174  }
2175 
2176  // We have the seed faces (faces with nPatchFaceLayers != maxLevel)
2177  // Now do a meshwave across the patch where we pick up neighbours
2178  // of seed faces.
2179  // Note: quite inefficient. Could probably be coded better.
2180 
2181  const labelListList& pointFaces = pp.pointFaces();
2182 
2183  label nLevels = gMax(patchNLayers);
2184 
2185  // flag neighbouring patch faces with number of layers to grow
2186  for (label ilevel = 1; ilevel < nLevels; ilevel++)
2187  {
2188  label nBuffer;
2189 
2190  if (ilevel == 1)
2191  {
2192  nBuffer = nBufferCellsNoExtrude - 1;
2193  }
2194  else
2195  {
2196  nBuffer = nBufferCellsNoExtrude;
2197  }
2198 
2199  for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
2200  {
2201  labelList tempCounter(nPatchFaceLayers);
2202 
2203  boolList foundNeighbour(pp.nPoints(), false);
2204 
2205  forAll(pp.meshPoints(), patchPointI)
2206  {
2207  forAll(pointFaces[patchPointI], pointFaceI)
2208  {
2209  label faceI = pointFaces[patchPointI][pointFaceI];
2210 
2211  if
2212  (
2213  nPatchFaceLayers[faceI] != -1
2214  && maxLevel[faceI] > 0
2215  )
2216  {
2217  foundNeighbour[patchPointI] = true;
2218  break;
2219  }
2220  }
2221  }
2222 
2224  (
2225  mesh,
2226  pp.meshPoints(),
2227  foundNeighbour,
2228  orEqOp<bool>(),
2229  false, // null value
2230  false // no separation
2231  );
2232 
2233  forAll(pp.meshPoints(), patchPointI)
2234  {
2235  if (foundNeighbour[patchPointI])
2236  {
2237  forAll(pointFaces[patchPointI], pointFaceI)
2238  {
2239  label faceI = pointFaces[patchPointI][pointFaceI];
2240  if
2241  (
2242  nPatchFaceLayers[faceI] == -1
2243  && maxLevel[faceI] > 0
2244  && ilevel < maxLevel[faceI]
2245  )
2246  {
2247  tempCounter[faceI] = ilevel;
2248  }
2249  }
2250  }
2251  }
2252  nPatchFaceLayers = tempCounter;
2253  }
2254  }
2255 
2256  forAll(pp.localFaces(), patchFaceI)
2257  {
2258  if (nPatchFaceLayers[patchFaceI] == -1)
2259  {
2260  nPatchFaceLayers[patchFaceI] = maxLevel[patchFaceI];
2261  }
2262  }
2263 
2264  forAll(pp.meshPoints(), patchPointI)
2265  {
2266  if (extrudeStatus[patchPointI] != NOEXTRUDE)
2267  {
2268  forAll(pointFaces[patchPointI], pointFaceI)
2269  {
2270  label face = pointFaces[patchPointI][pointFaceI];
2271  nPatchPointLayers[patchPointI] = max
2272  (
2273  nPatchPointLayers[patchPointI],
2274  nPatchFaceLayers[face]
2275  );
2276  }
2277  }
2278  else
2279  {
2280  nPatchPointLayers[patchPointI] = 0;
2281  }
2282  }
2284  (
2285  mesh,
2286  pp.meshPoints(),
2287  nPatchPointLayers,
2288  maxEqOp<label>(),
2289  0, // null value
2290  false // no separation
2291  );
2292  }
2293 }
2294 
2295 
2296 // Does any of the cells use a face from faces?
2297 bool Foam::autoLayerDriver::cellsUseFace
2298 (
2299  const polyMesh& mesh,
2300  const labelList& cellLabels,
2301  const labelHashSet& faces
2302 )
2303 {
2304  forAll(cellLabels, i)
2305  {
2306  const cell& cFaces = mesh.cells()[cellLabels[i]];
2307 
2308  forAll(cFaces, cFaceI)
2309  {
2310  if (faces.found(cFaces[cFaceI]))
2311  {
2312  return true;
2313  }
2314  }
2315  }
2316  return false;
2317 }
2318 
2319 
2320 // Checks the newly added cells and locally unmarks points so they
2321 // will not get extruded next time round. Returns global number of unmarked
2322 // points (0 if all was fine)
2323 Foam::label Foam::autoLayerDriver::checkAndUnmark
2324 (
2325  const addPatchCellLayer& addLayer,
2326  const dictionary& meshQualityDict,
2327  const indirectPrimitivePatch& pp,
2328  const fvMesh& newMesh,
2329 
2330  pointField& patchDisp,
2331  labelList& patchNLayers,
2332  List<extrudeMode>& extrudeStatus
2333 )
2334 {
2335  // Check the resulting mesh for errors
2336  Info<< nl << "Checking mesh with layer ..." << endl;
2337  faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
2338  motionSmoother::checkMesh(false, newMesh, meshQualityDict, wrongFaces);
2339  Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
2340  << " illegal faces"
2341  << " (concave, zero area or negative cell pyramid volume)"
2342  << endl;
2343 
2344  // Undo local extrusion if
2345  // - any of the added cells in error
2346 
2347  label nChanged = 0;
2348 
2349  // Get all cells in the layer.
2350  labelListList addedCells
2351  (
2353  (
2354  newMesh,
2355  addLayer.layerFaces()
2356  )
2357  );
2358 
2359  // Check if any of the faces in error uses any face of an added cell
2360  forAll(addedCells, oldPatchFaceI)
2361  {
2362  // Get the cells (in newMesh labels) per old patch face (in mesh
2363  // labels)
2364  const labelList& fCells = addedCells[oldPatchFaceI];
2365 
2366  if (cellsUseFace(newMesh, fCells, wrongFaces))
2367  {
2368  // Unmark points on old mesh
2369  if
2370  (
2371  unmarkExtrusion
2372  (
2373  pp.localFaces()[oldPatchFaceI],
2374  patchDisp,
2375  patchNLayers,
2376  extrudeStatus
2377  )
2378  )
2379  {
2380  nChanged++;
2381  }
2382  }
2383  }
2384 
2385  return returnReduce(nChanged, sumOp<label>());
2386 }
2387 
2388 
2389 //- Count global number of extruded faces
2390 Foam::label Foam::autoLayerDriver::countExtrusion
2391 (
2392  const indirectPrimitivePatch& pp,
2393  const List<extrudeMode>& extrudeStatus
2394 )
2395 {
2396  // Count number of extruded patch faces
2397  label nExtruded = 0;
2398  {
2399  const faceList& localFaces = pp.localFaces();
2400 
2401  forAll(localFaces, i)
2402  {
2403  const face& localFace = localFaces[i];
2404 
2405  forAll(localFace, fp)
2406  {
2407  if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
2408  {
2409  nExtruded++;
2410  break;
2411  }
2412  }
2413  }
2414  }
2415 
2416  return returnReduce(nExtruded, sumOp<label>());
2417 }
2418 
2419 
2420 // Collect layer faces and layer cells into bools for ease of handling
2421 void Foam::autoLayerDriver::getLayerCellsFaces
2422 (
2423  const polyMesh& mesh,
2424  const addPatchCellLayer& addLayer,
2425  boolList& flaggedCells,
2426  boolList& flaggedFaces
2427 )
2428 {
2429  flaggedCells.setSize(mesh.nCells());
2430  flaggedCells = false;
2431  flaggedFaces.setSize(mesh.nFaces());
2432  flaggedFaces = false;
2433 
2434  // Mark all faces in the layer
2435  const labelListList& layerFaces = addLayer.layerFaces();
2436 
2437  // Mark all cells in the layer.
2438  labelListList addedCells(addPatchCellLayer::addedCells(mesh, layerFaces));
2439 
2440  forAll(addedCells, oldPatchFaceI)
2441  {
2442  const labelList& added = addedCells[oldPatchFaceI];
2443 
2444  forAll(added, i)
2445  {
2446  flaggedCells[added[i]] = true;
2447  }
2448  }
2449 
2450  forAll(layerFaces, oldPatchFaceI)
2451  {
2452  const labelList& layer = layerFaces[oldPatchFaceI];
2453 
2454  if (layer.size())
2455  {
2456  for (label i = 1; i < layer.size()-1; i++)
2457  {
2458  flaggedFaces[layer[i]] = true;
2459  }
2460  }
2461  }
2462 }
2463 
2464 
2465 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2466 
2467 Foam::autoLayerDriver::autoLayerDriver(meshRefinement& meshRefiner)
2468 :
2469  meshRefiner_(meshRefiner)
2470 {}
2471 
2472 
2473 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2474 
2475 void Foam::autoLayerDriver::mergePatchFacesUndo
2477  const layerParameters& layerParams,
2478  const dictionary& motionDict
2479 )
2480 {
2481  scalar minCos = Foam::cos
2482  (
2483  layerParams.featureAngle()
2484  * mathematicalConstant::pi/180.0
2485  );
2486 
2487  scalar concaveCos = Foam::cos
2488  (
2489  layerParams.concaveAngle()
2490  * mathematicalConstant::pi/180.0
2491  );
2492 
2493  Info<< nl
2494  << "Merging all faces of a cell" << nl
2495  << "---------------------------" << nl
2496  << " - which are on the same patch" << nl
2497  << " - which make an angle < " << layerParams.featureAngle()
2498  << " degrees"
2499  << nl
2500  << " (cos:" << minCos << ')' << nl
2501  << " - as long as the resulting face doesn't become concave"
2502  << " by more than "
2503  << layerParams.concaveAngle() << " degrees" << nl
2504  << " (0=straight, 180=fully concave)" << nl
2505  << endl;
2506 
2507  label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
2508 
2509  nChanged += mergeEdgesUndo(minCos, motionDict);
2510 }
2511 
2512 
2515  const layerParameters& layerParams,
2516  const dictionary& motionDict,
2517  const labelList& patchIDs,
2518  const label nAllowableErrors,
2519  decompositionMethod& decomposer,
2520  fvMeshDistribute& distributor
2521 )
2522 {
2523  fvMesh& mesh = meshRefiner_.mesh();
2524 
2526  (
2528  (
2529  mesh,
2530  patchIDs
2531  )
2532  );
2533 
2534  // Construct iterative mesh mover.
2535  Info<< "Constructing mesh displacer ..." << endl;
2536 
2537  autoPtr<motionSmoother> meshMover
2538  (
2539  new motionSmoother
2540  (
2541  mesh,
2542  pp(),
2543  patchIDs,
2545  (
2546  pointMesh::New(mesh),
2547  patchIDs
2548  ),
2549  motionDict
2550  )
2551  );
2552 
2553 
2554  // Point-wise extrusion data
2555  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2556 
2557  // Displacement for all pp.localPoints.
2558  vectorField patchDisp(pp().nPoints(), vector::one);
2559 
2560  // Number of layers for all pp.localPoints. Note: only valid if
2561  // extrudeStatus = EXTRUDE.
2562  labelList patchNLayers(pp().nPoints(), 0);
2563 
2564  // Whether to add edge for all pp.localPoints.
2565  List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
2566 
2567  {
2568  // Get number of layer per point from number of layers per patch
2569  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2570 
2571  setNumLayers
2572  (
2573  layerParams.numLayers(), // per patch the num layers
2574  meshMover().adaptPatchIDs(),// patches that are being moved
2575  pp, // indirectpatch for all faces moving
2576 
2577  patchDisp,
2578  patchNLayers,
2579  extrudeStatus
2580  );
2581 
2582  // Precalculate mesh edge labels for patch edges
2583  labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
2584 
2585  // Disable extrusion on non-manifold points
2586  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2587 
2588  handleNonManifolds
2589  (
2590  pp,
2591  meshEdges,
2592 
2593  patchDisp,
2594  patchNLayers,
2595  extrudeStatus
2596  );
2597 
2598  // Disable extrusion on feature angles
2599  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2600 
2601  handleFeatureAngle
2602  (
2603  pp,
2604  meshEdges,
2605  layerParams.featureAngle()*mathematicalConstant::pi/180.0,
2606 
2607  patchDisp,
2608  patchNLayers,
2609  extrudeStatus
2610  );
2611 
2612  // Disable extrusion on warped faces
2613  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2614 
2615  // Undistorted edge length
2616  const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
2617  const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
2618 
2619  handleWarpedFaces
2620  (
2621  pp,
2622  layerParams.maxFaceThicknessRatio(),
2623  edge0Len,
2624  cellLevel,
2625 
2626  patchDisp,
2627  patchNLayers,
2628  extrudeStatus
2629  );
2630 
2633  //
2634  //handleMultiplePatchFaces
2635  //(
2636  // pp,
2637  //
2638  // patchDisp,
2639  // patchNLayers,
2640  // extrudeStatus
2641  //);
2642 
2643 
2644  // Grow out region of non-extrusion
2645  for (label i = 0; i < layerParams.nGrow(); i++)
2646  {
2647  growNoExtrusion
2648  (
2649  pp,
2650  patchDisp,
2651  patchNLayers,
2652  extrudeStatus
2653  );
2654  }
2655  }
2656 
2657 
2658  // Undistorted edge length
2659  const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
2660  const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
2661 
2662  // Determine (wanted) point-wise layer thickness and expansion ratio
2663  scalarField thickness(pp().nPoints());
2664  scalarField minThickness(pp().nPoints());
2665  scalarField expansionRatio(pp().nPoints());
2666  calculateLayerThickness
2667  (
2668  pp,
2669  meshMover().adaptPatchIDs(),
2670  layerParams.expansionRatio(),
2671 
2672  layerParams.relativeSizes(), // thickness relative to cellsize?
2673  layerParams.finalLayerThickness(), // wanted thicknes
2674  layerParams.minThickness(), // minimum thickness
2675 
2676  cellLevel,
2677  patchNLayers,
2678  edge0Len,
2679 
2680  thickness,
2681  minThickness,
2682  expansionRatio
2683  );
2684 
2685 
2686  // Print a bit
2687  {
2688  const polyBoundaryMesh& patches = mesh.boundaryMesh();
2689 
2690  // Find maximum length of a patch name, for a nicer output
2691  label maxPatchNameLen = 0;
2692  forAll(meshMover().adaptPatchIDs(), i)
2693  {
2694  label patchI = meshMover().adaptPatchIDs()[i];
2695  word patchName = patches[patchI].name();
2696  maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
2697  }
2698 
2699  Info<< nl
2700  << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
2701  << setw(0) << " faces layers avg thickness[m]" << nl
2702  << setf(ios_base::left) << setw(maxPatchNameLen) << " "
2703  << setw(0) << " near-wall overall" << nl
2704  << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
2705  << setw(0) << " ----- ------ --------- -------" << endl;
2706 
2707  forAll(meshMover().adaptPatchIDs(), i)
2708  {
2709  label patchI = meshMover().adaptPatchIDs()[i];
2710 
2711  const labelList& meshPoints = patches[patchI].meshPoints();
2712 
2713  //scalar maxThickness = -VGREAT;
2714  //scalar minThickness = VGREAT;
2715  scalar sumThickness = 0;
2716  scalar sumNearWallThickness = 0;
2717 
2718  forAll(meshPoints, patchPointI)
2719  {
2720  label ppPointI = pp().meshPointMap()[meshPoints[patchPointI]];
2721 
2722  //maxThickness = max(maxThickness, thickness[ppPointI]);
2723  //minThickness = min(minThickness, thickness[ppPointI]);
2724  sumThickness += thickness[ppPointI];
2725 
2726  label nLay = patchNLayers[ppPointI];
2727  if (nLay > 0)
2728  {
2729  if (expansionRatio[ppPointI] == 1)
2730  {
2731  sumNearWallThickness += thickness[ppPointI]/nLay;
2732  }
2733  else
2734  {
2735  scalar s =
2736  (1.0-pow(expansionRatio[ppPointI], nLay))
2737  / (1.0-expansionRatio[ppPointI]);
2738  sumNearWallThickness += thickness[ppPointI]/s;
2739  }
2740  }
2741  }
2742 
2743  label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
2744 
2745  // For empty patches, totNPoints is 0.
2746  scalar avgThickness = 0;
2747  scalar avgNearWallThickness = 0;
2748 
2749  if (totNPoints > 0)
2750  {
2751  //reduce(maxThickness, maxOp<scalar>());
2752  //reduce(minThickness, minOp<scalar>());
2753  avgThickness =
2754  returnReduce(sumThickness, sumOp<scalar>())
2755  / totNPoints;
2756  avgNearWallThickness =
2757  returnReduce(sumNearWallThickness, sumOp<scalar>())
2758  / totNPoints;
2759  }
2760 
2761  Info<< setf(ios_base::left) << setw(maxPatchNameLen)
2762  << patches[patchI].name() << setprecision(3)
2763  << " " << setw(8)
2764  << returnReduce(patches[patchI].size(), sumOp<scalar>())
2765  << " " << setw(6) << layerParams.numLayers()[patchI]
2766  << " " << setw(8) << avgNearWallThickness
2767  << " " << setw(8) << avgThickness
2768  //<< " " << setw(8) << minThickness
2769  //<< " " << setw(8) << maxThickness
2770  << endl;
2771  }
2772  Info<< endl;
2773  }
2774 
2775 
2776  // Calculate wall to medial axis distance for smoothing displacement
2777  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2778 
2779  pointScalarField pointMedialDist
2780  (
2781  IOobject
2782  (
2783  "pointMedialDist",
2784  meshRefiner_.timeName(),
2785  mesh,
2788  false
2789  ),
2790  meshMover().pMesh(),
2791  dimensionedScalar("pointMedialDist", dimless, 0.0)
2792  );
2793 
2794  pointVectorField dispVec
2795  (
2796  IOobject
2797  (
2798  "dispVec",
2799  meshRefiner_.timeName(),
2800  mesh,
2803  false
2804  ),
2805  meshMover().pMesh(),
2807  );
2808 
2809  pointScalarField medialRatio
2810  (
2811  IOobject
2812  (
2813  "medialRatio",
2814  meshRefiner_.timeName(),
2815  mesh,
2818  false
2819  ),
2820  meshMover().pMesh(),
2821  dimensionedScalar("medialRatio", dimless, 0.0)
2822  );
2823 
2824  // Setup information for medial axis smoothing. Calculates medial axis
2825  // and a smoothed displacement direction.
2826  // - pointMedialDist : distance to medial axis
2827  // - dispVec : normalised direction of nearest displacement
2828  // - medialRatio : ratio of medial distance to wall distance.
2829  // (1 at wall, 0 at medial axis)
2830  medialAxisSmoothingInfo
2831  (
2832  meshMover,
2833  layerParams.nSmoothNormals(),
2834  layerParams.nSmoothSurfaceNormals(),
2835  layerParams.minMedianAxisAngleCos(),
2836 
2837  dispVec,
2838  medialRatio,
2839  pointMedialDist
2840  );
2841 
2842 
2843 
2844  // Saved old points
2845  pointField oldPoints(mesh.points());
2846 
2847  // Last set of topology changes. (changing mesh clears out polyTopoChange)
2848  polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
2849 
2850  boolList flaggedCells;
2851  boolList flaggedFaces;
2852 
2853  for (label iteration = 0; iteration < layerParams.nLayerIter(); iteration++)
2854  {
2855  Info<< nl
2856  << "Layer addition iteration " << iteration << nl
2857  << "--------------------------" << endl;
2858 
2859 
2860  // Unset the extrusion at the pp.
2861  const dictionary& meshQualityDict =
2862  (
2863  iteration < layerParams.nRelaxedIter()
2864  ? motionDict
2865  : motionDict.subDict("relaxed")
2866  );
2867 
2868  if (iteration >= layerParams.nRelaxedIter())
2869  {
2870  Info<< "Switched to relaxed meshQuality constraints." << endl;
2871  }
2872 
2873 
2874 
2875  // Make sure displacement is equal on both sides of coupled patches.
2876  syncPatchDisplacement
2877  (
2878  meshMover,
2879  minThickness,
2880  patchDisp,
2881  patchNLayers,
2882  extrudeStatus
2883  );
2884 
2885  // Displacement acc. to pointnormals
2886  getPatchDisplacement
2887  (
2888  meshMover,
2889  thickness,
2890  minThickness,
2891  patchDisp,
2892  patchNLayers,
2893  extrudeStatus
2894  );
2895 
2896  // Shrink mesh by displacement value first.
2897  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2898 
2899  {
2900  pointField oldPatchPos(pp().localPoints());
2901 
2903  //shrinkMeshDistance
2904  //(
2905  // debug,
2906  // meshMover,
2907  // -patchDisp, // Shrink in opposite direction of addedPoints
2908  // layerParams.nSmoothDisp(),
2909  // layerParams.nSnap()
2910  //);
2911 
2912  // Medial axis based shrinking
2913  shrinkMeshMedialDistance
2914  (
2915  meshMover(),
2916  meshQualityDict,
2917 
2918  layerParams.nSmoothThickness(),
2919  layerParams.maxThicknessToMedialRatio(),
2920  nAllowableErrors,
2921  layerParams.nSnap(),
2922  layerParams.layerTerminationCos(),
2923 
2924  thickness,
2925  minThickness,
2926 
2927  dispVec,
2928  medialRatio,
2929  pointMedialDist,
2930 
2931  extrudeStatus,
2932  patchDisp,
2933  patchNLayers
2934  );
2935 
2936  // Update patchDisp (since not all might have been honoured)
2937  patchDisp = oldPatchPos - pp().localPoints();
2938  }
2939 
2940  // Truncate displacements that are too small (this will do internal
2941  // ones, coupled ones have already been truncated by syncPatch)
2942  faceSet dummySet(mesh, "wrongPatchFaces", 0);
2943  truncateDisplacement
2944  (
2945  meshMover(),
2946  minThickness,
2947  dummySet,
2948  patchDisp,
2949  patchNLayers,
2950  extrudeStatus
2951  );
2952 
2953 
2954  // Dump to .obj file for debugging.
2955  if (debug)
2956  {
2957  dumpDisplacement
2958  (
2959  mesh.time().path()/"layer",
2960  pp(),
2961  patchDisp,
2962  extrudeStatus
2963  );
2964 
2965  const_cast<Time&>(mesh.time())++;
2966  Info<< "Writing shrunk mesh to " << meshRefiner_.timeName() << endl;
2967 
2968  // See comment in autoSnapDriver why we should not remove meshPhi
2969  // using mesh.clearPout().
2970 
2971  mesh.write();
2972  }
2973 
2974 
2975  // Mesh topo change engine
2976  polyTopoChange meshMod(mesh);
2977 
2978  // Grow layer of cells on to patch. Handles zero sized displacement.
2979  addPatchCellLayer addLayer(mesh);
2980 
2981  // Determine per point/per face number of layers to extrude. Also
2982  // handles the slow termination of layers when going switching layers
2983 
2984  labelList nPatchPointLayers(pp().nPoints(),-1);
2985  labelList nPatchFaceLayers(pp().localFaces().size(),-1);
2986  setupLayerInfoTruncation
2987  (
2988  meshMover(),
2989  patchNLayers,
2990  extrudeStatus,
2991  layerParams.nBufferCellsNoExtrude(),
2992  nPatchPointLayers,
2993  nPatchFaceLayers
2994  );
2995 
2996  // Calculate displacement for first layer for addPatchLayer.
2997  // (first layer = layer of cells next to the original mesh)
2998  vectorField firstDisp(patchNLayers.size(), vector::zero);
2999 
3000  forAll(patchNLayers, i)
3001  {
3002  if (patchNLayers[i] > 0)
3003  {
3004  if (expansionRatio[i] == 1.0)
3005  {
3006  firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
3007  }
3008  else
3009  {
3010  label nLay = nPatchPointLayers[i];
3011  scalar h =
3012  pow(expansionRatio[i], nLay - 1)
3013  * (1.0 - expansionRatio[i])
3014  / (1.0 - pow(expansionRatio[i], nLay));
3015  firstDisp[i] = h*patchDisp[i];
3016  }
3017  }
3018  }
3019 
3020  scalarField invExpansionRatio = 1.0 / expansionRatio;
3021 
3022  // Add topo regardless of whether extrudeStatus is extruderemove.
3023  // Not add layer if patchDisp is zero.
3024  addLayer.setRefinement
3025  (
3026  invExpansionRatio,
3027  pp(),
3028  nPatchFaceLayers, // layers per face
3029  nPatchPointLayers, // layers per point
3030  firstDisp, // thickness of layer nearest internal mesh
3031  meshMod
3032  );
3033 
3034  if (debug)
3035  {
3036  const_cast<Time&>(mesh.time())++;
3037  }
3038 
3039  // Store mesh changes for if mesh is correct.
3040  savedMeshMod = meshMod;
3041 
3042 
3043  // With the stored topo changes we create a new mesh so we can
3044  // undo if necessary.
3045 
3046  autoPtr<fvMesh> newMeshPtr;
3047  autoPtr<mapPolyMesh> map = meshMod.makeMesh
3048  (
3049  newMeshPtr,
3050  IOobject
3051  (
3052  //mesh.name()+"_layer",
3053  mesh.name(),
3054  static_cast<polyMesh&>(mesh).instance(),
3055  mesh.time(), // register with runTime
3056  static_cast<polyMesh&>(mesh).readOpt(),
3057  static_cast<polyMesh&>(mesh).writeOpt()
3058  ), // io params from original mesh but new name
3059  mesh, // original mesh
3060  true // parallel sync
3061  );
3062  fvMesh& newMesh = newMeshPtr();
3063 
3064  //?neccesary? Update fields
3065  newMesh.updateMesh(map);
3066 
3067  if (meshRefiner_.overwrite())
3068  {
3069  newMesh.setInstance(meshRefiner_.oldInstance());
3070  }
3071 
3072  // Update numbering on addLayer:
3073  // - cell/point labels to be newMesh.
3074  // - patchFaces to remain in oldMesh order.
3075  addLayer.updateMesh
3076  (
3077  map,
3078  identity(pp().size()),
3079  identity(pp().nPoints())
3080  );
3081 
3082  // Collect layer faces and cells for outside loop.
3083  getLayerCellsFaces
3084  (
3085  newMesh,
3086  addLayer,
3087  flaggedCells,
3088  flaggedFaces
3089  );
3090 
3091 
3092  if (debug)
3093  {
3094  Info<< "Writing layer mesh to " << meshRefiner_.timeName() << endl;
3095  newMesh.write();
3096  cellSet addedCellSet
3097  (
3098  newMesh,
3099  "addedCells",
3100  findIndices(flaggedCells, true)
3101  );
3102  Info<< "Writing "
3103  << returnReduce(addedCellSet.size(), sumOp<label>())
3104  << " added cells to cellSet "
3105  << addedCellSet.name() << endl;
3106  addedCellSet.write();
3107 
3108  faceSet layerFacesSet
3109  (
3110  newMesh,
3111  "layerFaces",
3112  findIndices(flaggedCells, true)
3113  );
3114  Info<< "Writing "
3115  << returnReduce(layerFacesSet.size(), sumOp<label>())
3116  << " faces inside added layer to faceSet "
3117  << layerFacesSet.name() << endl;
3118  layerFacesSet.write();
3119  }
3120 
3121 
3122  label nTotChanged = checkAndUnmark
3123  (
3124  addLayer,
3125  meshQualityDict,
3126  pp(),
3127  newMesh,
3128 
3129  patchDisp,
3130  patchNLayers,
3131  extrudeStatus
3132  );
3133 
3134  Info<< "Extruding " << countExtrusion(pp, extrudeStatus)
3135  << " out of " << returnReduce(pp().size(), sumOp<label>())
3136  << " faces. Removed extrusion at " << nTotChanged << " faces."
3137  << endl;
3138 
3139  if (nTotChanged == 0)
3140  {
3141  break;
3142  }
3143 
3144  // Reset mesh points and start again
3145  meshMover().movePoints(oldPoints);
3146  meshMover().correct();
3147 
3148  Info<< endl;
3149  }
3150 
3151 
3152  // At this point we have a (shrunk) mesh and a set of topology changes
3153  // which will make a valid mesh with layer. Apply these changes to the
3154  // current mesh.
3155 
3156  // Apply the stored topo changes to the current mesh.
3157  autoPtr<mapPolyMesh> map = savedMeshMod.changeMesh(mesh, false);
3158 
3159  // Update fields
3160  mesh.updateMesh(map);
3161 
3162  // Move mesh (since morphing does not do this)
3163  if (map().hasMotionPoints())
3164  {
3165  mesh.movePoints(map().preMotionPoints());
3166  }
3167  else
3168  {
3169  // Delete mesh volumes.
3170  mesh.clearOut();
3171  }
3172 
3173  if (meshRefiner_.overwrite())
3174  {
3175  mesh.setInstance(meshRefiner_.oldInstance());
3176  }
3177 
3178  meshRefiner_.updateMesh(map, labelList(0));
3179 
3180 
3181  // Do final balancing
3182  // ~~~~~~~~~~~~~~~~~~
3183 
3184  if (Pstream::parRun())
3185  {
3186  Info<< nl
3187  << "Doing final balancing" << nl
3188  << "---------------------" << nl
3189  << endl;
3190 
3191  if (debug)
3192  {
3193  const_cast<Time&>(mesh.time())++;
3194  }
3195 
3196  // Balance. No restriction on face zones and baffles.
3197  autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
3198  (
3199  false,
3200  false,
3201  scalarField(mesh.nCells(), 1.0),
3202  decomposer,
3203  distributor
3204  );
3205 
3206  // Re-distribute flag of layer faces and cells
3207  map().distributeCellData(flaggedCells);
3208  map().distributeFaceData(flaggedFaces);
3209  }
3210 
3211 
3212  // Write mesh
3213  // ~~~~~~~~~~
3214 
3215  cellSet addedCellSet(mesh, "addedCells", findIndices(flaggedCells, true));
3216  Info<< "Writing "
3217  << returnReduce(addedCellSet.size(), sumOp<label>())
3218  << " added cells to cellSet "
3219  << addedCellSet.name() << endl;
3220  addedCellSet.write();
3221 
3222  faceSet layerFacesSet(mesh, "layerFaces", findIndices(flaggedFaces, true));
3223  Info<< "Writing "
3224  << returnReduce(layerFacesSet.size(), sumOp<label>())
3225  << " faces inside added layer to faceSet "
3226  << layerFacesSet.name() << endl;
3227  layerFacesSet.write();
3228 }
3229 
3230 
3233  const dictionary& shrinkDict,
3234  const dictionary& motionDict,
3235  const layerParameters& layerParams,
3236  const bool preBalance,
3237  decompositionMethod& decomposer,
3238  fvMeshDistribute& distributor
3239 )
3240 {
3241  const fvMesh& mesh = meshRefiner_.mesh();
3242 
3243  Info<< nl
3244  << "Shrinking and layer addition phase" << nl
3245  << "----------------------------------" << nl
3246  << endl;
3247 
3248  Info<< "Using mesh parameters " << motionDict << nl << endl;
3249 
3250  // Merge coplanar boundary faces
3251  mergePatchFacesUndo(layerParams, motionDict);
3252 
3253  // Per patch the number of layers (0 if no layer)
3254  const labelList& numLayers = layerParams.numLayers();
3255 
3256  // Patches that need to get a layer
3257  DynamicList<label> patchIDs(numLayers.size());
3258  label nFacesWithLayers = 0;
3259  forAll(numLayers, patchI)
3260  {
3261  if (numLayers[patchI] > 0)
3262  {
3263  patchIDs.append(patchI);
3264  nFacesWithLayers += mesh.boundaryMesh()[patchI].size();
3265  }
3266  }
3267  patchIDs.shrink();
3268 
3269  if (returnReduce(nFacesWithLayers, sumOp<label>()) == 0)
3270  {
3271  Info<< nl << "No layers to generate ..." << endl;
3272  }
3273  else
3274  {
3275  // Check that outside of mesh is not multiply connected.
3276  checkMeshManifold();
3277 
3278  // Check initial mesh
3279  Info<< "Checking initial mesh ..." << endl;
3280  labelHashSet wrongFaces(mesh.nFaces()/100);
3281  motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
3282  const label nInitErrors = returnReduce
3283  (
3284  wrongFaces.size(),
3285  sumOp<label>()
3286  );
3287 
3288  Info<< "Detected " << nInitErrors << " illegal faces"
3289  << " (concave, zero area or negative cell pyramid volume)"
3290  << endl;
3291 
3292 
3293  // Balance
3294  if (Pstream::parRun() && preBalance)
3295  {
3296  Info<< nl
3297  << "Doing initial balancing" << nl
3298  << "-----------------------" << nl
3299  << endl;
3300 
3301  scalarField cellWeights(mesh.nCells(), 1);
3302  forAll(numLayers, patchI)
3303  {
3304  if (numLayers[patchI] > 0)
3305  {
3306  const polyPatch& pp = mesh.boundaryMesh()[patchI];
3307  forAll(pp.faceCells(), i)
3308  {
3309  cellWeights[pp.faceCells()[i]] += numLayers[patchI];
3310  }
3311  }
3312  }
3313 
3314  // Balance mesh (and meshRefinement). No restriction on face zones
3315  // and baffles.
3316  autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
3317  (
3318  false,
3319  false,
3320  cellWeights,
3321  decomposer,
3322  distributor
3323  );
3324 
3325  //{
3326  // globalIndex globalCells(mesh.nCells());
3327  //
3328  // Info<< "** Distribution after balancing:" << endl;
3329  // for (label procI = 0; procI < Pstream::nProcs(); procI++)
3330  // {
3331  // Info<< " " << procI << '\t'
3332  // << globalCells.localSize(procI) << endl;
3333  // }
3334  // Info<< endl;
3335  //}
3336  }
3337 
3338 
3339  // Do all topo changes
3340  addLayers
3341  (
3342  layerParams,
3343  motionDict,
3344  patchIDs,
3345  nInitErrors,
3346  decomposer,
3347  distributor
3348  );
3349  }
3350 }
3351 
3352 
3353 // ************************ vim: set sw=4 sts=4 et: ************************ //