snappyLayerDriver.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) 2011-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 Description
25  All to do with adding cell layers
26 
27 \*----------------------------------------------------------------------------*/
28 
29 #include "snappyLayerDriver.H"
30 #include "fvMesh.H"
31 #include "Time.H"
32 #include "meshRefinement.H"
33 #include "removePoints.H"
34 #include "pointFields.H"
35 #include "motionSmoother.H"
36 #include "unitConversion.H"
37 #include "pointSet.H"
38 #include "faceSet.H"
39 #include "cellSet.H"
40 #include "polyTopoChange.H"
41 #include "mapPolyMesh.H"
42 #include "addPatchCellLayer.H"
43 #include "mapDistributePolyMesh.H"
44 #include "OBJstream.H"
45 #include "layerParameters.H"
46 #include "combineFaces.H"
47 #include "IOmanip.H"
48 #include "globalIndex.H"
49 #include "DynamicField.H"
50 #include "PatchTools.H"
51 #include "slipPointPatchFields.H"
57 #include "localPointRegion.H"
58 
60 #include "medialAxisMeshMover.H"
61 #include "scalarIOField.H"
62 
63 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
64 
65 namespace Foam
66 {
67 
68 defineTypeNameAndDebug(snappyLayerDriver, 0);
69 
70 } // End namespace Foam
71 
72 
73 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
74 
75 // For debugging: Dump displacement to .obj files
76 void Foam::snappyLayerDriver::dumpDisplacement
77 (
78  const fileName& prefix,
79  const indirectPrimitivePatch& pp,
80  const vectorField& patchDisp,
81  const List<extrudeMode>& extrudeStatus
82 )
83 {
84  OBJstream dispStr(prefix + "_disp.obj");
85  Info<< "Writing all displacements to " << dispStr.name() << endl;
86 
87  forAll(patchDisp, patchPointi)
88  {
89  const point& pt = pp.localPoints()[patchPointi];
90  dispStr.write(linePointRef(pt, pt + patchDisp[patchPointi]));
91  }
92 
93 
94  OBJstream illStr(prefix + "_illegal.obj");
95  Info<< "Writing invalid displacements to " << illStr.name() << endl;
96 
97  forAll(patchDisp, patchPointi)
98  {
99  if (extrudeStatus[patchPointi] != EXTRUDE)
100  {
101  const point& pt = pp.localPoints()[patchPointi];
102  illStr.write(linePointRef(pt, pt + patchDisp[patchPointi]));
103  }
104  }
105 }
106 
107 
108 Foam::tmp<Foam::scalarField> Foam::snappyLayerDriver::avgPointData
109 (
110  const indirectPrimitivePatch& pp,
111  const scalarField& pointFld
112 )
113 {
114  tmp<scalarField> tfaceFld(new scalarField(pp.size(), 0.0));
115  scalarField& faceFld = tfaceFld.ref();
116 
117  forAll(pp.localFaces(), facei)
118  {
119  const face& f = pp.localFaces()[facei];
120  if (f.size())
121  {
122  forAll(f, fp)
123  {
124  faceFld[facei] += pointFld[f[fp]];
125  }
126  faceFld[facei] /= f.size();
127  }
128  }
129  return tfaceFld;
130 }
131 
132 
133 // Check that primitivePatch is not multiply connected. Collect non-manifold
134 // points in pointSet.
135 void Foam::snappyLayerDriver::checkManifold
136 (
137  const indirectPrimitivePatch& fp,
138  pointSet& nonManifoldPoints
139 )
140 {
141  // Check for non-manifold points (surface pinched at point)
142  fp.checkPointManifold(false, &nonManifoldPoints);
143 
144  // Check for edge-faces (surface pinched at edge)
145  const labelListList& edgeFaces = fp.edgeFaces();
146 
147  forAll(edgeFaces, edgeI)
148  {
149  const labelList& eFaces = edgeFaces[edgeI];
150 
151  if (eFaces.size() > 2)
152  {
153  const edge& e = fp.edges()[edgeI];
154 
155  nonManifoldPoints.insert(fp.meshPoints()[e[0]]);
156  nonManifoldPoints.insert(fp.meshPoints()[e[1]]);
157  }
158  }
159 }
160 
161 
162 void Foam::snappyLayerDriver::checkMeshManifold() const
163 {
164  const fvMesh& mesh = meshRefiner_.mesh();
165 
166  Info<< nl << "Checking mesh manifoldness ..." << endl;
167 
168  // Get all outside faces
169  labelList outsideFaces(mesh.nFaces() - mesh.nInternalFaces());
170 
171  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
172  {
173  outsideFaces[facei - mesh.nInternalFaces()] = facei;
174  }
175 
176  pointSet nonManifoldPoints
177  (
178  mesh,
179  "nonManifoldPoints",
180  mesh.nPoints() / 100
181  );
182 
183  // Build primitivePatch out of faces and check it for problems.
184  checkManifold
185  (
187  (
188  IndirectList<face>(mesh.faces(), outsideFaces),
189  mesh.points()
190  ),
191  nonManifoldPoints
192  );
193 
194  label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
195 
196  if (nNonManif > 0)
197  {
198  Info<< "Outside of mesh is multiply connected across edges or"
199  << " points." << nl
200  << "This is not a fatal error but might cause some unexpected"
201  << " behaviour." << nl
202  //<< "Writing " << nNonManif
203  //<< " points where this happens to pointSet "
204  //<< nonManifoldPoints.name()
205  << endl;
206 
207  //nonManifoldPoints.instance() = meshRefiner_.timeName();
208  //nonManifoldPoints.write();
209  }
210  Info<< endl;
211 }
212 
213 
214 
215 // Unset extrusion on point. Returns true if anything unset.
216 bool Foam::snappyLayerDriver::unmarkExtrusion
217 (
218  const label patchPointi,
219  pointField& patchDisp,
220  labelList& patchNLayers,
221  List<extrudeMode>& extrudeStatus
222 )
223 {
224  if (extrudeStatus[patchPointi] == EXTRUDE)
225  {
226  extrudeStatus[patchPointi] = NOEXTRUDE;
227  patchNLayers[patchPointi] = 0;
228  patchDisp[patchPointi] = Zero;
229  return true;
230  }
231  else if (extrudeStatus[patchPointi] == EXTRUDEREMOVE)
232  {
233  extrudeStatus[patchPointi] = NOEXTRUDE;
234  patchNLayers[patchPointi] = 0;
235  patchDisp[patchPointi] = Zero;
236  return true;
237  }
238  else
239  {
240  return false;
241  }
242 }
243 
244 
245 // Unset extrusion on face. Returns true if anything unset.
246 bool Foam::snappyLayerDriver::unmarkExtrusion
247 (
248  const face& localFace,
249  pointField& patchDisp,
250  labelList& patchNLayers,
251  List<extrudeMode>& extrudeStatus
252 )
253 {
254  bool unextruded = false;
255 
256  forAll(localFace, fp)
257  {
258  if
259  (
260  unmarkExtrusion
261  (
262  localFace[fp],
263  patchDisp,
264  patchNLayers,
265  extrudeStatus
266  )
267  )
268  {
269  unextruded = true;
270  }
271  }
272  return unextruded;
273 }
274 
275 
276 // No extrusion at non-manifold points.
277 void Foam::snappyLayerDriver::handleNonManifolds
278 (
279  const indirectPrimitivePatch& pp,
280  const labelList& meshEdges,
281  const labelListList& edgeGlobalFaces,
282  pointField& patchDisp,
283  labelList& patchNLayers,
284  List<extrudeMode>& extrudeStatus
285 ) const
286 {
287  const fvMesh& mesh = meshRefiner_.mesh();
288 
289  Info<< nl << "Handling non-manifold points ..." << endl;
290 
291  // Detect non-manifold points
292  Info<< nl << "Checking patch manifoldness ..." << endl;
293 
294  pointSet nonManifoldPoints(mesh, "nonManifoldPoints", pp.nPoints());
295 
296  // 1. Local check
297  checkManifold(pp, nonManifoldPoints);
298 
299  // 2. Remote check for boundary edges on coupled boundaries
300  forAll(edgeGlobalFaces, edgeI)
301  {
302  if
303  (
304  pp.edgeFaces()[edgeI].size() == 1
305  && edgeGlobalFaces[edgeI].size() > 2
306  )
307  {
308  // So boundary edges that are connected to more than 2 processors
309  // i.e. a non-manifold edge which is exactly on a processor
310  // boundary.
311  const edge& e = pp.edges()[edgeI];
312  nonManifoldPoints.insert(pp.meshPoints()[e[0]]);
313  nonManifoldPoints.insert(pp.meshPoints()[e[1]]);
314  }
315  }
316 
317  // 3. Remote check for end of layer across coupled boundaries
318  {
319  PackedBoolList isCoupledEdge(mesh.nEdges());
320 
321  const labelList& cpEdges = mesh.globalData().coupledPatchMeshEdges();
322  forAll(cpEdges, i)
323  {
324  isCoupledEdge[cpEdges[i]] = true;
325  }
327  (
328  mesh,
329  isCoupledEdge,
330  orEqOp<unsigned int>(),
331  0
332  );
333 
334  forAll(edgeGlobalFaces, edgeI)
335  {
336  label meshEdgeI = meshEdges[edgeI];
337 
338  if
339  (
340  pp.edgeFaces()[edgeI].size() == 1
341  && edgeGlobalFaces[edgeI].size() == 1
342  && isCoupledEdge[meshEdgeI]
343  )
344  {
345  // Edge of patch but no continuation across processor.
346  const edge& e = pp.edges()[edgeI];
347  //Pout<< "** Stopping extrusion on edge "
348  // << pp.localPoints()[e[0]]
349  // << pp.localPoints()[e[1]] << endl;
350  nonManifoldPoints.insert(pp.meshPoints()[e[0]]);
351  nonManifoldPoints.insert(pp.meshPoints()[e[1]]);
352  }
353  }
354  }
355 
356 
357 
358  label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
359 
360  Info<< "Outside of local patch is multiply connected across edges or"
361  << " points at " << nNonManif << " points." << endl;
362 
363  if (nNonManif > 0)
364  {
365  const labelList& meshPoints = pp.meshPoints();
366 
367  forAll(meshPoints, patchPointi)
368  {
369  if (nonManifoldPoints.found(meshPoints[patchPointi]))
370  {
371  unmarkExtrusion
372  (
373  patchPointi,
374  patchDisp,
375  patchNLayers,
376  extrudeStatus
377  );
378  }
379  }
380  }
381 
382  Info<< "Set displacement to zero for all " << nNonManif
383  << " non-manifold points" << endl;
384 }
385 
386 
387 // Parallel feature edge detection. Assumes non-manifold edges already handled.
388 void Foam::snappyLayerDriver::handleFeatureAngle
389 (
390  const indirectPrimitivePatch& pp,
391  const labelList& meshEdges,
392  const scalar minCos,
393  pointField& patchDisp,
394  labelList& patchNLayers,
395  List<extrudeMode>& extrudeStatus
396 ) const
397 {
398  const fvMesh& mesh = meshRefiner_.mesh();
399 
400  Info<< nl << "Handling feature edges ..." << endl;
401 
402  if (minCos < 1-SMALL)
403  {
404  // Normal component of normals of connected faces.
405  vectorField edgeNormal(mesh.nEdges(), point::max);
406 
407  const labelListList& edgeFaces = pp.edgeFaces();
408 
409  forAll(edgeFaces, edgeI)
410  {
411  const labelList& eFaces = pp.edgeFaces()[edgeI];
412 
413  label meshEdgeI = meshEdges[edgeI];
414 
415  forAll(eFaces, i)
416  {
417  nomalsCombine()
418  (
419  edgeNormal[meshEdgeI],
420  pp.faceNormals()[eFaces[i]]
421  );
422  }
423  }
424 
426  (
427  mesh,
428  edgeNormal,
429  nomalsCombine(),
430  point::max // null value
431  );
432 
433  autoPtr<OBJstream> str;
434  if (debug&meshRefinement::MESH)
435  {
436  str.reset
437  (
438  new OBJstream
439  (
440  mesh.time().path()
441  / "featureEdges_"
442  + meshRefiner_.timeName()
443  + ".obj"
444  )
445  );
446  Info<< "Writing feature edges to " << str().name() << endl;
447  }
448 
449  label nFeats = 0;
450 
451  // Now on coupled edges the edgeNormal will have been truncated and
452  // only be still be the old value where two faces have the same normal
453  forAll(edgeFaces, edgeI)
454  {
455  const labelList& eFaces = pp.edgeFaces()[edgeI];
456 
457  label meshEdgeI = meshEdges[edgeI];
458 
459  const vector& n = edgeNormal[meshEdgeI];
460 
461  if (n != point::max)
462  {
463  scalar cos = n & pp.faceNormals()[eFaces[0]];
464 
465  if (cos < minCos)
466  {
467  const edge& e = pp.edges()[edgeI];
468 
469  unmarkExtrusion
470  (
471  e[0],
472  patchDisp,
473  patchNLayers,
474  extrudeStatus
475  );
476  unmarkExtrusion
477  (
478  e[1],
479  patchDisp,
480  patchNLayers,
481  extrudeStatus
482  );
483 
484  nFeats++;
485 
486  if (str.valid())
487  {
488  const point& p0 = pp.localPoints()[e[0]];
489  const point& p1 = pp.localPoints()[e[1]];
490  str().write(linePointRef(p0, p1));
491  }
492  }
493  }
494  }
495 
496  Info<< "Set displacement to zero for points on "
497  << returnReduce(nFeats, sumOp<label>())
498  << " feature edges" << endl;
499  }
500 }
501 
502 
503 // No extrusion on cells with warped faces. Calculates the thickness of the
504 // layer and compares it to the space the warped face takes up. Disables
505 // extrusion if layer thickness is more than faceRatio of the thickness of
506 // the face.
507 void Foam::snappyLayerDriver::handleWarpedFaces
508 (
509  const indirectPrimitivePatch& pp,
510  const scalar faceRatio,
511  const scalar edge0Len,
512  const labelList& cellLevel,
513  pointField& patchDisp,
514  labelList& patchNLayers,
515  List<extrudeMode>& extrudeStatus
516 ) const
517 {
518  const fvMesh& mesh = meshRefiner_.mesh();
519 
520  Info<< nl << "Handling cells with warped patch faces ..." << nl;
521 
522  const pointField& points = mesh.points();
523 
524  label nWarpedFaces = 0;
525 
526  forAll(pp, i)
527  {
528  const face& f = pp[i];
529 
530  if (f.size() > 3)
531  {
532  label facei = pp.addressing()[i];
533 
534  label ownLevel = cellLevel[mesh.faceOwner()[facei]];
535  scalar edgeLen = edge0Len/(1<<ownLevel);
536 
537  // Normal distance to face centre plane
538  const point& fc = mesh.faceCentres()[facei];
539  const vector& fn = pp.faceNormals()[i];
540 
541  scalarField vProj(f.size());
542 
543  forAll(f, fp)
544  {
545  vector n = points[f[fp]] - fc;
546  vProj[fp] = (n & fn);
547  }
548 
549  // Get normal 'span' of face
550  scalar minVal = min(vProj);
551  scalar maxVal = max(vProj);
552 
553  if ((maxVal - minVal) > faceRatio * edgeLen)
554  {
555  if
556  (
557  unmarkExtrusion
558  (
559  pp.localFaces()[i],
560  patchDisp,
561  patchNLayers,
562  extrudeStatus
563  )
564  )
565  {
566  nWarpedFaces++;
567  }
568  }
569  }
570  }
571 
572  Info<< "Set displacement to zero on "
573  << returnReduce(nWarpedFaces, sumOp<label>())
574  << " warped faces since layer would be > " << faceRatio
575  << " of the size of the bounding box." << endl;
576 }
577 
578 
581 //void Foam::snappyLayerDriver::handleMultiplePatchFaces
582 //(
583 // const indirectPrimitivePatch& pp,
584 // pointField& patchDisp,
585 // labelList& patchNLayers,
586 // List<extrudeMode>& extrudeStatus
587 //) const
588 //{
589 // const fvMesh& mesh = meshRefiner_.mesh();
590 //
591 // Info<< nl << "Handling cells with multiple patch faces ..." << nl;
592 //
593 // const labelListList& pointFaces = pp.pointFaces();
594 //
595 // // Cells that should not get an extrusion layer
596 // cellSet multiPatchCells(mesh, "multiPatchCells", pp.size());
597 //
598 // // Detect points that use multiple faces on same cell.
599 // forAll(pointFaces, patchPointi)
600 // {
601 // const labelList& pFaces = pointFaces[patchPointi];
602 //
603 // labelHashSet pointCells(pFaces.size());
604 //
605 // forAll(pFaces, i)
606 // {
607 // label celli = mesh.faceOwner()[pp.addressing()[pFaces[i]]];
608 //
609 // if (!pointCells.insert(celli))
610 // {
611 // // Second or more occurrence of cell so cell has two or more
612 // // pp faces connected to this point.
613 // multiPatchCells.insert(celli);
614 // }
615 // }
616 // }
617 //
618 // label nMultiPatchCells = returnReduce
619 // (
620 // multiPatchCells.size(),
621 // sumOp<label>()
622 // );
623 //
624 // Info<< "Detected " << nMultiPatchCells
625 // << " cells with multiple (connected) patch faces." << endl;
626 //
627 // label nChanged = 0;
628 //
629 // if (nMultiPatchCells > 0)
630 // {
631 // multiPatchCells.instance() = meshRefiner_.timeName();
632 // Info<< "Writing " << nMultiPatchCells
633 // << " cells with multiple (connected) patch faces to cellSet "
634 // << multiPatchCells.objectPath() << endl;
635 // multiPatchCells.write();
636 //
637 //
638 // // Go through all points and remove extrusion on any cell in
639 // // multiPatchCells
640 // // (has to be done in separate loop since having one point on
641 // // multipatches has to reset extrusion on all points of cell)
642 //
643 // forAll(pointFaces, patchPointi)
644 // {
645 // if (extrudeStatus[patchPointi] != NOEXTRUDE)
646 // {
647 // const labelList& pFaces = pointFaces[patchPointi];
648 //
649 // forAll(pFaces, i)
650 // {
651 // label celli =
652 // mesh.faceOwner()[pp.addressing()[pFaces[i]]];
653 //
654 // if (multiPatchCells.found(celli))
655 // {
656 // if
657 // (
658 // unmarkExtrusion
659 // (
660 // patchPointi,
661 // patchDisp,
662 // patchNLayers,
663 // extrudeStatus
664 // )
665 // )
666 // {
667 // nChanged++;
668 // }
669 // }
670 // }
671 // }
672 // }
673 //
674 // reduce(nChanged, sumOp<label>());
675 // }
676 //
677 // Info<< "Prevented extrusion on " << nChanged
678 // << " points due to multiple patch faces." << nl << endl;
679 //}
680 
681 
682 void Foam::snappyLayerDriver::setNumLayers
683 (
684  const labelList& patchToNLayers,
685  const labelList& patchIDs,
686  const indirectPrimitivePatch& pp,
687  pointField& patchDisp,
688  labelList& patchNLayers,
689  List<extrudeMode>& extrudeStatus,
690  label& nAddedCells
691 ) const
692 {
693  const fvMesh& mesh = meshRefiner_.mesh();
694 
695  Info<< nl << "Handling points with inconsistent layer specification ..."
696  << endl;
697 
698  // Get for every point (really only necessary on patch external points)
699  // the max and min of any patch faces using it.
700  labelList maxLayers(patchNLayers.size(), labelMin);
701  labelList minLayers(patchNLayers.size(), labelMax);
702 
703  forAll(patchIDs, i)
704  {
705  label patchi = patchIDs[i];
706 
707  const labelList& meshPoints = mesh.boundaryMesh()[patchi].meshPoints();
708 
709  label wantedLayers = patchToNLayers[patchi];
710 
711  forAll(meshPoints, patchPointi)
712  {
713  label ppPointi = pp.meshPointMap()[meshPoints[patchPointi]];
714 
715  maxLayers[ppPointi] = max(wantedLayers, maxLayers[ppPointi]);
716  minLayers[ppPointi] = min(wantedLayers, minLayers[ppPointi]);
717  }
718  }
719 
721  (
722  mesh,
723  pp.meshPoints(),
724  maxLayers,
725  maxEqOp<label>(),
726  labelMin // null value
727  );
729  (
730  mesh,
731  pp.meshPoints(),
732  minLayers,
733  minEqOp<label>(),
734  labelMax // null value
735  );
736 
737  // Unmark any point with different min and max
738  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
739 
740  //label nConflicts = 0;
741 
742  forAll(maxLayers, i)
743  {
744  if (maxLayers[i] == labelMin || minLayers[i] == labelMax)
745  {
747  << "Patchpoint:" << i << " coord:" << pp.localPoints()[i]
748  << " maxLayers:" << maxLayers
749  << " minLayers:" << minLayers
750  << abort(FatalError);
751  }
752  else if (maxLayers[i] == minLayers[i])
753  {
754  // Ok setting.
755  patchNLayers[i] = maxLayers[i];
756  }
757  else
758  {
759  // Inconsistent num layers between patch faces using point
760  //if
761  //(
762  // unmarkExtrusion
763  // (
764  // i,
765  // patchDisp,
766  // patchNLayers,
767  // extrudeStatus
768  // )
769  //)
770  //{
771  // nConflicts++;
772  //}
773  patchNLayers[i] = maxLayers[i];
774  }
775  }
776 
777 
778  // Calculate number of cells to create
779  nAddedCells = 0;
780  forAll(pp.localFaces(), facei)
781  {
782  const face& f = pp.localFaces()[facei];
783 
784  // Get max of extrusion per point
785  label nCells = 0;
786  forAll(f, fp)
787  {
788  nCells = max(nCells, patchNLayers[f[fp]]);
789  }
790 
791  nAddedCells += nCells;
792  }
793  reduce(nAddedCells, sumOp<label>());
794 
795  //reduce(nConflicts, sumOp<label>());
796  //
797  //Info<< "Set displacement to zero for " << nConflicts
798  // << " points due to points being on multiple regions"
799  // << " with inconsistent nLayers specification." << endl;
800 }
801 
802 
803 // Construct pointVectorField with correct boundary conditions for adding
804 // layers
806 Foam::snappyLayerDriver::makeLayerDisplacementField
807 (
808  const pointMesh& pMesh,
809  const labelList& numLayers
810 )
811 {
812  // Construct displacement field.
813  const pointBoundaryMesh& pointPatches = pMesh.boundary();
814 
815  wordList patchFieldTypes
816  (
817  pointPatches.size(),
818  slipPointPatchVectorField::typeName
819  );
820  wordList actualPatchTypes(patchFieldTypes.size());
821  forAll(pointPatches, patchi)
822  {
823  actualPatchTypes[patchi] = pointPatches[patchi].type();
824  }
825 
826  forAll(numLayers, patchi)
827  {
828  // 0 layers: do not allow slip so fixedValue 0
829  // >0 layers: fixedValue which gets adapted
830  if (numLayers[patchi] == 0)
831  {
832  patchFieldTypes[patchi] =
833  zeroFixedValuePointPatchVectorField::typeName;
834  }
835  else if (numLayers[patchi] > 0)
836  {
837  patchFieldTypes[patchi] = fixedValuePointPatchVectorField::typeName;
838  }
839  }
840 
841  forAll(pointPatches, patchi)
842  {
843  if (isA<processorPointPatch>(pointPatches[patchi]))
844  {
845  patchFieldTypes[patchi] = calculatedPointPatchVectorField::typeName;
846  }
847  else if (isA<cyclicPointPatch>(pointPatches[patchi]))
848  {
849  patchFieldTypes[patchi] = cyclicSlipPointPatchVectorField::typeName;
850  }
851  }
852 
853 
854  const polyMesh& mesh = pMesh();
855 
856  // Note: time().timeName() instead of meshRefinement::timeName() since
857  // postprocessable field.
858  tmp<pointVectorField> tfld
859  (
860  new pointVectorField
861  (
862  IOobject
863  (
864  "pointDisplacement",
865  mesh.time().timeName(),
866  mesh,
869  ),
870  pMesh,
871  dimensionedVector("displacement", dimLength, Zero),
872  patchFieldTypes,
873  actualPatchTypes
874  )
875  );
876  return tfld;
877 }
878 
879 
880 void Foam::snappyLayerDriver::growNoExtrusion
881 (
882  const indirectPrimitivePatch& pp,
883  pointField& patchDisp,
884  labelList& patchNLayers,
885  List<extrudeMode>& extrudeStatus
886 ) const
887 {
888  Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
889 
890  List<extrudeMode> grownExtrudeStatus(extrudeStatus);
891 
892  const faceList& localFaces = pp.localFaces();
893 
894  label nGrown = 0;
895 
896  forAll(localFaces, facei)
897  {
898  const face& f = localFaces[facei];
899 
900  bool hasSqueeze = false;
901  forAll(f, fp)
902  {
903  if (extrudeStatus[f[fp]] == NOEXTRUDE)
904  {
905  hasSqueeze = true;
906  break;
907  }
908  }
909 
910  if (hasSqueeze)
911  {
912  // Squeeze all points of face
913  forAll(f, fp)
914  {
915  if
916  (
917  extrudeStatus[f[fp]] == EXTRUDE
918  && grownExtrudeStatus[f[fp]] != NOEXTRUDE
919  )
920  {
921  grownExtrudeStatus[f[fp]] = NOEXTRUDE;
922  nGrown++;
923  }
924  }
925  }
926  }
927 
928  extrudeStatus.transfer(grownExtrudeStatus);
929 
930 
931  // Synchronise since might get called multiple times.
932  // Use the fact that NOEXTRUDE is the minimum value.
933  {
934  labelList status(extrudeStatus.size());
935  forAll(status, i)
936  {
937  status[i] = extrudeStatus[i];
938  }
940  (
941  meshRefiner_.mesh(),
942  pp.meshPoints(),
943  status,
944  minEqOp<label>(),
945  labelMax // null value
946  );
947  forAll(status, i)
948  {
949  extrudeStatus[i] = extrudeMode(status[i]);
950  }
951  }
952 
953 
954  forAll(extrudeStatus, patchPointi)
955  {
956  if (extrudeStatus[patchPointi] == NOEXTRUDE)
957  {
958  patchDisp[patchPointi] = Zero;
959  patchNLayers[patchPointi] = 0;
960  }
961  }
962 
963  reduce(nGrown, sumOp<label>());
964 
965  Info<< "Set displacement to zero for an additional " << nGrown
966  << " points." << endl;
967 }
968 
969 
970 void Foam::snappyLayerDriver::determineSidePatches
971 (
972  const globalIndex& globalFaces,
973  const labelListList& edgeGlobalFaces,
974  const indirectPrimitivePatch& pp,
975 
976  labelList& sidePatchID
977 )
978 {
979  // Sometimes edges-to-be-extruded are on more than 2 processors.
980  // Work out which 2 hold the faces to be extruded and thus which procpatch
981  // the side-face should be in. As an additional complication this might
982  // mean that 2 procesors that were only edge-connected now suddenly need
983  // to become face-connected i.e. have a processor patch between them.
984 
985  fvMesh& mesh = meshRefiner_.mesh();
986 
987  // Determine sidePatchID. Any additional processor boundary gets added to
988  // patchToNbrProc,nbrProcToPatch and nPatches gets set to the new number
989  // of patches.
990  label nPatches;
991  Map<label> nbrProcToPatch;
992  Map<label> patchToNbrProc;
994  (
995  mesh,
996  globalFaces,
997  edgeGlobalFaces,
998  pp,
999 
1000  sidePatchID,
1001  nPatches,
1002  nbrProcToPatch,
1003  patchToNbrProc
1004  );
1005 
1006  label nOldPatches = mesh.boundaryMesh().size();
1007  label nAdded = returnReduce(nPatches-nOldPatches, sumOp<label>());
1008  Info<< nl << "Adding in total " << nAdded/2 << " inter-processor patches to"
1009  << " handle extrusion of non-manifold processor boundaries."
1010  << endl;
1011 
1012  if (nAdded > 0)
1013  {
1014  // We might not add patches in same order as in patchToNbrProc
1015  // so prepare to renumber sidePatchID
1016  Map<label> wantedToAddedPatch;
1017 
1018  for (label patchi = nOldPatches; patchi < nPatches; patchi++)
1019  {
1020  label nbrProci = patchToNbrProc[patchi];
1021  word name
1022  (
1024  );
1025 
1026  dictionary patchDict;
1027  patchDict.add("type", processorPolyPatch::typeName);
1028  patchDict.add("myProcNo", Pstream::myProcNo());
1029  patchDict.add("neighbProcNo", nbrProci);
1030  patchDict.add("nFaces", 0);
1031  patchDict.add("startFace", mesh.nFaces());
1032 
1033  label procPatchi = meshRefiner_.appendPatch
1034  (
1035  mesh,
1036  mesh.boundaryMesh().size(), // new patch index
1037  name,
1038  patchDict
1039  );
1040  wantedToAddedPatch.insert(patchi, procPatchi);
1041  }
1042 
1043  // Renumber sidePatchID
1044  forAll(sidePatchID, i)
1045  {
1046  label patchi = sidePatchID[i];
1047  Map<label>::const_iterator fnd = wantedToAddedPatch.find(patchi);
1048  if (fnd != wantedToAddedPatch.end())
1049  {
1050  sidePatchID[i] = fnd();
1051  }
1052  }
1053 
1054  mesh.clearOut();
1055  const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()).updateMesh();
1056  }
1057 }
1058 
1059 
1060 void Foam::snappyLayerDriver::calculateLayerThickness
1061 (
1062  const indirectPrimitivePatch& pp,
1063  const labelList& patchIDs,
1064  const layerParameters& layerParams,
1065  const labelList& cellLevel,
1066  const labelList& patchNLayers,
1067  const scalar edge0Len,
1068 
1069  scalarField& thickness,
1070  scalarField& minThickness,
1071  scalarField& expansionRatio
1072 ) const
1073 {
1074  const fvMesh& mesh = meshRefiner_.mesh();
1075  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1076 
1077 
1078  // Rework patch-wise layer parameters into minimum per point
1079  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1080  // Note: only layer parameters consistent with layer specification
1081  // method (see layerParameters) will be correct.
1082  scalarField firstLayerThickness(pp.nPoints(), GREAT);
1083  scalarField finalLayerThickness(pp.nPoints(), GREAT);
1084  scalarField totalThickness(pp.nPoints(), GREAT);
1085  scalarField expRatio(pp.nPoints(), GREAT);
1086 
1087  minThickness.setSize(pp.nPoints());
1088  minThickness = GREAT;
1089 
1090  forAll(patchIDs, i)
1091  {
1092  label patchi = patchIDs[i];
1093 
1094  const labelList& meshPoints = patches[patchi].meshPoints();
1095 
1096  forAll(meshPoints, patchPointi)
1097  {
1098  label ppPointi = pp.meshPointMap()[meshPoints[patchPointi]];
1099 
1100  firstLayerThickness[ppPointi] = min
1101  (
1102  firstLayerThickness[ppPointi],
1103  layerParams.firstLayerThickness()[patchi]
1104  );
1105  finalLayerThickness[ppPointi] = min
1106  (
1107  finalLayerThickness[ppPointi],
1108  layerParams.finalLayerThickness()[patchi]
1109  );
1110  totalThickness[ppPointi] = min
1111  (
1112  totalThickness[ppPointi],
1113  layerParams.thickness()[patchi]
1114  );
1115  expRatio[ppPointi] = min
1116  (
1117  expRatio[ppPointi],
1118  layerParams.expansionRatio()[patchi]
1119  );
1120  minThickness[ppPointi] = min
1121  (
1122  minThickness[ppPointi],
1123  layerParams.minThickness()[patchi]
1124  );
1125  }
1126  }
1127 
1129  (
1130  mesh,
1131  pp.meshPoints(),
1132  firstLayerThickness,
1133  minEqOp<scalar>(),
1134  GREAT // null value
1135  );
1137  (
1138  mesh,
1139  pp.meshPoints(),
1140  finalLayerThickness,
1141  minEqOp<scalar>(),
1142  GREAT // null value
1143  );
1145  (
1146  mesh,
1147  pp.meshPoints(),
1148  totalThickness,
1149  minEqOp<scalar>(),
1150  GREAT // null value
1151  );
1153  (
1154  mesh,
1155  pp.meshPoints(),
1156  expRatio,
1157  minEqOp<scalar>(),
1158  GREAT // null value
1159  );
1161  (
1162  mesh,
1163  pp.meshPoints(),
1164  minThickness,
1165  minEqOp<scalar>(),
1166  GREAT // null value
1167  );
1168 
1169 
1170  // Now the thicknesses are set according to the minimum of connected
1171  // patches.
1172 
1173 
1174  // Rework relative thickness into absolute
1175  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1176  // by multiplying with the internal cell size.
1177 
1178  if (layerParams.relativeSizes())
1179  {
1180  if
1181  (
1182  min(layerParams.minThickness()) < 0
1183  || max(layerParams.minThickness()) > 2
1184  )
1185  {
1187  << "Thickness should be factor of local undistorted cell size."
1188  << " Valid values are [0..2]." << nl
1189  << " minThickness:" << layerParams.minThickness()
1190  << exit(FatalError);
1191  }
1192 
1193 
1194  // Determine per point the max cell level of connected cells
1195  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1196 
1197  labelList maxPointLevel(pp.nPoints(), labelMin);
1198 
1199  forAll(pp, i)
1200  {
1201  label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
1202 
1203  const face& f = pp.localFaces()[i];
1204 
1205  forAll(f, fp)
1206  {
1207  maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
1208  }
1209  }
1210 
1212  (
1213  mesh,
1214  pp.meshPoints(),
1215  maxPointLevel,
1216  maxEqOp<label>(),
1217  labelMin // null value
1218  );
1219 
1220 
1221  forAll(maxPointLevel, pointi)
1222  {
1223  // Find undistorted edge size for this level.
1224  scalar edgeLen = edge0Len/(1<<maxPointLevel[pointi]);
1225  firstLayerThickness[pointi] *= edgeLen;
1226  finalLayerThickness[pointi] *= edgeLen;
1227  totalThickness[pointi] *= edgeLen;
1228  minThickness[pointi] *= edgeLen;
1229  }
1230  }
1231 
1232 
1233 
1234  // Rework thickness parameters into overall thickness
1235  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1236 
1237  forAll(firstLayerThickness, pointi)
1238  {
1239  thickness[pointi] = layerParams.layerThickness
1240  (
1241  patchNLayers[pointi],
1242  firstLayerThickness[pointi],
1243  finalLayerThickness[pointi],
1244  totalThickness[pointi],
1245  expRatio[pointi]
1246  );
1247 
1248  expansionRatio[pointi] = layerParams.layerExpansionRatio
1249  (
1250  patchNLayers[pointi],
1251  firstLayerThickness[pointi],
1252  finalLayerThickness[pointi],
1253  totalThickness[pointi],
1254  expRatio[pointi]
1255  );
1256  }
1257 
1258  //Info<< "calculateLayerThickness : min:" << gMin(thickness)
1259  // << " max:" << gMax(thickness) << endl;
1260 
1261  // Print a bit
1262  {
1263  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1264 
1265  int oldPrecision = Info().precision();
1266 
1267  // Find maximum length of a patch name, for a nicer output
1268  label maxPatchNameLen = 0;
1269  forAll(patchIDs, i)
1270  {
1271  label patchi = patchIDs[i];
1272  word patchName = patches[patchi].name();
1273  maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
1274  }
1275 
1276  Info<< nl
1277  << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
1278  << setw(0) << " faces layers avg thickness[m]" << nl
1279  << setf(ios_base::left) << setw(maxPatchNameLen) << " "
1280  << setw(0) << " near-wall overall" << nl
1281  << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
1282  << setw(0) << " ----- ------ --------- -------" << endl;
1283 
1284 
1285  const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
1286 
1287  forAll(patchIDs, i)
1288  {
1289  label patchi = patchIDs[i];
1290 
1291  const labelList& meshPoints = patches[patchi].meshPoints();
1292 
1293  scalar sumThickness = 0;
1294  scalar sumNearWallThickness = 0;
1295  label nMasterPoints = 0;
1296 
1297  forAll(meshPoints, patchPointi)
1298  {
1299  label meshPointi = meshPoints[patchPointi];
1300  if (isMasterPoint[meshPointi])
1301  {
1302  label ppPointi = pp.meshPointMap()[meshPointi];
1303 
1304  sumThickness += thickness[ppPointi];
1305  sumNearWallThickness += layerParams.firstLayerThickness
1306  (
1307  patchNLayers[ppPointi],
1308  firstLayerThickness[ppPointi],
1309  finalLayerThickness[ppPointi],
1310  thickness[ppPointi],
1311  expansionRatio[ppPointi]
1312  );
1313  nMasterPoints++;
1314  }
1315  }
1316 
1317  label totNPoints = returnReduce(nMasterPoints, sumOp<label>());
1318 
1319  // For empty patches, totNPoints is 0.
1320  scalar avgThickness = 0;
1321  scalar avgNearWallThickness = 0;
1322 
1323  if (totNPoints > 0)
1324  {
1325  avgThickness =
1326  returnReduce(sumThickness, sumOp<scalar>())
1327  / totNPoints;
1328  avgNearWallThickness =
1329  returnReduce(sumNearWallThickness, sumOp<scalar>())
1330  / totNPoints;
1331  }
1332 
1333  Info<< setf(ios_base::left) << setw(maxPatchNameLen)
1334  << patches[patchi].name() << setprecision(3)
1335  << " " << setw(8)
1336  << returnReduce(patches[patchi].size(), sumOp<scalar>())
1337  << " " << setw(6) << layerParams.numLayers()[patchi]
1338  << " " << setw(8) << avgNearWallThickness
1339  << " " << setw(8) << avgThickness
1340  << endl;
1341  }
1342  Info<< setprecision(oldPrecision) << endl;
1343  }
1344 }
1345 
1346 
1347 // Synchronize displacement among coupled patches.
1348 void Foam::snappyLayerDriver::syncPatchDisplacement
1349 (
1350  const indirectPrimitivePatch& pp,
1351  const scalarField& minThickness,
1352  pointField& patchDisp,
1353  labelList& patchNLayers,
1354  List<extrudeMode>& extrudeStatus
1355 ) const
1356 {
1357  const fvMesh& mesh = meshRefiner_.mesh();
1358  const labelList& meshPoints = pp.meshPoints();
1359 
1360  label nChangedTotal = 0;
1361 
1362  while (true)
1363  {
1364  label nChanged = 0;
1365 
1366  // Sync displacement (by taking min)
1368  (
1369  mesh,
1370  meshPoints,
1371  patchDisp,
1372  minMagSqrEqOp<vector>(),
1373  point::rootMax // null value
1374  );
1375 
1376  // Unmark if displacement too small
1377  forAll(patchDisp, i)
1378  {
1379  if (mag(patchDisp[i]) < minThickness[i])
1380  {
1381  if
1382  (
1383  unmarkExtrusion
1384  (
1385  i,
1386  patchDisp,
1387  patchNLayers,
1388  extrudeStatus
1389  )
1390  )
1391  {
1392  nChanged++;
1393  }
1394  }
1395  }
1396 
1397  labelList syncPatchNLayers(patchNLayers);
1398 
1400  (
1401  mesh,
1402  meshPoints,
1403  syncPatchNLayers,
1404  minEqOp<label>(),
1405  labelMax // null value
1406  );
1407 
1408  // Reset if differs
1409  // 1. take max
1410  forAll(syncPatchNLayers, i)
1411  {
1412  if (syncPatchNLayers[i] != patchNLayers[i])
1413  {
1414  if
1415  (
1416  unmarkExtrusion
1417  (
1418  i,
1419  patchDisp,
1420  patchNLayers,
1421  extrudeStatus
1422  )
1423  )
1424  {
1425  nChanged++;
1426  }
1427  }
1428  }
1429 
1431  (
1432  mesh,
1433  meshPoints,
1434  syncPatchNLayers,
1435  maxEqOp<label>(),
1436  labelMin // null value
1437  );
1438 
1439  // Reset if differs
1440  // 2. take min
1441  forAll(syncPatchNLayers, i)
1442  {
1443  if (syncPatchNLayers[i] != patchNLayers[i])
1444  {
1445  if
1446  (
1447  unmarkExtrusion
1448  (
1449  i,
1450  patchDisp,
1451  patchNLayers,
1452  extrudeStatus
1453  )
1454  )
1455  {
1456  nChanged++;
1457  }
1458  }
1459  }
1460  nChangedTotal += nChanged;
1461 
1462  if (!returnReduce(nChanged, sumOp<label>()))
1463  {
1464  break;
1465  }
1466  }
1467 
1468  //Info<< "Prevented extrusion on "
1469  // << returnReduce(nChangedTotal, sumOp<label>())
1470  // << " coupled patch points during syncPatchDisplacement." << endl;
1471 }
1472 
1473 
1474 // Calculate displacement vector for all patch points. Uses pointNormal.
1475 // Checks that displaced patch point would be visible from all centres
1476 // of the faces using it.
1477 // extrudeStatus is both input and output and gives the status of each
1478 // patch point.
1479 void Foam::snappyLayerDriver::getPatchDisplacement
1480 (
1481  const indirectPrimitivePatch& pp,
1482  const scalarField& thickness,
1483  const scalarField& minThickness,
1484  pointField& patchDisp,
1485  labelList& patchNLayers,
1486  List<extrudeMode>& extrudeStatus
1487 ) const
1488 {
1489  Info<< nl << "Determining displacement for added points"
1490  << " according to pointNormal ..." << endl;
1491 
1492  const fvMesh& mesh = meshRefiner_.mesh();
1493  const vectorField& faceNormals = pp.faceNormals();
1494  const labelListList& pointFaces = pp.pointFaces();
1495  const pointField& localPoints = pp.localPoints();
1496 
1497  // Determine pointNormal
1498  // ~~~~~~~~~~~~~~~~~~~~~
1499 
1500  pointField pointNormals(PatchTools::pointNormals(mesh, pp));
1501 
1502 
1503  // Determine local length scale on patch
1504  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1505 
1506  // Start off from same thickness everywhere (except where no extrusion)
1507  patchDisp = thickness*pointNormals;
1508 
1509 
1510  label nNoVisNormal = 0;
1511  label nExtrudeRemove = 0;
1512 
1513 
1514  // Check if no extrude possible.
1515  forAll(pointNormals, patchPointi)
1516  {
1517  label meshPointi = pp.meshPoints()[patchPointi];
1518 
1519  if (extrudeStatus[patchPointi] == NOEXTRUDE)
1520  {
1521  // Do not use unmarkExtrusion; forcibly set to zero extrusion.
1522  patchNLayers[patchPointi] = 0;
1523  patchDisp[patchPointi] = Zero;
1524  }
1525  else
1526  {
1527  // Get normal
1528  const vector& n = pointNormals[patchPointi];
1529 
1530  if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointi]))
1531  {
1532  if (debug&meshRefinement::ATTRACTION)
1533  {
1534  Pout<< "No valid normal for point " << meshPointi
1535  << ' ' << pp.points()[meshPointi]
1536  << "; setting displacement to "
1537  << patchDisp[patchPointi]
1538  << endl;
1539  }
1540 
1541  extrudeStatus[patchPointi] = EXTRUDEREMOVE;
1542  nNoVisNormal++;
1543  }
1544  }
1545  }
1546 
1547  // At illegal points make displacement average of new neighbour positions
1548  forAll(extrudeStatus, patchPointi)
1549  {
1550  if (extrudeStatus[patchPointi] == EXTRUDEREMOVE)
1551  {
1552  point avg(Zero);
1553  label nPoints = 0;
1554 
1555  const labelList& pEdges = pp.pointEdges()[patchPointi];
1556 
1557  forAll(pEdges, i)
1558  {
1559  label edgeI = pEdges[i];
1560 
1561  label otherPointi = pp.edges()[edgeI].otherVertex(patchPointi);
1562 
1563  if (extrudeStatus[otherPointi] != NOEXTRUDE)
1564  {
1565  avg += localPoints[otherPointi] + patchDisp[otherPointi];
1566  nPoints++;
1567  }
1568  }
1569 
1570  if (nPoints > 0)
1571  {
1572  if (debug&meshRefinement::ATTRACTION)
1573  {
1574  Pout<< "Displacement at illegal point "
1575  << localPoints[patchPointi]
1576  << " set to "
1577  << (avg / nPoints - localPoints[patchPointi])
1578  << endl;
1579  }
1580 
1581  patchDisp[patchPointi] =
1582  avg / nPoints
1583  - localPoints[patchPointi];
1584 
1585  nExtrudeRemove++;
1586  }
1587  else
1588  {
1589  // All surrounding points are not extruded. Leave patchDisp
1590  // intact.
1591  }
1592  }
1593  }
1594 
1595  Info<< "Detected " << returnReduce(nNoVisNormal, sumOp<label>())
1596  << " points with point normal pointing through faces." << nl
1597  << "Reset displacement at "
1598  << returnReduce(nExtrudeRemove, sumOp<label>())
1599  << " points to average of surrounding points." << endl;
1600 
1601  // Make sure displacement is equal on both sides of coupled patches.
1602  syncPatchDisplacement
1603  (
1604  pp,
1605  minThickness,
1606  patchDisp,
1607  patchNLayers,
1608  extrudeStatus
1609  );
1610 
1611  Info<< endl;
1612 }
1613 
1614 
1615 bool Foam::snappyLayerDriver::sameEdgeNeighbour
1616 (
1617  const labelListList& globalEdgeFaces,
1618  const label myGlobalFacei,
1619  const label nbrGlobFacei,
1620  const label edgeI
1621 ) const
1622 {
1623  const labelList& eFaces = globalEdgeFaces[edgeI];
1624  if (eFaces.size() == 2)
1625  {
1626  return edge(myGlobalFacei, nbrGlobFacei) == edge(eFaces[0], eFaces[1]);
1627  }
1628  else
1629  {
1630  return false;
1631  }
1632 }
1633 
1634 
1635 void Foam::snappyLayerDriver::getVertexString
1636 (
1637  const indirectPrimitivePatch& pp,
1638  const labelListList& globalEdgeFaces,
1639  const label facei,
1640  const label edgeI,
1641  const label myGlobFacei,
1642  const label nbrGlobFacei,
1643  DynamicList<label>& vertices
1644 ) const
1645 {
1646  const labelList& fEdges = pp.faceEdges()[facei];
1647  label fp = findIndex(fEdges, edgeI);
1648 
1649  if (fp == -1)
1650  {
1652  << "problem." << abort(FatalError);
1653  }
1654 
1655  // Search back
1656  label startFp = fp;
1657 
1658  forAll(fEdges, i)
1659  {
1660  label prevFp = fEdges.rcIndex(startFp);
1661  if
1662  (
1663  !sameEdgeNeighbour
1664  (
1665  globalEdgeFaces,
1666  myGlobFacei,
1667  nbrGlobFacei,
1668  fEdges[prevFp]
1669  )
1670  )
1671  {
1672  break;
1673  }
1674  startFp = prevFp;
1675  }
1676 
1677  label endFp = fp;
1678  forAll(fEdges, i)
1679  {
1680  label nextFp = fEdges.fcIndex(endFp);
1681  if
1682  (
1683  !sameEdgeNeighbour
1684  (
1685  globalEdgeFaces,
1686  myGlobFacei,
1687  nbrGlobFacei,
1688  fEdges[nextFp]
1689  )
1690  )
1691  {
1692  break;
1693  }
1694  endFp = nextFp;
1695  }
1696 
1697  const face& f = pp.localFaces()[facei];
1698  vertices.clear();
1699  fp = startFp;
1700  while (fp != endFp)
1701  {
1702  vertices.append(f[fp]);
1703  fp = f.fcIndex(fp);
1704  }
1705  vertices.append(f[fp]);
1706  fp = f.fcIndex(fp);
1707  vertices.append(f[fp]);
1708 }
1709 
1710 
1711 // Truncates displacement
1712 // - for all patchFaces in the faceset displacement gets set to zero
1713 // - all displacement < minThickness gets set to zero
1714 Foam::label Foam::snappyLayerDriver::truncateDisplacement
1715 (
1716  const globalIndex& globalFaces,
1717  const labelListList& edgeGlobalFaces,
1718  const indirectPrimitivePatch& pp,
1719  const scalarField& minThickness,
1720  const faceSet& illegalPatchFaces,
1721  pointField& patchDisp,
1722  labelList& patchNLayers,
1723  List<extrudeMode>& extrudeStatus
1724 ) const
1725 {
1726  const fvMesh& mesh = meshRefiner_.mesh();
1727 
1728  label nChanged = 0;
1729 
1730  const Map<label>& meshPointMap = pp.meshPointMap();
1731 
1732  forAllConstIter(faceSet, illegalPatchFaces, iter)
1733  {
1734  label facei = iter.key();
1735 
1736  if (mesh.isInternalFace(facei))
1737  {
1739  << "Faceset " << illegalPatchFaces.name()
1740  << " contains internal face " << facei << nl
1741  << "It should only contain patch faces" << abort(FatalError);
1742  }
1743 
1744  const face& f = mesh.faces()[facei];
1745 
1746 
1747  forAll(f, fp)
1748  {
1749  if (meshPointMap.found(f[fp]))
1750  {
1751  label patchPointi = meshPointMap[f[fp]];
1752 
1753  if (extrudeStatus[patchPointi] != NOEXTRUDE)
1754  {
1755  unmarkExtrusion
1756  (
1757  patchPointi,
1758  patchDisp,
1759  patchNLayers,
1760  extrudeStatus
1761  );
1762  nChanged++;
1763  }
1764  }
1765  }
1766  }
1767 
1768  forAll(patchDisp, patchPointi)
1769  {
1770  if (mag(patchDisp[patchPointi]) < minThickness[patchPointi])
1771  {
1772  if
1773  (
1774  unmarkExtrusion
1775  (
1776  patchPointi,
1777  patchDisp,
1778  patchNLayers,
1779  extrudeStatus
1780  )
1781  )
1782  {
1783  nChanged++;
1784  }
1785  }
1786  else if (extrudeStatus[patchPointi] == NOEXTRUDE)
1787  {
1788  // Make sure displacement is 0. Should already be so but ...
1789  patchDisp[patchPointi] = Zero;
1790  patchNLayers[patchPointi] = 0;
1791  }
1792  }
1793 
1794 
1795  const faceList& localFaces = pp.localFaces();
1796 
1797  while (true)
1798  {
1799  syncPatchDisplacement
1800  (
1801  pp,
1802  minThickness,
1803  patchDisp,
1804  patchNLayers,
1805  extrudeStatus
1806  );
1807 
1808 
1809  // Pinch
1810  // ~~~~~
1811 
1812  // Make sure that a face doesn't have two non-consecutive areas
1813  // not extruded (e.g. quad where vertex 0 and 2 are not extruded
1814  // but 1 and 3 are) since this gives topological errors.
1815 
1816  label nPinched = 0;
1817 
1818  forAll(localFaces, i)
1819  {
1820  const face& localF = localFaces[i];
1821 
1822  // Count number of transitions from unsnapped to snapped.
1823  label nTrans = 0;
1824 
1825  extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
1826 
1827  forAll(localF, fp)
1828  {
1829  extrudeMode fpMode = extrudeStatus[localF[fp]];
1830 
1831  if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
1832  {
1833  nTrans++;
1834  }
1835  prevMode = fpMode;
1836  }
1837 
1838  if (nTrans > 1)
1839  {
1840  // Multiple pinches. Reset whole face as unextruded.
1841  if
1842  (
1843  unmarkExtrusion
1844  (
1845  localF,
1846  patchDisp,
1847  patchNLayers,
1848  extrudeStatus
1849  )
1850  )
1851  {
1852  nPinched++;
1853  nChanged++;
1854  }
1855  }
1856  }
1857 
1858  reduce(nPinched, sumOp<label>());
1859 
1860  Info<< "truncateDisplacement : Unextruded " << nPinched
1861  << " faces due to non-consecutive vertices being extruded." << endl;
1862 
1863 
1864  // Butterfly
1865  // ~~~~~~~~~
1866 
1867  // Make sure that a string of edges becomes a single face so
1868  // not a butterfly. Occassionally an 'edge' will have a single dangling
1869  // vertex due to face combining. These get extruded as a single face
1870  // (with a dangling vertex) so make sure this extrusion forms a single
1871  // shape.
1872  // - continuous i.e. no butterfly:
1873  // + +
1874  // |\ /|
1875  // | \ / |
1876  // +--+--+
1877  // - extrudes from all but the endpoints i.e. no partial
1878  // extrude
1879  // +
1880  // /|
1881  // / |
1882  // +--+--+
1883  // The common error topology is a pinch somewhere in the middle
1884  label nButterFly = 0;
1885  {
1886  DynamicList<label> stringedVerts;
1887  forAll(pp.edges(), edgeI)
1888  {
1889  const labelList& globFaces = edgeGlobalFaces[edgeI];
1890 
1891  if (globFaces.size() == 2)
1892  {
1893  label myFacei = pp.edgeFaces()[edgeI][0];
1894  label myGlobalFacei = globalFaces.toGlobal
1895  (
1896  pp.addressing()[myFacei]
1897  );
1898  label nbrGlobalFacei =
1899  (
1900  globFaces[0] != myGlobalFacei
1901  ? globFaces[0]
1902  : globFaces[1]
1903  );
1904  getVertexString
1905  (
1906  pp,
1907  edgeGlobalFaces,
1908  myFacei,
1909  edgeI,
1910  myGlobalFacei,
1911  nbrGlobalFacei,
1912  stringedVerts
1913  );
1914 
1915  if
1916  (
1917  extrudeStatus[stringedVerts[0]] != NOEXTRUDE
1918  || extrudeStatus[stringedVerts.last()] != NOEXTRUDE
1919  )
1920  {
1921  // Any pinch in the middle
1922  bool pinch = false;
1923  for (label i = 1; i < stringedVerts.size()-1; i++)
1924  {
1925  if
1926  (
1927  extrudeStatus[stringedVerts[i]] == NOEXTRUDE
1928  )
1929  {
1930  pinch = true;
1931  break;
1932  }
1933  }
1934  if (pinch)
1935  {
1936  forAll(stringedVerts, i)
1937  {
1938  if
1939  (
1940  unmarkExtrusion
1941  (
1942  stringedVerts[i],
1943  patchDisp,
1944  patchNLayers,
1945  extrudeStatus
1946  )
1947  )
1948  {
1949  nButterFly++;
1950  nChanged++;
1951  }
1952  }
1953  }
1954  }
1955  }
1956  }
1957  }
1958 
1959  reduce(nButterFly, sumOp<label>());
1960 
1961  Info<< "truncateDisplacement : Unextruded " << nButterFly
1962  << " faces due to stringed edges with inconsistent extrusion."
1963  << endl;
1964 
1965 
1966 
1967  // Consistent number of layers
1968  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1969 
1970  // Make sure that a face has consistent number of layers for all
1971  // its vertices.
1972 
1973  label nDiffering = 0;
1974 
1975  //forAll(localFaces, i)
1976  //{
1977  // const face& localF = localFaces[i];
1978  //
1979  // label numLayers = -1;
1980  //
1981  // forAll(localF, fp)
1982  // {
1983  // if (patchNLayers[localF[fp]] > 0)
1984  // {
1985  // if (numLayers == -1)
1986  // {
1987  // numLayers = patchNLayers[localF[fp]];
1988  // }
1989  // else if (numLayers != patchNLayers[localF[fp]])
1990  // {
1991  // // Differing number of layers
1992  // if
1993  // (
1994  // unmarkExtrusion
1995  // (
1996  // localF,
1997  // patchDisp,
1998  // patchNLayers,
1999  // extrudeStatus
2000  // )
2001  // )
2002  // {
2003  // nDiffering++;
2004  // nChanged++;
2005  // }
2006  // break;
2007  // }
2008  // }
2009  // }
2010  //}
2011  //
2012  //reduce(nDiffering, sumOp<label>());
2013  //
2014  //Info<< "truncateDisplacement : Unextruded " << nDiffering
2015  // << " faces due to having differing number of layers." << endl;
2016 
2017  if (nPinched+nButterFly+nDiffering == 0)
2018  {
2019  break;
2020  }
2021  }
2022 
2023  return nChanged;
2024 }
2025 
2026 
2027 // Setup layer information (at points and faces) to modify mesh topology in
2028 // regions where layer mesh terminates.
2029 void Foam::snappyLayerDriver::setupLayerInfoTruncation
2030 (
2031  const indirectPrimitivePatch& pp,
2032  const labelList& patchNLayers,
2033  const List<extrudeMode>& extrudeStatus,
2034  const label nBufferCellsNoExtrude,
2035  labelList& nPatchPointLayers,
2036  labelList& nPatchFaceLayers
2037 ) const
2038 {
2039  Info<< nl << "Setting up information for layer truncation ..." << endl;
2040 
2041  const fvMesh& mesh = meshRefiner_.mesh();
2042 
2043  if (nBufferCellsNoExtrude < 0)
2044  {
2045  Info<< nl << "Performing no layer truncation."
2046  << " nBufferCellsNoExtrude set to less than 0 ..." << endl;
2047 
2048  // Face layers if any point gets extruded
2049  forAll(pp.localFaces(), patchFacei)
2050  {
2051  const face& f = pp.localFaces()[patchFacei];
2052 
2053  forAll(f, fp)
2054  {
2055  if (patchNLayers[f[fp]] > 0)
2056  {
2057  nPatchFaceLayers[patchFacei] = patchNLayers[f[fp]];
2058  break;
2059  }
2060  }
2061  }
2062  nPatchPointLayers = patchNLayers;
2063 
2064  // Set any unset patch face layers
2065  forAll(nPatchFaceLayers, patchFacei)
2066  {
2067  if (nPatchFaceLayers[patchFacei] == -1)
2068  {
2069  nPatchFaceLayers[patchFacei] = 0;
2070  }
2071  }
2072  }
2073  else
2074  {
2075  // Determine max point layers per face.
2076  labelList maxLevel(pp.size(), 0);
2077 
2078  forAll(pp.localFaces(), patchFacei)
2079  {
2080  const face& f = pp.localFaces()[patchFacei];
2081 
2082  // find patch faces where layer terminates (i.e contains extrude
2083  // and noextrude points).
2084 
2085  bool noExtrude = false;
2086  label mLevel = 0;
2087 
2088  forAll(f, fp)
2089  {
2090  if (extrudeStatus[f[fp]] == NOEXTRUDE)
2091  {
2092  noExtrude = true;
2093  }
2094  mLevel = max(mLevel, patchNLayers[f[fp]]);
2095  }
2096 
2097  if (mLevel > 0)
2098  {
2099  // So one of the points is extruded. Check if all are extruded
2100  // or is a mix.
2101 
2102  if (noExtrude)
2103  {
2104  nPatchFaceLayers[patchFacei] = 1;
2105  maxLevel[patchFacei] = mLevel;
2106  }
2107  else
2108  {
2109  maxLevel[patchFacei] = mLevel;
2110  }
2111  }
2112  }
2113 
2114  // We have the seed faces (faces with nPatchFaceLayers != maxLevel)
2115  // Now do a meshwave across the patch where we pick up neighbours
2116  // of seed faces.
2117  // Note: quite inefficient. Could probably be coded better.
2118 
2119  const labelListList& pointFaces = pp.pointFaces();
2120 
2121  label nLevels = gMax(patchNLayers);
2122 
2123  // flag neighbouring patch faces with number of layers to grow
2124  for (label ilevel = 1; ilevel < nLevels; ilevel++)
2125  {
2126  label nBuffer;
2127 
2128  if (ilevel == 1)
2129  {
2130  nBuffer = nBufferCellsNoExtrude - 1;
2131  }
2132  else
2133  {
2134  nBuffer = nBufferCellsNoExtrude;
2135  }
2136 
2137  for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
2138  {
2139  labelList tempCounter(nPatchFaceLayers);
2140 
2141  boolList foundNeighbour(pp.nPoints(), false);
2142 
2143  forAll(pp.meshPoints(), patchPointi)
2144  {
2145  forAll(pointFaces[patchPointi], pointFacei)
2146  {
2147  label facei = pointFaces[patchPointi][pointFacei];
2148 
2149  if
2150  (
2151  nPatchFaceLayers[facei] != -1
2152  && maxLevel[facei] > 0
2153  )
2154  {
2155  foundNeighbour[patchPointi] = true;
2156  break;
2157  }
2158  }
2159  }
2160 
2162  (
2163  mesh,
2164  pp.meshPoints(),
2165  foundNeighbour,
2166  orEqOp<bool>(),
2167  false // null value
2168  );
2169 
2170  forAll(pp.meshPoints(), patchPointi)
2171  {
2172  if (foundNeighbour[patchPointi])
2173  {
2174  forAll(pointFaces[patchPointi], pointFacei)
2175  {
2176  label facei = pointFaces[patchPointi][pointFacei];
2177  if
2178  (
2179  nPatchFaceLayers[facei] == -1
2180  && maxLevel[facei] > 0
2181  && ilevel < maxLevel[facei]
2182  )
2183  {
2184  tempCounter[facei] = ilevel;
2185  }
2186  }
2187  }
2188  }
2189  nPatchFaceLayers = tempCounter;
2190  }
2191  }
2192 
2193  forAll(pp.localFaces(), patchFacei)
2194  {
2195  if (nPatchFaceLayers[patchFacei] == -1)
2196  {
2197  nPatchFaceLayers[patchFacei] = maxLevel[patchFacei];
2198  }
2199  }
2200 
2201  forAll(pp.meshPoints(), patchPointi)
2202  {
2203  if (extrudeStatus[patchPointi] != NOEXTRUDE)
2204  {
2205  forAll(pointFaces[patchPointi], pointFacei)
2206  {
2207  label face = pointFaces[patchPointi][pointFacei];
2208  nPatchPointLayers[patchPointi] = max
2209  (
2210  nPatchPointLayers[patchPointi],
2211  nPatchFaceLayers[face]
2212  );
2213  }
2214  }
2215  else
2216  {
2217  nPatchPointLayers[patchPointi] = 0;
2218  }
2219  }
2221  (
2222  mesh,
2223  pp.meshPoints(),
2224  nPatchPointLayers,
2225  maxEqOp<label>(),
2226  label(0) // null value
2227  );
2228  }
2229 }
2230 
2231 
2232 // Does any of the cells use a face from faces?
2233 bool Foam::snappyLayerDriver::cellsUseFace
2234 (
2235  const polyMesh& mesh,
2236  const labelList& cellLabels,
2237  const labelHashSet& faces
2238 )
2239 {
2240  forAll(cellLabels, i)
2241  {
2242  const cell& cFaces = mesh.cells()[cellLabels[i]];
2243 
2244  forAll(cFaces, cFacei)
2245  {
2246  if (faces.found(cFaces[cFacei]))
2247  {
2248  return true;
2249  }
2250  }
2251  }
2252  return false;
2253 }
2254 
2255 
2256 // Checks the newly added cells and locally unmarks points so they
2257 // will not get extruded next time round. Returns global number of unmarked
2258 // points (0 if all was fine)
2259 Foam::label Foam::snappyLayerDriver::checkAndUnmark
2260 (
2261  const addPatchCellLayer& addLayer,
2262  const dictionary& meshQualityDict,
2263  const bool additionalReporting,
2264  const List<labelPair>& baffles,
2265  const indirectPrimitivePatch& pp,
2266  const fvMesh& newMesh,
2267 
2268  pointField& patchDisp,
2269  labelList& patchNLayers,
2270  List<extrudeMode>& extrudeStatus
2271 )
2272 {
2273  // Check the resulting mesh for errors
2274  Info<< nl << "Checking mesh with layer ..." << endl;
2275  faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
2277  (
2278  false,
2279  newMesh,
2280  meshQualityDict,
2281  identity(newMesh.nFaces()),
2282  baffles,
2283  wrongFaces
2284  );
2285  Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
2286  << " illegal faces"
2287  << " (concave, zero area or negative cell pyramid volume)"
2288  << endl;
2289 
2290  // Undo local extrusion if
2291  // - any of the added cells in error
2292 
2293  label nChanged = 0;
2294 
2295  // Get all cells in the layer.
2296  labelListList addedCells
2297  (
2299  (
2300  newMesh,
2301  addLayer.layerFaces()
2302  )
2303  );
2304 
2305  // Check if any of the faces in error uses any face of an added cell
2306  // - if additionalReporting print the few remaining areas for ease of
2307  // finding out where the problems are.
2308 
2309  const label nReportMax = 10;
2310  DynamicField<point> disabledFaceCentres(nReportMax);
2311 
2312  forAll(addedCells, oldPatchFacei)
2313  {
2314  // Get the cells (in newMesh labels) per old patch face (in mesh
2315  // labels)
2316  const labelList& fCells = addedCells[oldPatchFacei];
2317 
2318  if (cellsUseFace(newMesh, fCells, wrongFaces))
2319  {
2320  // Unmark points on old mesh
2321  if
2322  (
2323  unmarkExtrusion
2324  (
2325  pp.localFaces()[oldPatchFacei],
2326  patchDisp,
2327  patchNLayers,
2328  extrudeStatus
2329  )
2330  )
2331  {
2332  if (additionalReporting && (nChanged < nReportMax))
2333  {
2334  disabledFaceCentres.append
2335  (
2336  pp.faceCentres()[oldPatchFacei]
2337  );
2338  }
2339 
2340  nChanged++;
2341  }
2342  }
2343  }
2344 
2345 
2346  label nChangedTotal = returnReduce(nChanged, sumOp<label>());
2347 
2348  if (additionalReporting)
2349  {
2350  // Limit the number of points to be printed so that
2351  // not too many points are reported when running in parallel
2352  // Not accurate, i.e. not always nReportMax points are written,
2353  // but this estimation avoid some communication here.
2354  // The important thing, however, is that when only a few faces
2355  // are disabled, their coordinates are printed, and this should be
2356  // the case
2357  label nReportLocal = nChanged;
2358  if (nChangedTotal > nReportMax)
2359  {
2360  nReportLocal = min
2361  (
2362  max(nChangedTotal / Pstream::nProcs(), 1),
2363  min
2364  (
2365  nChanged,
2366  max(nReportMax / Pstream::nProcs(), 1)
2367  )
2368  );
2369  }
2370 
2371  if (nReportLocal)
2372  {
2373  Pout<< "Checked mesh with layers. Disabled extrusion at " << endl;
2374  for (label i=0; i < nReportLocal; i++)
2375  {
2376  Pout<< " " << disabledFaceCentres[i] << endl;
2377  }
2378  }
2379 
2380  label nReportTotal = returnReduce(nReportLocal, sumOp<label>());
2381 
2382  if (nReportTotal < nChangedTotal)
2383  {
2384  Info<< "Suppressed disabled extrusion message for other "
2385  << nChangedTotal - nReportTotal << " faces." << endl;
2386  }
2387  }
2388 
2389  return nChangedTotal;
2390 }
2391 
2392 
2393 Foam::label Foam::snappyLayerDriver::countExtrusion
2394 (
2395  const indirectPrimitivePatch& pp,
2396  const List<extrudeMode>& extrudeStatus
2397 )
2398 {
2399  // Count number of extruded patch faces
2400  label nExtruded = 0;
2401  {
2402  const faceList& localFaces = pp.localFaces();
2403 
2404  forAll(localFaces, i)
2405  {
2406  const face& localFace = localFaces[i];
2407 
2408  forAll(localFace, fp)
2409  {
2410  if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
2411  {
2412  nExtruded++;
2413  break;
2414  }
2415  }
2416  }
2417  }
2418 
2419  return returnReduce(nExtruded, sumOp<label>());
2420 }
2421 
2422 
2423 // Collect layer faces and layer cells into mesh fields for ease of handling
2424 void Foam::snappyLayerDriver::getLayerCellsFaces
2425 (
2426  const polyMesh& mesh,
2427  const addPatchCellLayer& addLayer,
2428  const scalarField& oldRealThickness,
2429 
2430  labelList& cellNLayers,
2431  scalarField& faceRealThickness
2432 )
2433 {
2434  cellNLayers.setSize(mesh.nCells());
2435  cellNLayers = 0;
2436  faceRealThickness.setSize(mesh.nFaces());
2437  faceRealThickness = 0;
2438 
2439  // Mark all faces in the layer
2440  const labelListList& layerFaces = addLayer.layerFaces();
2441 
2442  // Mark all cells in the layer.
2443  labelListList addedCells(addPatchCellLayer::addedCells(mesh, layerFaces));
2444 
2445  forAll(addedCells, oldPatchFacei)
2446  {
2447  const labelList& added = addedCells[oldPatchFacei];
2448 
2449  const labelList& layer = layerFaces[oldPatchFacei];
2450 
2451  if (layer.size())
2452  {
2453  forAll(added, i)
2454  {
2455  cellNLayers[added[i]] = layer.size()-1;
2456  }
2457  }
2458  }
2459 
2460  forAll(layerFaces, oldPatchFacei)
2461  {
2462  const labelList& layer = layerFaces[oldPatchFacei];
2463  const scalar realThickness = oldRealThickness[oldPatchFacei];
2464 
2465  if (layer.size())
2466  {
2467  // Layer contains both original boundary face and new boundary
2468  // face so is nLayers+1. Leave out old internal face.
2469  for (label i = 1; i < layer.size(); i++)
2470  {
2471  faceRealThickness[layer[i]] = realThickness;
2472  }
2473  }
2474  }
2475 }
2476 
2477 
2478 void Foam::snappyLayerDriver::printLayerData
2479 (
2480  const fvMesh& mesh,
2481  const labelList& patchIDs,
2482  const labelList& cellNLayers,
2483  const scalarField& faceWantedThickness,
2484  const scalarField& faceRealThickness
2485 ) const
2486 {
2487  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
2488 
2489  int oldPrecision = Info().precision();
2490 
2491  // Find maximum length of a patch name, for a nicer output
2492  label maxPatchNameLen = 0;
2493  forAll(patchIDs, i)
2494  {
2495  label patchi = patchIDs[i];
2496  word patchName = pbm[patchi].name();
2497  maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
2498  }
2499 
2500  Info<< nl
2501  << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
2502  << setw(0) << " faces layers overall thickness" << nl
2503  << setf(ios_base::left) << setw(maxPatchNameLen) << " "
2504  << setw(0) << " [m] [%]" << nl
2505  << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
2506  << setw(0) << " ----- ------ --- ---" << endl;
2507 
2508 
2509  forAll(patchIDs, i)
2510  {
2511  label patchi = patchIDs[i];
2512  const polyPatch& pp = pbm[patchi];
2513 
2514  label sumSize = pp.size();
2515 
2516  // Number of layers
2517  const labelList& faceCells = pp.faceCells();
2518  label sumNLayers = 0;
2519  forAll(faceCells, i)
2520  {
2521  sumNLayers += cellNLayers[faceCells[i]];
2522  }
2523 
2524  // Thickness
2525  scalarField::subField patchWanted = pbm[patchi].patchSlice
2526  (
2527  faceWantedThickness
2528  );
2529  scalarField::subField patchReal = pbm[patchi].patchSlice
2530  (
2531  faceRealThickness
2532  );
2533 
2534  scalar sumRealThickness = sum(patchReal);
2535  scalar sumFraction = 0;
2536  forAll(patchReal, i)
2537  {
2538  if (patchWanted[i] > VSMALL)
2539  {
2540  sumFraction += (patchReal[i]/patchWanted[i]);
2541  }
2542  }
2543 
2544 
2545  reduce(sumSize, sumOp<label>());
2546  reduce(sumNLayers, sumOp<label>());
2547  reduce(sumRealThickness, sumOp<scalar>());
2548  reduce(sumFraction, sumOp<scalar>());
2549 
2550 
2551  scalar avgLayers = 0;
2552  scalar avgReal = 0;
2553  scalar avgFraction = 0;
2554  if (sumSize > 0)
2555  {
2556  avgLayers = scalar(sumNLayers)/sumSize;
2557  avgReal = sumRealThickness/sumSize;
2558  avgFraction = sumFraction/sumSize;
2559  }
2560 
2561  Info<< setf(ios_base::left) << setw(maxPatchNameLen)
2562  << pbm[patchi].name() << setprecision(3)
2563  << " " << setw(8) << sumSize
2564  << " " << setw(8) << avgLayers
2565  << " " << setw(8) << avgReal
2566  << " " << setw(8) << 100*avgFraction
2567  << endl;
2568  }
2569  Info<< setprecision(oldPrecision) << endl;
2570 }
2571 
2572 
2573 bool Foam::snappyLayerDriver::writeLayerData
2574 (
2575  const fvMesh& mesh,
2576  const labelList& patchIDs,
2577  const labelList& cellNLayers,
2578  const scalarField& faceWantedThickness,
2579  const scalarField& faceRealThickness
2580 ) const
2581 {
2582  bool allOk = true;
2583 
2585  {
2586  {
2587  label nAdded = 0;
2588  forAll(cellNLayers, celli)
2589  {
2590  if (cellNLayers[celli] > 0)
2591  {
2592  nAdded++;
2593  }
2594  }
2595  cellSet addedCellSet(mesh, "addedCells", nAdded);
2596  forAll(cellNLayers, celli)
2597  {
2598  if (cellNLayers[celli] > 0)
2599  {
2600  addedCellSet.insert(celli);
2601  }
2602  }
2603  addedCellSet.instance() = meshRefiner_.timeName();
2604  Info<< "Writing "
2605  << returnReduce(addedCellSet.size(), sumOp<label>())
2606  << " added cells to cellSet "
2607  << addedCellSet.name() << endl;
2608  bool ok = addedCellSet.write();
2609  allOk = allOk & ok;
2610  }
2611  {
2612  label nAdded = 0;
2613  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
2614  {
2615  if (faceRealThickness[facei] > 0)
2616  {
2617  nAdded++;
2618  }
2619  }
2620 
2621  faceSet layerFacesSet(mesh, "layerFaces", nAdded);
2622  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
2623  {
2624  if (faceRealThickness[facei] > 0)
2625  {
2626  layerFacesSet.insert(facei);
2627  }
2628  }
2629  layerFacesSet.instance() = meshRefiner_.timeName();
2630  Info<< "Writing "
2631  << returnReduce(layerFacesSet.size(), sumOp<label>())
2632  << " faces inside added layer to faceSet "
2633  << layerFacesSet.name() << endl;
2634  bool ok = layerFacesSet.write();
2635  allOk = allOk & ok;
2636  }
2637  }
2638 
2640  {
2641  Info<< nl << "Writing fields with layer information:" << incrIndent
2642  << endl;
2643  {
2645  (
2646  IOobject
2647  (
2648  "nSurfaceLayers",
2649  mesh.time().timeName(),
2650  mesh,
2653  false
2654  ),
2655  mesh,
2656  dimensionedScalar("zero", dimless, 0),
2657  fixedValueFvPatchScalarField::typeName
2658  );
2659  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
2660 
2661  volScalarField::Boundary& fldBf =
2662  fld.boundaryFieldRef();
2663 
2664  forAll(patchIDs, i)
2665  {
2666  label patchi = patchIDs[i];
2667  const polyPatch& pp = pbm[patchi];
2668  const labelList& faceCells = pp.faceCells();
2669  scalarField pfld(faceCells.size());
2670  forAll(faceCells, i)
2671  {
2672  pfld[i] = cellNLayers[faceCells[i]];
2673  }
2674  fldBf[patchi] == pfld;
2675  }
2676  Info<< indent << fld.name() << " : actual number of layers"
2677  << endl;
2678  bool ok = fld.write();
2679  allOk = allOk & ok;
2680  }
2681  {
2683  (
2684  IOobject
2685  (
2686  "thickness",
2687  mesh.time().timeName(),
2688  mesh,
2691  false
2692  ),
2693  mesh,
2694  dimensionedScalar("zero", dimless, 0),
2695  fixedValueFvPatchScalarField::typeName
2696  );
2697 
2698  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
2699 
2700  volScalarField::Boundary& fldBf =
2701  fld.boundaryFieldRef();
2702 
2703  forAll(patchIDs, i)
2704  {
2705  label patchi = patchIDs[i];
2706  fldBf[patchi] == pbm[patchi].patchSlice
2707  (
2708  faceRealThickness
2709  );
2710  }
2711  Info<< indent << fld.name() << " : overall layer thickness"
2712  << endl;
2713  bool ok = fld.write();
2714  allOk = allOk & ok;
2715  }
2716  {
2718  (
2719  IOobject
2720  (
2721  "thicknessFraction",
2722  mesh.time().timeName(),
2723  mesh,
2726  false
2727  ),
2728  mesh,
2729  dimensionedScalar("zero", dimless, 0),
2730  fixedValueFvPatchScalarField::typeName
2731  );
2732 
2733  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
2734 
2735  volScalarField::Boundary& fldBf =
2736  fld.boundaryFieldRef();
2737 
2738  forAll(patchIDs, i)
2739  {
2740  label patchi = patchIDs[i];
2741 
2742  scalarField::subField patchWanted = pbm[patchi].patchSlice
2743  (
2744  faceWantedThickness
2745  );
2746  scalarField::subField patchReal = pbm[patchi].patchSlice
2747  (
2748  faceRealThickness
2749  );
2750 
2751  // Convert patchReal to relavtive thickness
2752  scalarField pfld(patchReal.size(), 0.0);
2753  forAll(patchReal, i)
2754  {
2755  if (patchWanted[i] > VSMALL)
2756  {
2757  pfld[i] = patchReal[i]/patchWanted[i];
2758  }
2759  }
2760 
2761  fldBf[patchi] == pfld;
2762  }
2763  Info<< indent << fld.name()
2764  << " : overall layer thickness (fraction"
2765  << " of desired thickness)" << endl;
2766  bool ok = fld.write();
2767  allOk = allOk & ok;
2768  }
2769  Info<< decrIndent<< endl;
2770  }
2771 
2772  //if (meshRefinement::outputLevel() & meshRefinement::OUTPUTLAYERINFO)
2773  {
2774  printLayerData
2775  (
2776  mesh,
2777  patchIDs,
2778  cellNLayers,
2779  faceWantedThickness,
2780  faceRealThickness
2781  );
2782  }
2783 
2784  return allOk;
2785 }
2786 
2787 
2788 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2789 
2790 Foam::snappyLayerDriver::snappyLayerDriver
2792  meshRefinement& meshRefiner,
2793  const labelList& globalToMasterPatch,
2794  const labelList& globalToSlavePatch
2795 )
2796 :
2797  meshRefiner_(meshRefiner),
2798  globalToMasterPatch_(globalToMasterPatch),
2799  globalToSlavePatch_(globalToSlavePatch)
2800 {}
2801 
2802 
2803 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2804 
2807  const layerParameters& layerParams,
2808  const dictionary& motionDict
2809 )
2810 {
2811  // Clip to 30 degrees. Not helpful!
2812  //scalar planarAngle = min(30.0, layerParams.featureAngle());
2813  scalar planarAngle = layerParams.featureAngle();
2814  scalar minCos = Foam::cos(degToRad(planarAngle));
2815 
2816  scalar concaveCos = Foam::cos(degToRad(layerParams.concaveAngle()));
2817 
2818  Info<< nl
2819  << "Merging all faces of a cell" << nl
2820  << "---------------------------" << nl
2821  << " - which are on the same patch" << nl
2822  << " - which make an angle < " << planarAngle
2823  << " degrees"
2824  << nl
2825  << " (cos:" << minCos << ')' << nl
2826  << " - as long as the resulting face doesn't become concave"
2827  << " by more than "
2828  << layerParams.concaveAngle() << " degrees" << nl
2829  << " (0=straight, 180=fully concave)" << nl
2830  << endl;
2831 
2832  const fvMesh& mesh = meshRefiner_.mesh();
2833 
2835 
2836  labelList duplicateFace(mesh.nFaces(), -1);
2837  forAll(couples, i)
2838  {
2839  const labelPair& cpl = couples[i];
2840  duplicateFace[cpl[0]] = cpl[1];
2841  duplicateFace[cpl[1]] = cpl[0];
2842  }
2843 
2844  label nChanged = meshRefiner_.mergePatchFacesUndo
2845  (
2846  minCos,
2847  concaveCos,
2848  meshRefiner_.meshedPatches(),
2849  motionDict,
2850  duplicateFace
2851  );
2852 
2853  nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
2854 }
2855 
2856 
2859  const layerParameters& layerParams,
2860  const dictionary& motionDict,
2861  const labelList& patchIDs,
2862  const label nAllowableErrors,
2863  decompositionMethod& decomposer,
2864  fvMeshDistribute& distributor
2865 )
2866 {
2867  fvMesh& mesh = meshRefiner_.mesh();
2868 
2869  // Create baffles (pairs of faces that share the same points)
2870  // Baffles stored as owner and neighbour face that have been created.
2871  List<labelPair> baffles;
2872  meshRefiner_.createZoneBaffles
2873  (
2874  globalToMasterPatch_,
2875  globalToSlavePatch_,
2876  baffles
2877  );
2878 
2879  if (debug&meshRefinement::MESH)
2880  {
2881  const_cast<Time&>(mesh.time())++;
2882  Info<< "Writing baffled mesh to time "
2883  << meshRefiner_.timeName() << endl;
2884  meshRefiner_.write
2885  (
2888  (
2891  ),
2892  mesh.time().path()/meshRefiner_.timeName()
2893  );
2894  }
2895 
2896 
2898  (
2900  (
2901  mesh,
2902  patchIDs
2903  )
2904  );
2905 
2906 
2907  // Global face indices engine
2908  const globalIndex globalFaces(mesh.nFaces());
2909 
2910  // Determine extrudePatch.edgeFaces in global numbering (so across
2911  // coupled patches). This is used only to string up edges between coupled
2912  // faces (all edges between same (global)face indices get extruded).
2913  labelListList edgeGlobalFaces
2914  (
2916  (
2917  mesh,
2918  globalFaces,
2919  pp
2920  )
2921  );
2922 
2923  // Determine patches for extruded boundary edges. Calculates if any
2924  // additional processor patches need to be constructed.
2925 
2926  labelList sidePatchID;
2927  determineSidePatches
2928  (
2929  globalFaces,
2930  edgeGlobalFaces,
2931  pp,
2932 
2933  sidePatchID
2934  );
2935 
2936 
2937  // Point-wise extrusion data
2938  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2939 
2940  // Displacement for all pp.localPoints.
2941  vectorField patchDisp(pp().nPoints(), vector::one);
2942 
2943  // Number of layers for all pp.localPoints. Note: only valid if
2944  // extrudeStatus = EXTRUDE.
2945  labelList patchNLayers(pp().nPoints(), 0);
2946 
2947  // Ideal number of cells added
2948  label nIdealTotAddedCells = 0;
2949 
2950  // Whether to add edge for all pp.localPoints.
2951  List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
2952 
2953 
2954  {
2955  // Get number of layer per point from number of layers per patch
2956  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2957 
2958  setNumLayers
2959  (
2960  layerParams.numLayers(), // per patch the num layers
2961  patchIDs, // patches that are being moved
2962  pp, // indirectpatch for all faces moving
2963 
2964  patchDisp,
2965  patchNLayers,
2966  extrudeStatus,
2967  nIdealTotAddedCells
2968  );
2969 
2970  // Precalculate mesh edge labels for patch edges
2971  labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
2972 
2973  // Disable extrusion on non-manifold points
2974  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2975 
2976  handleNonManifolds
2977  (
2978  pp,
2979  meshEdges,
2980  edgeGlobalFaces,
2981 
2982  patchDisp,
2983  patchNLayers,
2984  extrudeStatus
2985  );
2986 
2987  // Disable extrusion on feature angles
2988  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2989 
2990  handleFeatureAngle
2991  (
2992  pp,
2993  meshEdges,
2994  degToRad(layerParams.featureAngle()),
2995 
2996  patchDisp,
2997  patchNLayers,
2998  extrudeStatus
2999  );
3000 
3001  // Disable extrusion on warped faces
3002  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3003  // It is hard to calculate some length scale if not in relative
3004  // mode so disable this check.
3005  if (layerParams.relativeSizes())
3006  {
3007  // Undistorted edge length
3008  const scalar edge0Len =
3009  meshRefiner_.meshCutter().level0EdgeLength();
3010  const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
3011 
3012  handleWarpedFaces
3013  (
3014  pp,
3015  layerParams.maxFaceThicknessRatio(),
3016  edge0Len,
3017  cellLevel,
3018 
3019  patchDisp,
3020  patchNLayers,
3021  extrudeStatus
3022  );
3023  }
3024 
3027  //
3028  //handleMultiplePatchFaces
3029  //(
3030  // pp,
3031  //
3032  // patchDisp,
3033  // patchNLayers,
3034  // extrudeStatus
3035  //);
3036 
3037 
3038  // Grow out region of non-extrusion
3039  for (label i = 0; i < layerParams.nGrow(); i++)
3040  {
3041  growNoExtrusion
3042  (
3043  pp,
3044  patchDisp,
3045  patchNLayers,
3046  extrudeStatus
3047  );
3048  }
3049  }
3050 
3051 
3052  // Undistorted edge length
3053  const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
3054  const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
3055 
3056  // Determine (wanted) point-wise overall layer thickness and expansion
3057  // ratio
3058  scalarField thickness(pp().nPoints());
3059  scalarIOField minThickness
3060  (
3061  IOobject
3062  (
3063  "minThickness",
3064  meshRefiner_.timeName(),
3065  mesh,
3067  ),
3068  pp().nPoints()
3069  );
3070  scalarField expansionRatio(pp().nPoints());
3071  calculateLayerThickness
3072  (
3073  pp,
3074  patchIDs,
3075  layerParams,
3076  cellLevel,
3077  patchNLayers,
3078  edge0Len,
3079 
3080  thickness,
3081  minThickness,
3082  expansionRatio
3083  );
3084 
3085 
3086 
3087  // Overall displacement field
3088  pointVectorField displacement
3089  (
3090  makeLayerDisplacementField
3091  (
3092  pointMesh::New(mesh),
3093  layerParams.numLayers()
3094  )
3095  );
3096 
3097  // Allocate run-time selectable mesh mover
3098  autoPtr<externalDisplacementMeshMover> medialAxisMoverPtr;
3099  {
3100  // Set up controls for meshMover
3101  dictionary combinedDict(layerParams.dict());
3102  // Add mesh quality constraints
3103  combinedDict.merge(motionDict);
3104  // Where to get minThickness from
3105  combinedDict.add("minThicknessName", minThickness.name());
3106 
3107  // Take over patchDisp as boundary conditions on displacement
3108  // pointVectorField
3109  medialAxisMoverPtr = externalDisplacementMeshMover::New
3110  (
3111  layerParams.meshShrinker(),
3112  combinedDict,
3113  baffles,
3114  displacement
3115  );
3116  }
3117 
3118 
3119  // Saved old points
3120  pointField oldPoints(mesh.points());
3121 
3122  // Current set of topology changes. (changing mesh clears out
3123  // polyTopoChange)
3124  polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
3125  // Per cell 0 or number of layers in the cell column it is part of
3126  labelList cellNLayers;
3127  // Per face actual overall layer thickness
3128  scalarField faceRealThickness;
3129  // Per face wanted overall layer thickness
3130  scalarField faceWantedThickness(mesh.nFaces(), 0.0);
3131  {
3132  UIndirectList<scalar>(faceWantedThickness, pp().addressing()) =
3133  avgPointData(pp, thickness);
3134  }
3135 
3136 
3137  for (label iteration = 0; iteration < layerParams.nLayerIter(); iteration++)
3138  {
3139  Info<< nl
3140  << "Layer addition iteration " << iteration << nl
3141  << "--------------------------" << endl;
3142 
3143 
3144  // Unset the extrusion at the pp.
3145  const dictionary& meshQualityDict =
3146  (
3147  iteration < layerParams.nRelaxedIter()
3148  ? motionDict
3149  : motionDict.subDict("relaxed")
3150  );
3151 
3152  if (iteration >= layerParams.nRelaxedIter())
3153  {
3154  Info<< "Switched to relaxed meshQuality constraints." << endl;
3155  }
3156 
3157 
3158 
3159  // Make sure displacement is equal on both sides of coupled patches.
3160  syncPatchDisplacement
3161  (
3162  pp,
3163  minThickness,
3164  patchDisp,
3165  patchNLayers,
3166  extrudeStatus
3167  );
3168 
3169  // Displacement acc. to pointnormals
3170  getPatchDisplacement
3171  (
3172  pp,
3173  thickness,
3174  minThickness,
3175  patchDisp,
3176  patchNLayers,
3177  extrudeStatus
3178  );
3179 
3180  // Shrink mesh by displacement value first.
3181  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3182 
3183  {
3184  pointField oldPatchPos(pp().localPoints());
3185 
3186  // Take over patchDisp into pointDisplacement field and
3187  // adjust both for multi-patch constraints
3189  (
3190  patchIDs,
3191  pp,
3192  patchDisp,
3193  displacement
3194  );
3195 
3196 
3197  // Move mesh
3198  // ~~~~~~~~~
3199 
3200  // Set up controls for meshMover
3201  dictionary combinedDict(layerParams.dict());
3202  // Add standard quality constraints
3203  combinedDict.merge(motionDict);
3204  // Add relaxed constraints (overrides standard ones)
3205  combinedDict.merge(meshQualityDict);
3206  // Where to get minThickness from
3207  combinedDict.add("minThicknessName", minThickness.name());
3208 
3209  labelList checkFaces(identity(mesh.nFaces()));
3210  medialAxisMoverPtr().move
3211  (
3212  combinedDict,
3213  nAllowableErrors,
3214  checkFaces
3215  );
3216 
3217  pp().movePoints(mesh.points());
3218 
3219  // Update patchDisp (since not all might have been honoured)
3220  patchDisp = oldPatchPos - pp().localPoints();
3221  }
3222 
3223  // Truncate displacements that are too small (this will do internal
3224  // ones, coupled ones have already been truncated by
3225  // syncPatchDisplacement)
3226  faceSet dummySet(mesh, "wrongPatchFaces", 0);
3227  truncateDisplacement
3228  (
3229  globalFaces,
3230  edgeGlobalFaces,
3231  pp,
3232  minThickness,
3233  dummySet,
3234  patchDisp,
3235  patchNLayers,
3236  extrudeStatus
3237  );
3238 
3239 
3240  // Dump to .obj file for debugging.
3241  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
3242  {
3243  dumpDisplacement
3244  (
3245  mesh.time().path()/"layer_" + meshRefiner_.timeName(),
3246  pp(),
3247  patchDisp,
3248  extrudeStatus
3249  );
3250 
3251  const_cast<Time&>(mesh.time())++;
3252  Info<< "Writing shrunk mesh to time "
3253  << meshRefiner_.timeName() << endl;
3254 
3255  // See comment in snappySnapDriver why we should not remove meshPhi
3256  // using mesh.clearOut().
3257 
3258  meshRefiner_.write
3259  (
3262  (
3265  ),
3266  mesh.time().path()/meshRefiner_.timeName()
3267  );
3268  }
3269 
3270 
3271  // Mesh topo change engine
3272  polyTopoChange meshMod(mesh);
3273 
3274  // Grow layer of cells on to patch. Handles zero sized displacement.
3275  addPatchCellLayer addLayer(mesh);
3276 
3277  // Determine per point/per face number of layers to extrude. Also
3278  // handles the slow termination of layers when going switching layers
3279 
3280  labelList nPatchPointLayers(pp().nPoints(), -1);
3281  labelList nPatchFaceLayers(pp().size(), -1);
3282  setupLayerInfoTruncation
3283  (
3284  pp,
3285  patchNLayers,
3286  extrudeStatus,
3287  layerParams.nBufferCellsNoExtrude(),
3288  nPatchPointLayers,
3289  nPatchFaceLayers
3290  );
3291 
3292  // Calculate displacement for final layer for addPatchLayer.
3293  // (layer of cells next to the original mesh)
3294  vectorField finalDisp(patchNLayers.size(), Zero);
3295 
3296  forAll(nPatchPointLayers, i)
3297  {
3298  scalar ratio = layerParams.finalLayerThicknessRatio
3299  (
3300  nPatchPointLayers[i],
3301  expansionRatio[i]
3302  );
3303  finalDisp[i] = ratio*patchDisp[i];
3304  }
3305 
3306 
3307  const scalarField invExpansionRatio(1.0/expansionRatio);
3308 
3309  // Add topo regardless of whether extrudeStatus is extruderemove.
3310  // Not add layer if patchDisp is zero.
3311  addLayer.setRefinement
3312  (
3313  globalFaces,
3314  edgeGlobalFaces,
3315 
3316  invExpansionRatio,
3317  pp(),
3318  sidePatchID, // boundary patch for extruded boundary edges
3319  labelList(0), // exposed patchIDs, not used for adding layers
3320  nPatchFaceLayers, // layers per face
3321  nPatchPointLayers, // layers per point
3322  finalDisp, // thickness of layer nearest internal mesh
3323  meshMod
3324  );
3325 
3326  if (debug)
3327  {
3328  const_cast<Time&>(mesh.time())++;
3329  }
3330 
3331  // Store mesh changes for if mesh is correct.
3332  savedMeshMod = meshMod;
3333 
3334 
3335  // With the stored topo changes we create a new mesh so we can
3336  // undo if neccesary.
3337 
3338  autoPtr<fvMesh> newMeshPtr;
3339  autoPtr<mapPolyMesh> map = meshMod.makeMesh
3340  (
3341  newMeshPtr,
3342  IOobject
3343  (
3344  //mesh.name()+"_layer",
3345  mesh.name(),
3346  static_cast<polyMesh&>(mesh).instance(),
3347  mesh.time(), // register with runTime
3349  static_cast<polyMesh&>(mesh).writeOpt()
3350  ), // io params from original mesh but new name
3351  mesh, // original mesh
3352  true // parallel sync
3353  );
3354  fvMesh& newMesh = newMeshPtr();
3355 
3356  //?neccesary? Update fields
3357  newMesh.updateMesh(map);
3358 
3359  newMesh.setInstance(meshRefiner_.timeName());
3360 
3361  // Update numbering on addLayer:
3362  // - cell/point labels to be newMesh.
3363  // - patchFaces to remain in oldMesh order.
3364  addLayer.updateMesh
3365  (
3366  map,
3367  identity(pp().size()),
3368  identity(pp().nPoints())
3369  );
3370 
3371  // Update numbering of baffles
3372  List<labelPair> newMeshBaffles(baffles.size());
3373  forAll(baffles, i)
3374  {
3375  const labelPair& p = baffles[i];
3376  newMeshBaffles[i][0] = map().reverseFaceMap()[p[0]];
3377  newMeshBaffles[i][1] = map().reverseFaceMap()[p[1]];
3378  }
3379 
3380  // Collect layer faces and cells for outside loop.
3381  getLayerCellsFaces
3382  (
3383  newMesh,
3384  addLayer,
3385  avgPointData(pp, mag(patchDisp))(), // current thickness
3386 
3387  cellNLayers,
3388  faceRealThickness
3389  );
3390 
3391 
3392  // Count number of added cells
3393  label nAddedCells = 0;
3394  forAll(cellNLayers, celli)
3395  {
3396  if (cellNLayers[celli] > 0)
3397  {
3398  nAddedCells++;
3399  }
3400  }
3401 
3402 
3403  if (debug&meshRefinement::MESH)
3404  {
3405  Info<< "Writing layer mesh to time " << meshRefiner_.timeName()
3406  << endl;
3407  newMesh.write();
3408 
3409  cellSet addedCellSet(newMesh, "addedCells", nAddedCells);
3410  forAll(cellNLayers, celli)
3411  {
3412  if (cellNLayers[celli] > 0)
3413  {
3414  addedCellSet.insert(celli);
3415  }
3416  }
3417  addedCellSet.instance() = meshRefiner_.timeName();
3418  Info<< "Writing "
3419  << returnReduce(addedCellSet.size(), sumOp<label>())
3420  << " added cells to cellSet " << addedCellSet.name() << endl;
3421  addedCellSet.write();
3422 
3423  faceSet layerFacesSet(newMesh, "layerFaces", newMesh.nFaces()/100);
3424  for (label facei = 0; facei < newMesh.nInternalFaces(); facei++)
3425  {
3426  if (faceRealThickness[facei] > 0)
3427  {
3428  layerFacesSet.insert(facei);
3429  }
3430  }
3431  layerFacesSet.instance() = meshRefiner_.timeName();
3432  Info<< "Writing "
3433  << returnReduce(layerFacesSet.size(), sumOp<label>())
3434  << " faces inside added layer to faceSet "
3435  << layerFacesSet.name() << endl;
3436  layerFacesSet.write();
3437  }
3438 
3439 
3440  label nTotChanged = checkAndUnmark
3441  (
3442  addLayer,
3443  meshQualityDict,
3444  layerParams.additionalReporting(),
3445  newMeshBaffles,
3446  pp(),
3447  newMesh,
3448 
3449  patchDisp,
3450  patchNLayers,
3451  extrudeStatus
3452  );
3453 
3454  label nTotExtruded = countExtrusion(pp, extrudeStatus);
3455  label nTotFaces = returnReduce(pp().size(), sumOp<label>());
3456  label nTotAddedCells = returnReduce(nAddedCells, sumOp<label>());
3457 
3458  Info<< "Extruding " << nTotExtruded
3459  << " out of " << nTotFaces
3460  << " faces (" << 100.0*nTotExtruded/nTotFaces << "%)."
3461  << " Removed extrusion at " << nTotChanged << " faces."
3462  << endl
3463  << "Added " << nTotAddedCells << " out of " << nIdealTotAddedCells
3464  << " cells (" << 100.0*nTotAddedCells/nIdealTotAddedCells << "%)."
3465  << endl;
3466 
3467  if (nTotChanged == 0)
3468  {
3469  break;
3470  }
3471 
3472  // Reset mesh points and start again
3473  mesh.movePoints(oldPoints);
3474  pp().movePoints(mesh.points());
3475 
3476  // Grow out region of non-extrusion
3477  for (label i = 0; i < layerParams.nGrow(); i++)
3478  {
3479  growNoExtrusion
3480  (
3481  pp,
3482  patchDisp,
3483  patchNLayers,
3484  extrudeStatus
3485  );
3486  }
3487 
3488  Info<< endl;
3489  }
3490 
3491 
3492  // At this point we have a (shrunk) mesh and a set of topology changes
3493  // which will make a valid mesh with layer. Apply these changes to the
3494  // current mesh.
3495 
3496  // Apply the stored topo changes to the current mesh.
3497  autoPtr<mapPolyMesh> map = savedMeshMod.changeMesh(mesh, false);
3498 
3499  // Hack to remove meshPhi - mapped incorrectly. TBD.
3500  mesh.clearOut();
3501 
3502  // Update fields
3503  mesh.updateMesh(map);
3504 
3505  // Move mesh (since morphing does not do this)
3506  if (map().hasMotionPoints())
3507  {
3508  mesh.movePoints(map().preMotionPoints());
3509  }
3510  else
3511  {
3512  // Delete mesh volumes.
3513  mesh.clearOut();
3514  }
3515 
3516  // Reset the instance for if in overwrite mode
3517  mesh.setInstance(meshRefiner_.timeName());
3518 
3519  meshRefiner_.updateMesh(map, labelList(0));
3520 
3521  // Update numbering of faceWantedThickness
3522  meshRefinement::updateList(map().faceMap(), scalar(0), faceWantedThickness);
3523 
3524  // Update numbering on baffles
3525  forAll(baffles, i)
3526  {
3527  labelPair& p = baffles[i];
3528  p[0] = map().reverseFaceMap()[p[0]];
3529  p[1] = map().reverseFaceMap()[p[1]];
3530  }
3531 
3532 
3533  label nBaffles = returnReduce(baffles.size(), sumOp<label>());
3534  if (nBaffles > 0)
3535  {
3536  // Merge any baffles
3537  Info<< "Converting " << nBaffles
3538  << " baffles back into zoned faces ..."
3539  << endl;
3540 
3541  autoPtr<mapPolyMesh> map = meshRefiner_.mergeBaffles(baffles);
3542 
3543  inplaceReorder(map().reverseCellMap(), cellNLayers);
3544  inplaceReorder(map().reverseFaceMap(), faceWantedThickness);
3545  inplaceReorder(map().reverseFaceMap(), faceRealThickness);
3546 
3547  Info<< "Converted baffles in = "
3548  << meshRefiner_.mesh().time().cpuTimeIncrement()
3549  << " s\n" << nl << endl;
3550  }
3551 
3552 
3553  // Do final balancing
3554  // ~~~~~~~~~~~~~~~~~~
3555 
3556  if (Pstream::parRun())
3557  {
3558  Info<< nl
3559  << "Doing final balancing" << nl
3560  << "---------------------" << nl
3561  << endl;
3562 
3563  if (debug)
3564  {
3565  const_cast<Time&>(mesh.time())++;
3566  }
3567 
3568  // Balance. No restriction on face zones and baffles.
3569  autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
3570  (
3571  false,
3572  false,
3573  scalarField(mesh.nCells(), 1.0),
3574  decomposer,
3575  distributor
3576  );
3577 
3578  // Re-distribute flag of layer faces and cells
3579  map().distributeCellData(cellNLayers);
3580  map().distributeFaceData(faceWantedThickness);
3581  map().distributeFaceData(faceRealThickness);
3582  }
3583 
3584 
3585  // Write mesh data
3586  // ~~~~~~~~~~~~~~~
3587 
3588  writeLayerData
3589  (
3590  mesh,
3591  patchIDs,
3592  cellNLayers,
3593  faceWantedThickness,
3594  faceRealThickness
3595  );
3596 }
3597 
3598 
3601  const dictionary& shrinkDict,
3602  const dictionary& motionDict,
3603  const layerParameters& layerParams,
3604  const bool preBalance,
3605  decompositionMethod& decomposer,
3606  fvMeshDistribute& distributor
3607 )
3608 {
3609  const fvMesh& mesh = meshRefiner_.mesh();
3610 
3611  Info<< nl
3612  << "Shrinking and layer addition phase" << nl
3613  << "----------------------------------" << nl
3614  << endl;
3615 
3616  Info<< "Using mesh parameters " << motionDict << nl << endl;
3617 
3618  // Merge coplanar boundary faces
3619  mergePatchFacesUndo(layerParams, motionDict);
3620 
3621  // Per patch the number of layers (-1 or 0 if no layer)
3622  const labelList& numLayers = layerParams.numLayers();
3623 
3624  // Patches that need to get a layer
3625  DynamicList<label> patchIDs(numLayers.size());
3626  label nFacesWithLayers = 0;
3627  forAll(numLayers, patchi)
3628  {
3629  if (numLayers[patchi] > 0)
3630  {
3631  const polyPatch& pp = mesh.boundaryMesh()[patchi];
3632 
3633  if (!pp.coupled())
3634  {
3635  patchIDs.append(patchi);
3636  nFacesWithLayers += mesh.boundaryMesh()[patchi].size();
3637  }
3638  else
3639  {
3641  << "Ignoring layers on coupled patch " << pp.name()
3642  << endl;
3643  }
3644  }
3645  }
3646  patchIDs.shrink();
3647 
3648  if (returnReduce(nFacesWithLayers, sumOp<label>()) == 0)
3649  {
3650  Info<< nl << "No layers to generate ..." << endl;
3651  }
3652  else
3653  {
3654  // Check that outside of mesh is not multiply connected.
3655  checkMeshManifold();
3656 
3657  // Check initial mesh
3658  Info<< "Checking initial mesh ..." << endl;
3659  labelHashSet wrongFaces(mesh.nFaces()/100);
3660  motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
3661  const label nInitErrors = returnReduce
3662  (
3663  wrongFaces.size(),
3664  sumOp<label>()
3665  );
3666 
3667  Info<< "Detected " << nInitErrors << " illegal faces"
3668  << " (concave, zero area or negative cell pyramid volume)"
3669  << endl;
3670 
3671 
3672  // Balance
3673  if (Pstream::parRun() && preBalance)
3674  {
3675  Info<< nl
3676  << "Doing initial balancing" << nl
3677  << "-----------------------" << nl
3678  << endl;
3679 
3680  scalarField cellWeights(mesh.nCells(), 1);
3681  forAll(numLayers, patchi)
3682  {
3683  if (numLayers[patchi] > 0)
3684  {
3685  const polyPatch& pp = mesh.boundaryMesh()[patchi];
3686  forAll(pp.faceCells(), i)
3687  {
3688  cellWeights[pp.faceCells()[i]] += numLayers[patchi];
3689  }
3690  }
3691  }
3692 
3693  // Balance mesh (and meshRefinement). Restrict faceZones to
3694  // be on internal faces only since they will be converted into
3695  // baffles.
3696  autoPtr<mapDistributePolyMesh> map = meshRefiner_.balance
3697  (
3698  true, //false, // keepZoneFaces
3699  false,
3700  cellWeights,
3701  decomposer,
3702  distributor
3703  );
3704  }
3705 
3706 
3707  // Do all topo changes
3708  addLayers
3709  (
3710  layerParams,
3711  motionDict,
3712  patchIDs,
3713  nInitErrors,
3714  decomposer,
3715  distributor
3716  );
3717  }
3718 }
3719 
3720 
3721 // ************************************************************************* //
label nPatches
Definition: readKivaGrid.H:402
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void clearOut()
Clear all geometry and addressing.
Definition: fvMesh.C:230
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
const labelList & numLayers() const
How many layers to add:
scalar featureAngle() const
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
Simple container to keep together layer specific information.
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
label nLayerIter() const
Number of overall layer addition iterations.
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:223
A list of face labels.
Definition: faceSet.H:48
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
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.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
const dictionary & dict() const
Unit conversion functions.
static autoPtr< externalDisplacementMeshMover > New(const word &type, const dictionary &dict, const List< labelPair > &baffles, pointVectorField &pointDisplacement)
Return a reference to the selected meshMover model.
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 dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:633
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:417
List< face > faceList
Definition: faceListFwd.H:43
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
static word newName(const label myProcNo, const label neighbProcNo)
Return the name of a processorPolyPatch.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static writeType writeLevel()
Get/set write level.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
bool relativeSizes() const
Are size parameters relative to inner cell size or.
HashTable< label, label, Hash< label > >::const_iterator const_iterator
Definition: Map.H:59
void addLayers(const layerParameters &layerParams, const dictionary &motionDict, const labelList &patchIDs, const label nAllowableErrors, decompositionMethod &decomposer, fvMeshDistribute &distributor)
Add cell layers.
label nRelaxedIter() const
Number of iterations after which relaxed motion rules.
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
label nCells() const
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
static const pointMesh & New(const polyMesh &mesh)
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
static void updateList(const labelList &newToOld, const T &nullValue, List< T > &elems)
Helper: reorder list according to map.
Omanip< int > setprecision(const int i)
Definition: IOmanip.H:205
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
virtual void updateMesh(const mapPolyMesh &mpm)
Update mesh corresponding to the given map.
Definition: fvMesh.C:800
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
dynamicFvMesh & mesh
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< ' ';}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< ' ';}gmvFile<< nl;forAll(lagrangianScalarNames, i){const word &name=lagrangianScalarNames[i];IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
scalar maxFaceThicknessRatio() const
Stop layer growth on highly warped cells.
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
Adds layers of cells to outside of polyPatch. Can optionally create stand-alone extruded mesh (addToM...
void doLayers(const dictionary &shrinkDict, const dictionary &motionDict, const layerParameters &layerParams, const bool preBalance, decompositionMethod &decomposer, fvMeshDistribute &distributor)
Add layers according to the dictionary settings.
label nPoints
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const Switch & additionalReporting() const
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
label nGrow() const
If points get not extruded do nGrow layers of connected faces.
Smanip< ios_base::fmtflags > setf(const ios_base::fmtflags flags)
Definition: IOmanip.H:164
static const label labelMax
Definition: label.H:62
List< label > labelList
A List of labels.
Definition: labelList.H:56
const labelListList & pointEdges() const
static const zero Zero
Definition: zero.H:91
Abstract base class for decomposition.
const word & name() const
Return name.
labelListList addedCells() const
Added cells given current mesh & layerfaces.
bool visNormal(const vector &n, const vectorField &faceNormals, const labelList &faceLabels)
Check if n is in same direction as normals of all faceLabels.
Definition: meshTools.C:33
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static tmp< pointField > pointNormals(const polyMesh &, const PrimitivePatch< Face, FaceList, PointField, PointType > &)
Return parallel consistent point normals for patches using mesh points.
const word & meshShrinker() const
Type of mesh shrinker.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
Istream and Ostream manipulators taking arguments.
static const char nl
Definition: Ostream.H:262
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
static PackedBoolList getMasterPoints(const polyMesh &)
Get per point whether it is uncoupled or a master of a.
Definition: syncTools.C:65
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:237
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static autoPtr< indirectPrimitivePatch > makePatch(const polyMesh &, const labelList &)
Create patch from set of patches.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
scalar finalLayerThicknessRatio(const label nLayers, const scalar expansionRatio) const
Determine ratio of final layer thickness to.
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: fvMesh.C:721
void updateMesh(const mapPolyMesh &, const labelList &faceMap, const labelList &pointMap)
Update any locally stored mesh information. Gets additional.
label nFaces() const
List< word > wordList
A List of words.
Definition: fileName.H:54
void setSize(const label)
Reset size of List.
Definition: List.C:295
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:393
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:399
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:315
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A collection of cell labels.
Definition: cellSet.H:48
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
A List with indirect addressing.
Definition: fvMatrix.H:106
Direct mesh changes based on v1.3 polyTopoChange syntax.
scalar concaveAngle() const
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
Definition: pointFields.H:50
virtual bool write() const
Write using setting from DB.
void mergePatchFacesUndo(const layerParameters &layerParams, const dictionary &motionDict)
Merge patch faces on same cell.
messageStream Info
SubField< scalar > subField
Declare type of subField.
Definition: Field.H:89
bool merge(const dictionary &)
Merge entries from the given dictionary.
Definition: dictionary.C:1013
dimensioned< scalar > mag(const dimensioned< Type > &)
static labelListList globalEdgeFaces(const polyMesh &, const globalIndex &globalFaces, const indirectPrimitivePatch &pp)
Per patch edge the pp faces (in global indices) using it. Uses.
Field< vector > vectorField
Specialisation of Field<T> for vector.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
virtual Ostream & write(const token &)=0
Write next token to stream.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
volScalarField & p
A class for managing temporary objects.
Definition: PtrList.H:54
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
fileName path() const
Return path.
Definition: Time.H:269
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:230
void setRefinement(const globalIndex &globalFaces, const labelListList &globalEdgeFaces, const scalarField &expansionRatio, const indirectPrimitivePatch &pp, const labelList &sidePatchID, const labelList &exposedPatchID, const labelList &nFaceLayers, const labelList &nPointLayers, const vectorField &firstLayerDisp, polyTopoChange &meshMod)
Play commands into polyTopoChange to create layers on top.
label nBufferCellsNoExtrude() const
Create buffer region for new layer terminations.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
static void calcSidePatch(const polyMesh &, const globalIndex &globalFaces, const labelListList &globalEdgeFaces, const indirectPrimitivePatch &pp, labelList &sidePatchID, label &nPatches, Map< label > &nbrProcToPatch, Map< label > &patchToNbrProc)
Boundary edges get extruded into boundary faces. Determine patch.
const fileName & instance() const
Definition: IOobject.H:337
label nInternalFaces() const
A primitive field of type <T> with automated input and output.
Definition: IOField.H:50
virtual bool write() const
Write mesh using IO settings from time.
Definition: fvMesh.C:870
const word & name() const
Return reference to name.
Definition: fvMesh.H:257
static const label labelMin
Definition: label.H:61
const word & name() const
Return name.
Definition: IOobject.H:260
Namespace for OpenFOAM.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29