snappyLayerDriver.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 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 "polyTopoChangeMap.H"
42 #include "addPatchCellLayer.H"
43 #include "polyDistributionMap.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  normalsCombine()
418  (
419  edgeNormal[meshEdgeI],
420  pp.faceNormals()[eFaces[i]]
421  );
422  }
423  }
424 
426  (
427  mesh,
428  edgeNormal,
429  normalsCombine(),
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  // Note: time().timeName() instead of meshRefinement::timeName() since
855  // postprocessable field.
856  tmp<pointVectorField> tfld
857  (
859  (
860  "pointDisplacement",
861  pMesh,
863  patchFieldTypes,
864  actualPatchTypes
865  )
866  );
867  return tfld;
868 }
869 
870 
871 void Foam::snappyLayerDriver::growNoExtrusion
872 (
873  const indirectPrimitivePatch& pp,
874  pointField& patchDisp,
875  labelList& patchNLayers,
876  List<extrudeMode>& extrudeStatus
877 ) const
878 {
879  Info<< nl << "Growing non-extrusion points by one layer ..." << endl;
880 
881  List<extrudeMode> grownExtrudeStatus(extrudeStatus);
882 
883  const faceList& localFaces = pp.localFaces();
884 
885  label nGrown = 0;
886 
887  forAll(localFaces, facei)
888  {
889  const face& f = localFaces[facei];
890 
891  bool hasSqueeze = false;
892  forAll(f, fp)
893  {
894  if (extrudeStatus[f[fp]] == NOEXTRUDE)
895  {
896  hasSqueeze = true;
897  break;
898  }
899  }
900 
901  if (hasSqueeze)
902  {
903  // Squeeze all points of face
904  forAll(f, fp)
905  {
906  if
907  (
908  extrudeStatus[f[fp]] == EXTRUDE
909  && grownExtrudeStatus[f[fp]] != NOEXTRUDE
910  )
911  {
912  grownExtrudeStatus[f[fp]] = NOEXTRUDE;
913  nGrown++;
914  }
915  }
916  }
917  }
918 
919  extrudeStatus.transfer(grownExtrudeStatus);
920 
921 
922  // Synchronise since might get called multiple times.
923  // Use the fact that NOEXTRUDE is the minimum value.
924  {
925  labelList status(extrudeStatus.size());
926  forAll(status, i)
927  {
928  status[i] = extrudeStatus[i];
929  }
931  (
932  meshRefiner_.mesh(),
933  pp.meshPoints(),
934  status,
935  minEqOp<label>(),
936  labelMax // null value
937  );
938  forAll(status, i)
939  {
940  extrudeStatus[i] = extrudeMode(status[i]);
941  }
942  }
943 
944 
945  forAll(extrudeStatus, patchPointi)
946  {
947  if (extrudeStatus[patchPointi] == NOEXTRUDE)
948  {
949  patchDisp[patchPointi] = Zero;
950  patchNLayers[patchPointi] = 0;
951  }
952  }
953 
954  reduce(nGrown, sumOp<label>());
955 
956  Info<< "Set displacement to zero for an additional " << nGrown
957  << " points." << endl;
958 }
959 
960 
961 void Foam::snappyLayerDriver::determineSidePatches
962 (
963  const globalIndex& globalFaces,
964  const labelListList& edgeGlobalFaces,
965  const indirectPrimitivePatch& pp,
966 
967  labelList& sidePatchID
968 )
969 {
970  // Sometimes edges-to-be-extruded are on more than 2 processors.
971  // Work out which 2 hold the faces to be extruded and thus which procpatch
972  // the side-face should be in. As an additional complication this might
973  // mean that 2 processors that were only edge-connected now suddenly need
974  // to become face-connected i.e. have a processor patch between them.
975 
976  fvMesh& mesh = meshRefiner_.mesh();
977 
978  // Determine sidePatchID. Any additional processor boundary gets added to
979  // patchToNbrProc,nbrProcToPatch and nPatches gets set to the new number
980  // of patches.
981  label nPatches;
982  Map<label> nbrProcToPatch;
983  Map<label> patchToNbrProc;
985  (
986  mesh,
987  globalFaces,
988  edgeGlobalFaces,
989  pp,
990 
991  sidePatchID,
992  nPatches,
993  nbrProcToPatch,
994  patchToNbrProc
995  );
996 
997  label nOldPatches = mesh.boundaryMesh().size();
998  label nAdded = returnReduce(nPatches-nOldPatches, sumOp<label>());
999  Info<< nl << "Adding in total " << nAdded/2 << " inter-processor patches to"
1000  << " handle extrusion of non-manifold processor boundaries."
1001  << endl;
1002 
1003  if (nAdded > 0)
1004  {
1005  // We might not add patches in same order as in patchToNbrProc
1006  // so prepare to renumber sidePatchID
1007  Map<label> wantedToAddedPatch;
1008 
1009  for (label patchi = nOldPatches; patchi < nPatches; patchi++)
1010  {
1011  const label nbrProci = patchToNbrProc[patchi];
1012  const label procPatchi = mesh.boundaryMesh().size();
1013  const processorPolyPatch pp
1014  (
1015  0, // size
1016  0, // start
1017  procPatchi, // index
1018  mesh.boundaryMesh(),
1020  nbrProci
1021  );
1022  mesh.addPatch
1023  (
1024  procPatchi, // new patch index
1025  pp,
1026  dictionary(), // patchField dict
1028  false
1029  );
1030  wantedToAddedPatch.insert(patchi, procPatchi);
1031  }
1032 
1033  // Renumber sidePatchID
1034  forAll(sidePatchID, i)
1035  {
1036  label patchi = sidePatchID[i];
1037  Map<label>::const_iterator fnd = wantedToAddedPatch.find(patchi);
1038  if (fnd != wantedToAddedPatch.end())
1039  {
1040  sidePatchID[i] = fnd();
1041  }
1042  }
1043 
1044  mesh.clearOut();
1045  const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()).topoChange();
1046  }
1047 }
1048 
1049 
1050 void Foam::snappyLayerDriver::calculateLayerThickness
1051 (
1052  const indirectPrimitivePatch& pp,
1053  const labelList& patchIDs,
1054  const layerParameters& layerParams,
1055  const labelList& cellLevel,
1056  const labelList& patchNLayers,
1057  const scalar edge0Len,
1058 
1059  scalarField& thickness,
1060  scalarField& minThickness,
1061  scalarField& expansionRatio
1062 ) const
1063 {
1064  const fvMesh& mesh = meshRefiner_.mesh();
1065  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1066 
1067 
1068  // Rework patch-wise layer parameters into minimum per point
1069  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1070  // Note: only layer parameters consistent with layer specification
1071  // method (see layerParameters) will be correct.
1072  scalarField firstLayerThickness(pp.nPoints(), great);
1073  scalarField finalLayerThickness(pp.nPoints(), great);
1074  scalarField totalThickness(pp.nPoints(), great);
1075  scalarField expRatio(pp.nPoints(), great);
1076 
1077  minThickness.setSize(pp.nPoints());
1078  minThickness = great;
1079 
1080  forAll(patchIDs, i)
1081  {
1082  label patchi = patchIDs[i];
1083 
1084  const labelList& meshPoints = patches[patchi].meshPoints();
1085 
1086  forAll(meshPoints, patchPointi)
1087  {
1088  label ppPointi = pp.meshPointMap()[meshPoints[patchPointi]];
1089 
1090  firstLayerThickness[ppPointi] = min
1091  (
1092  firstLayerThickness[ppPointi],
1093  layerParams.firstLayerThickness()[patchi]
1094  );
1095  finalLayerThickness[ppPointi] = min
1096  (
1097  finalLayerThickness[ppPointi],
1098  layerParams.finalLayerThickness()[patchi]
1099  );
1100  totalThickness[ppPointi] = min
1101  (
1102  totalThickness[ppPointi],
1103  layerParams.thickness()[patchi]
1104  );
1105  expRatio[ppPointi] = min
1106  (
1107  expRatio[ppPointi],
1108  layerParams.expansionRatio()[patchi]
1109  );
1110  minThickness[ppPointi] = min
1111  (
1112  minThickness[ppPointi],
1113  layerParams.minThickness()[patchi]
1114  );
1115  }
1116  }
1117 
1119  (
1120  mesh,
1121  pp.meshPoints(),
1122  firstLayerThickness,
1123  minEqOp<scalar>(),
1124  great // null value
1125  );
1127  (
1128  mesh,
1129  pp.meshPoints(),
1130  finalLayerThickness,
1131  minEqOp<scalar>(),
1132  great // null value
1133  );
1135  (
1136  mesh,
1137  pp.meshPoints(),
1138  totalThickness,
1139  minEqOp<scalar>(),
1140  great // null value
1141  );
1143  (
1144  mesh,
1145  pp.meshPoints(),
1146  expRatio,
1147  minEqOp<scalar>(),
1148  great // null value
1149  );
1151  (
1152  mesh,
1153  pp.meshPoints(),
1154  minThickness,
1155  minEqOp<scalar>(),
1156  great // null value
1157  );
1158 
1159 
1160  // Now the thicknesses are set according to the minimum of connected
1161  // patches.
1162 
1163 
1164  // Rework relative thickness into absolute
1165  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1166  // by multiplying with the internal cell size.
1167 
1168  if (layerParams.relativeSizes())
1169  {
1170  if
1171  (
1172  min(layerParams.minThickness()) < 0
1173  || max(layerParams.minThickness()) > 2
1174  )
1175  {
1177  << "Thickness should be factor of local undistorted cell size."
1178  << " Valid values are [0..2]." << nl
1179  << " minThickness:" << layerParams.minThickness()
1180  << exit(FatalError);
1181  }
1182 
1183 
1184  // Determine per point the max cell level of connected cells
1185  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1186 
1187  labelList maxPointLevel(pp.nPoints(), labelMin);
1188 
1189  forAll(pp, i)
1190  {
1191  label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
1192 
1193  const face& f = pp.localFaces()[i];
1194 
1195  forAll(f, fp)
1196  {
1197  maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
1198  }
1199  }
1200 
1202  (
1203  mesh,
1204  pp.meshPoints(),
1205  maxPointLevel,
1206  maxEqOp<label>(),
1207  labelMin // null value
1208  );
1209 
1210 
1211  forAll(maxPointLevel, pointi)
1212  {
1213  // Find undistorted edge size for this level.
1214  scalar edgeLen = edge0Len/(1<<maxPointLevel[pointi]);
1215  firstLayerThickness[pointi] *= edgeLen;
1216  finalLayerThickness[pointi] *= edgeLen;
1217  totalThickness[pointi] *= edgeLen;
1218  minThickness[pointi] *= edgeLen;
1219  }
1220  }
1221 
1222 
1223 
1224  // Rework thickness parameters into overall thickness
1225  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1226 
1227  forAll(firstLayerThickness, pointi)
1228  {
1229  thickness[pointi] = layerParams.layerThickness
1230  (
1231  patchNLayers[pointi],
1232  firstLayerThickness[pointi],
1233  finalLayerThickness[pointi],
1234  totalThickness[pointi],
1235  expRatio[pointi]
1236  );
1237 
1238  expansionRatio[pointi] = layerParams.layerExpansionRatio
1239  (
1240  patchNLayers[pointi],
1241  firstLayerThickness[pointi],
1242  finalLayerThickness[pointi],
1243  totalThickness[pointi],
1244  expRatio[pointi]
1245  );
1246  }
1247 
1248  // Info<< "calculateLayerThickness : min:" << gMin(thickness)
1249  // << " max:" << gMax(thickness) << endl;
1250 
1251  // Print a bit
1252  {
1253  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1254 
1255  int oldPrecision = Info().precision();
1256 
1257  // Find maximum length of a patch name, for a nicer output
1258  label maxPatchNameLen = 0;
1259  forAll(patchIDs, i)
1260  {
1261  label patchi = patchIDs[i];
1262  word patchName = patches[patchi].name();
1263  maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
1264  }
1265 
1266  Info<< nl
1267  << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
1268  << setw(0) << " faces layers avg thickness[m]" << nl
1269  << setf(ios_base::left) << setw(maxPatchNameLen) << " "
1270  << setw(0) << " near-wall overall" << nl
1271  << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
1272  << setw(0) << " ----- ------ --------- -------" << endl;
1273 
1274 
1275  const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
1276 
1277  forAll(patchIDs, i)
1278  {
1279  label patchi = patchIDs[i];
1280 
1281  const labelList& meshPoints = patches[patchi].meshPoints();
1282 
1283  scalar sumThickness = 0;
1284  scalar sumNearWallThickness = 0;
1285  label nMasterPoints = 0;
1286 
1287  forAll(meshPoints, patchPointi)
1288  {
1289  label meshPointi = meshPoints[patchPointi];
1290  if (isMasterPoint[meshPointi])
1291  {
1292  label ppPointi = pp.meshPointMap()[meshPointi];
1293 
1294  sumThickness += thickness[ppPointi];
1295  sumNearWallThickness += layerParams.firstLayerThickness
1296  (
1297  patchNLayers[ppPointi],
1298  firstLayerThickness[ppPointi],
1299  finalLayerThickness[ppPointi],
1300  thickness[ppPointi],
1301  expansionRatio[ppPointi]
1302  );
1303  nMasterPoints++;
1304  }
1305  }
1306 
1307  label totNPoints = returnReduce(nMasterPoints, sumOp<label>());
1308 
1309  // For empty patches, totNPoints is 0.
1310  scalar avgThickness = 0;
1311  scalar avgNearWallThickness = 0;
1312 
1313  if (totNPoints > 0)
1314  {
1315  avgThickness =
1316  returnReduce(sumThickness, sumOp<scalar>())
1317  / totNPoints;
1318  avgNearWallThickness =
1319  returnReduce(sumNearWallThickness, sumOp<scalar>())
1320  / totNPoints;
1321  }
1322 
1323  Info<< setf(ios_base::left) << setw(maxPatchNameLen)
1324  << patches[patchi].name() << setprecision(3)
1325  << " " << setw(8)
1326  << returnReduce(patches[patchi].size(), sumOp<scalar>())
1327  << " " << setw(6) << layerParams.numLayers()[patchi]
1328  << " " << setw(8) << avgNearWallThickness
1329  << " " << setw(8) << avgThickness
1330  << endl;
1331  }
1332  Info<< setprecision(oldPrecision) << endl;
1333  }
1334 }
1335 
1336 
1337 // Synchronise displacement among coupled patches.
1338 void Foam::snappyLayerDriver::syncPatchDisplacement
1339 (
1340  const indirectPrimitivePatch& pp,
1341  const scalarField& minThickness,
1342  pointField& patchDisp,
1343  labelList& patchNLayers,
1344  List<extrudeMode>& extrudeStatus
1345 ) const
1346 {
1347  const fvMesh& mesh = meshRefiner_.mesh();
1348  const labelList& meshPoints = pp.meshPoints();
1349 
1350  label nChangedTotal = 0;
1351 
1352  while (true)
1353  {
1354  label nChanged = 0;
1355 
1356  // Sync displacement (by taking min)
1358  (
1359  mesh,
1360  meshPoints,
1361  patchDisp,
1362  minMagSqrEqOp<vector>(),
1363  point::rootMax // null value
1364  );
1365 
1366  // Unmark if displacement too small
1367  forAll(patchDisp, i)
1368  {
1369  if (mag(patchDisp[i]) < minThickness[i])
1370  {
1371  if
1372  (
1373  unmarkExtrusion
1374  (
1375  i,
1376  patchDisp,
1377  patchNLayers,
1378  extrudeStatus
1379  )
1380  )
1381  {
1382  nChanged++;
1383  }
1384  }
1385  }
1386 
1387  labelList syncPatchNLayers(patchNLayers);
1388 
1390  (
1391  mesh,
1392  meshPoints,
1393  syncPatchNLayers,
1394  minEqOp<label>(),
1395  labelMax // null value
1396  );
1397 
1398  // Reset if differs
1399  // 1. take max
1400  forAll(syncPatchNLayers, i)
1401  {
1402  if (syncPatchNLayers[i] != patchNLayers[i])
1403  {
1404  if
1405  (
1406  unmarkExtrusion
1407  (
1408  i,
1409  patchDisp,
1410  patchNLayers,
1411  extrudeStatus
1412  )
1413  )
1414  {
1415  nChanged++;
1416  }
1417  }
1418  }
1419 
1421  (
1422  mesh,
1423  meshPoints,
1424  syncPatchNLayers,
1425  maxEqOp<label>(),
1426  labelMin // null value
1427  );
1428 
1429  // Reset if differs
1430  // 2. take min
1431  forAll(syncPatchNLayers, i)
1432  {
1433  if (syncPatchNLayers[i] != patchNLayers[i])
1434  {
1435  if
1436  (
1437  unmarkExtrusion
1438  (
1439  i,
1440  patchDisp,
1441  patchNLayers,
1442  extrudeStatus
1443  )
1444  )
1445  {
1446  nChanged++;
1447  }
1448  }
1449  }
1450  nChangedTotal += nChanged;
1451 
1452  if (!returnReduce(nChanged, sumOp<label>()))
1453  {
1454  break;
1455  }
1456  }
1457 
1458  // Info<< "Prevented extrusion on "
1459  // << returnReduce(nChangedTotal, sumOp<label>())
1460  // << " coupled patch points during syncPatchDisplacement." << endl;
1461 }
1462 
1463 
1464 // Calculate displacement vector for all patch points. Uses pointNormal.
1465 // Checks that displaced patch point would be visible from all centres
1466 // of the faces using it.
1467 // extrudeStatus is both input and output and gives the status of each
1468 // patch point.
1469 void Foam::snappyLayerDriver::getPatchDisplacement
1470 (
1471  const indirectPrimitivePatch& pp,
1472  const scalarField& thickness,
1473  const scalarField& minThickness,
1474  pointField& patchDisp,
1475  labelList& patchNLayers,
1476  List<extrudeMode>& extrudeStatus
1477 ) const
1478 {
1479  Info<< nl << "Determining displacement for added points"
1480  << " according to pointNormal ..." << endl;
1481 
1482  const fvMesh& mesh = meshRefiner_.mesh();
1483  const vectorField& faceNormals = pp.faceNormals();
1484  const labelListList& pointFaces = pp.pointFaces();
1485  const pointField& localPoints = pp.localPoints();
1486 
1487  // Determine pointNormal
1488  // ~~~~~~~~~~~~~~~~~~~~~
1489 
1490  pointField pointNormals(PatchTools::pointNormals(mesh, pp));
1491 
1492 
1493  // Determine local length scale on patch
1494  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1495 
1496  // Start off from same thickness everywhere (except where no extrusion)
1497  patchDisp = thickness*pointNormals;
1498 
1499 
1500  label nNoVisNormal = 0;
1501  label nExtrudeRemove = 0;
1502 
1503 
1504  // Check if no extrude possible.
1505  forAll(pointNormals, patchPointi)
1506  {
1507  label meshPointi = pp.meshPoints()[patchPointi];
1508 
1509  if (extrudeStatus[patchPointi] == NOEXTRUDE)
1510  {
1511  // Do not use unmarkExtrusion; forcibly set to zero extrusion.
1512  patchNLayers[patchPointi] = 0;
1513  patchDisp[patchPointi] = Zero;
1514  }
1515  else
1516  {
1517  // Get normal
1518  const vector& n = pointNormals[patchPointi];
1519 
1520  if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointi]))
1521  {
1522  if (debug&meshRefinement::ATTRACTION)
1523  {
1524  Pout<< "No valid normal for point " << meshPointi
1525  << ' ' << pp.points()[meshPointi]
1526  << "; setting displacement to "
1527  << patchDisp[patchPointi]
1528  << endl;
1529  }
1530 
1531  extrudeStatus[patchPointi] = EXTRUDEREMOVE;
1532  nNoVisNormal++;
1533  }
1534  }
1535  }
1536 
1537  // At illegal points make displacement average of new neighbour positions
1538  forAll(extrudeStatus, patchPointi)
1539  {
1540  if (extrudeStatus[patchPointi] == EXTRUDEREMOVE)
1541  {
1542  point avg(Zero);
1543  label nPoints = 0;
1544 
1545  const labelList& pEdges = pp.pointEdges()[patchPointi];
1546 
1547  forAll(pEdges, i)
1548  {
1549  label edgei = pEdges[i];
1550 
1551  label otherPointi = pp.edges()[edgei].otherVertex(patchPointi);
1552 
1553  if (extrudeStatus[otherPointi] != NOEXTRUDE)
1554  {
1555  avg += localPoints[otherPointi] + patchDisp[otherPointi];
1556  nPoints++;
1557  }
1558  }
1559 
1560  if (nPoints > 0)
1561  {
1562  if (debug&meshRefinement::ATTRACTION)
1563  {
1564  Pout<< "Displacement at illegal point "
1565  << localPoints[patchPointi]
1566  << " set to "
1567  << (avg / nPoints - localPoints[patchPointi])
1568  << endl;
1569  }
1570 
1571  patchDisp[patchPointi] =
1572  avg / nPoints
1573  - localPoints[patchPointi];
1574 
1575  nExtrudeRemove++;
1576  }
1577  else
1578  {
1579  // All surrounding points are not extruded. Leave patchDisp
1580  // intact.
1581  }
1582  }
1583  }
1584 
1585  Info<< "Detected " << returnReduce(nNoVisNormal, sumOp<label>())
1586  << " points with point normal pointing through faces." << nl
1587  << "Reset displacement at "
1588  << returnReduce(nExtrudeRemove, sumOp<label>())
1589  << " points to average of surrounding points." << endl;
1590 
1591  // Make sure displacement is equal on both sides of coupled patches.
1592  syncPatchDisplacement
1593  (
1594  pp,
1595  minThickness,
1596  patchDisp,
1597  patchNLayers,
1598  extrudeStatus
1599  );
1600 
1601  Info<< endl;
1602 }
1603 
1604 
1605 bool Foam::snappyLayerDriver::sameEdgeNeighbour
1606 (
1607  const labelListList& globalEdgeFaces,
1608  const label myGlobalFacei,
1609  const label nbrGlobFacei,
1610  const label edgei
1611 ) const
1612 {
1613  const labelList& eFaces = globalEdgeFaces[edgei];
1614  if (eFaces.size() == 2)
1615  {
1616  return edge(myGlobalFacei, nbrGlobFacei) == edge(eFaces[0], eFaces[1]);
1617  }
1618  else
1619  {
1620  return false;
1621  }
1622 }
1623 
1624 
1625 void Foam::snappyLayerDriver::getVertexString
1626 (
1627  const indirectPrimitivePatch& pp,
1628  const labelListList& globalEdgeFaces,
1629  const label facei,
1630  const label edgei,
1631  const label myGlobFacei,
1632  const label nbrGlobFacei,
1633  DynamicList<label>& vertices
1634 ) const
1635 {
1636  const labelList& fEdges = pp.faceEdges()[facei];
1637  label fp = findIndex(fEdges, edgei);
1638 
1639  if (fp == -1)
1640  {
1642  << "problem." << abort(FatalError);
1643  }
1644 
1645  // Search back
1646  label startFp = fp;
1647 
1648  forAll(fEdges, i)
1649  {
1650  label prevFp = fEdges.rcIndex(startFp);
1651  if
1652  (
1653  !sameEdgeNeighbour
1654  (
1655  globalEdgeFaces,
1656  myGlobFacei,
1657  nbrGlobFacei,
1658  fEdges[prevFp]
1659  )
1660  )
1661  {
1662  break;
1663  }
1664  startFp = prevFp;
1665  }
1666 
1667  label endFp = fp;
1668  forAll(fEdges, i)
1669  {
1670  label nextFp = fEdges.fcIndex(endFp);
1671  if
1672  (
1673  !sameEdgeNeighbour
1674  (
1675  globalEdgeFaces,
1676  myGlobFacei,
1677  nbrGlobFacei,
1678  fEdges[nextFp]
1679  )
1680  )
1681  {
1682  break;
1683  }
1684  endFp = nextFp;
1685  }
1686 
1687  const face& f = pp.localFaces()[facei];
1688  vertices.clear();
1689  fp = startFp;
1690  while (fp != endFp)
1691  {
1692  vertices.append(f[fp]);
1693  fp = f.fcIndex(fp);
1694  }
1695  vertices.append(f[fp]);
1696  fp = f.fcIndex(fp);
1697  vertices.append(f[fp]);
1698 }
1699 
1700 
1701 // Truncates displacement
1702 // - for all patchFaces in the faceset displacement gets set to zero
1703 // - all displacement < minThickness gets set to zero
1704 Foam::label Foam::snappyLayerDriver::truncateDisplacement
1705 (
1706  const globalIndex& globalFaces,
1707  const labelListList& edgeGlobalFaces,
1708  const indirectPrimitivePatch& pp,
1709  const scalarField& minThickness,
1710  const faceSet& illegalPatchFaces,
1711  pointField& patchDisp,
1712  labelList& patchNLayers,
1713  List<extrudeMode>& extrudeStatus
1714 ) const
1715 {
1716  const fvMesh& mesh = meshRefiner_.mesh();
1717 
1718  label nChanged = 0;
1719 
1720  const Map<label>& meshPointMap = pp.meshPointMap();
1721 
1722  forAllConstIter(faceSet, illegalPatchFaces, iter)
1723  {
1724  label facei = iter.key();
1725 
1726  if (mesh.isInternalFace(facei))
1727  {
1729  << "Faceset " << illegalPatchFaces.name()
1730  << " contains internal face " << facei << nl
1731  << "It should only contain patch faces" << abort(FatalError);
1732  }
1733 
1734  const face& f = mesh.faces()[facei];
1735 
1736 
1737  forAll(f, fp)
1738  {
1739  if (meshPointMap.found(f[fp]))
1740  {
1741  label patchPointi = meshPointMap[f[fp]];
1742 
1743  if (extrudeStatus[patchPointi] != NOEXTRUDE)
1744  {
1745  unmarkExtrusion
1746  (
1747  patchPointi,
1748  patchDisp,
1749  patchNLayers,
1750  extrudeStatus
1751  );
1752  nChanged++;
1753  }
1754  }
1755  }
1756  }
1757 
1758  forAll(patchDisp, patchPointi)
1759  {
1760  if (mag(patchDisp[patchPointi]) < minThickness[patchPointi])
1761  {
1762  if
1763  (
1764  unmarkExtrusion
1765  (
1766  patchPointi,
1767  patchDisp,
1768  patchNLayers,
1769  extrudeStatus
1770  )
1771  )
1772  {
1773  nChanged++;
1774  }
1775  }
1776  else if (extrudeStatus[patchPointi] == NOEXTRUDE)
1777  {
1778  // Make sure displacement is 0. Should already be so but ...
1779  patchDisp[patchPointi] = Zero;
1780  patchNLayers[patchPointi] = 0;
1781  }
1782  }
1783 
1784 
1785  const faceList& localFaces = pp.localFaces();
1786 
1787  while (true)
1788  {
1789  syncPatchDisplacement
1790  (
1791  pp,
1792  minThickness,
1793  patchDisp,
1794  patchNLayers,
1795  extrudeStatus
1796  );
1797 
1798 
1799  // Pinch
1800  // ~~~~~
1801 
1802  // Make sure that a face doesn't have two non-consecutive areas
1803  // not extruded (e.g. quad where vertex 0 and 2 are not extruded
1804  // but 1 and 3 are) since this gives topological errors.
1805 
1806  label nPinched = 0;
1807 
1808  forAll(localFaces, i)
1809  {
1810  const face& localF = localFaces[i];
1811 
1812  // Count number of transitions from unsnapped to snapped.
1813  label nTrans = 0;
1814 
1815  extrudeMode prevMode = extrudeStatus[localF.prevLabel(0)];
1816 
1817  forAll(localF, fp)
1818  {
1819  extrudeMode fpMode = extrudeStatus[localF[fp]];
1820 
1821  if (prevMode == NOEXTRUDE && fpMode != NOEXTRUDE)
1822  {
1823  nTrans++;
1824  }
1825  prevMode = fpMode;
1826  }
1827 
1828  if (nTrans > 1)
1829  {
1830  // Multiple pinches. Reset whole face as unextruded.
1831  if
1832  (
1833  unmarkExtrusion
1834  (
1835  localF,
1836  patchDisp,
1837  patchNLayers,
1838  extrudeStatus
1839  )
1840  )
1841  {
1842  nPinched++;
1843  nChanged++;
1844  }
1845  }
1846  }
1847 
1848  reduce(nPinched, sumOp<label>());
1849 
1850  Info<< "truncateDisplacement : Unextruded " << nPinched
1851  << " faces due to non-consecutive vertices being extruded." << endl;
1852 
1853 
1854  // Butterfly
1855  // ~~~~~~~~~
1856 
1857  // Make sure that a string of edges becomes a single face so
1858  // not a butterfly. Occasionally an 'edge' will have a single dangling
1859  // vertex due to face combining. These get extruded as a single face
1860  // (with a dangling vertex) so make sure this extrusion forms a single
1861  // shape.
1862  // - continuous i.e. no butterfly:
1863  // + +
1864  // |\ /|
1865  // | \ / |
1866  // +--+--+
1867  // - extrudes from all but the endpoints i.e. no partial
1868  // extrude
1869  // +
1870  // /|
1871  // / |
1872  // +--+--+
1873  // The common error topology is a pinch somewhere in the middle
1874  label nButterFly = 0;
1875  {
1876  DynamicList<label> stringedVerts;
1877  forAll(pp.edges(), edgei)
1878  {
1879  const labelList& globFaces = edgeGlobalFaces[edgei];
1880 
1881  if (globFaces.size() == 2)
1882  {
1883  label myFacei = pp.edgeFaces()[edgei][0];
1884  label myGlobalFacei = globalFaces.toGlobal
1885  (
1886  pp.addressing()[myFacei]
1887  );
1888  label nbrGlobalFacei =
1889  (
1890  globFaces[0] != myGlobalFacei
1891  ? globFaces[0]
1892  : globFaces[1]
1893  );
1894  getVertexString
1895  (
1896  pp,
1897  edgeGlobalFaces,
1898  myFacei,
1899  edgei,
1900  myGlobalFacei,
1901  nbrGlobalFacei,
1902  stringedVerts
1903  );
1904 
1905  if
1906  (
1907  extrudeStatus[stringedVerts[0]] != NOEXTRUDE
1908  || extrudeStatus[stringedVerts.last()] != NOEXTRUDE
1909  )
1910  {
1911  // Any pinch in the middle
1912  bool pinch = false;
1913  for (label i = 1; i < stringedVerts.size()-1; i++)
1914  {
1915  if
1916  (
1917  extrudeStatus[stringedVerts[i]] == NOEXTRUDE
1918  )
1919  {
1920  pinch = true;
1921  break;
1922  }
1923  }
1924  if (pinch)
1925  {
1926  forAll(stringedVerts, i)
1927  {
1928  if
1929  (
1930  unmarkExtrusion
1931  (
1932  stringedVerts[i],
1933  patchDisp,
1934  patchNLayers,
1935  extrudeStatus
1936  )
1937  )
1938  {
1939  nButterFly++;
1940  nChanged++;
1941  }
1942  }
1943  }
1944  }
1945  }
1946  }
1947  }
1948 
1949  reduce(nButterFly, sumOp<label>());
1950 
1951  Info<< "truncateDisplacement : Unextruded " << nButterFly
1952  << " faces due to stringed edges with inconsistent extrusion."
1953  << endl;
1954 
1955 
1956 
1957  // Consistent number of layers
1958  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1959 
1960  // Make sure that a face has consistent number of layers for all
1961  // its vertices.
1962 
1963  label nDiffering = 0;
1964 
1965  // forAll(localFaces, i)
1966  //{
1967  // const face& localF = localFaces[i];
1968  //
1969  // label numLayers = -1;
1970  //
1971  // forAll(localF, fp)
1972  // {
1973  // if (patchNLayers[localF[fp]] > 0)
1974  // {
1975  // if (numLayers == -1)
1976  // {
1977  // numLayers = patchNLayers[localF[fp]];
1978  // }
1979  // else if (numLayers != patchNLayers[localF[fp]])
1980  // {
1981  // // Differing number of layers
1982  // if
1983  // (
1984  // unmarkExtrusion
1985  // (
1986  // localF,
1987  // patchDisp,
1988  // patchNLayers,
1989  // extrudeStatus
1990  // )
1991  // )
1992  // {
1993  // nDiffering++;
1994  // nChanged++;
1995  // }
1996  // break;
1997  // }
1998  // }
1999  // }
2000  //}
2001  //
2002  // reduce(nDiffering, sumOp<label>());
2003  //
2004  // Info<< "truncateDisplacement : Unextruded " << nDiffering
2005  // << " faces due to having differing number of layers." << endl;
2006 
2007  if (nPinched+nButterFly+nDiffering == 0)
2008  {
2009  break;
2010  }
2011  }
2012 
2013  return nChanged;
2014 }
2015 
2016 
2017 // Setup layer information (at points and faces) to modify mesh topology in
2018 // regions where layer mesh terminates.
2019 void Foam::snappyLayerDriver::setupLayerInfoTruncation
2020 (
2021  const indirectPrimitivePatch& pp,
2022  const labelList& patchNLayers,
2023  const List<extrudeMode>& extrudeStatus,
2024  const label nBufferCellsNoExtrude,
2025  labelList& nPatchPointLayers,
2026  labelList& nPatchFaceLayers
2027 ) const
2028 {
2029  Info<< nl << "Setting up information for layer truncation ..." << endl;
2030 
2031  const fvMesh& mesh = meshRefiner_.mesh();
2032 
2033  if (nBufferCellsNoExtrude < 0)
2034  {
2035  Info<< nl << "Performing no layer truncation."
2036  << " nBufferCellsNoExtrude set to less than 0 ..." << endl;
2037 
2038  // Face layers if any point gets extruded
2039  forAll(pp.localFaces(), patchFacei)
2040  {
2041  const face& f = pp.localFaces()[patchFacei];
2042 
2043  forAll(f, fp)
2044  {
2045  if (patchNLayers[f[fp]] > 0)
2046  {
2047  nPatchFaceLayers[patchFacei] = patchNLayers[f[fp]];
2048  break;
2049  }
2050  }
2051  }
2052  nPatchPointLayers = patchNLayers;
2053 
2054  // Set any unset patch face layers
2055  forAll(nPatchFaceLayers, patchFacei)
2056  {
2057  if (nPatchFaceLayers[patchFacei] == -1)
2058  {
2059  nPatchFaceLayers[patchFacei] = 0;
2060  }
2061  }
2062  }
2063  else
2064  {
2065  // Determine max point layers per face.
2066  labelList maxLevel(pp.size(), 0);
2067 
2068  forAll(pp.localFaces(), patchFacei)
2069  {
2070  const face& f = pp.localFaces()[patchFacei];
2071 
2072  // find patch faces where layer terminates (i.e contains extrude
2073  // and noextrude points).
2074 
2075  bool noExtrude = false;
2076  label mLevel = 0;
2077 
2078  forAll(f, fp)
2079  {
2080  if (extrudeStatus[f[fp]] == NOEXTRUDE)
2081  {
2082  noExtrude = true;
2083  }
2084  mLevel = max(mLevel, patchNLayers[f[fp]]);
2085  }
2086 
2087  if (mLevel > 0)
2088  {
2089  // So one of the points is extruded. Check if all are extruded
2090  // or is a mix.
2091 
2092  if (noExtrude)
2093  {
2094  nPatchFaceLayers[patchFacei] = 1;
2095  maxLevel[patchFacei] = mLevel;
2096  }
2097  else
2098  {
2099  maxLevel[patchFacei] = mLevel;
2100  }
2101  }
2102  }
2103 
2104  // We have the seed faces (faces with nPatchFaceLayers != maxLevel)
2105  // Now do a meshwave across the patch where we pick up neighbours
2106  // of seed faces.
2107  // Note: quite inefficient. Could probably be coded better.
2108 
2109  const labelListList& pointFaces = pp.pointFaces();
2110 
2111  label nLevels = gMax(patchNLayers);
2112 
2113  // flag neighbouring patch faces with number of layers to grow
2114  for (label ilevel = 1; ilevel < nLevels; ilevel++)
2115  {
2116  label nBuffer;
2117 
2118  if (ilevel == 1)
2119  {
2120  nBuffer = nBufferCellsNoExtrude - 1;
2121  }
2122  else
2123  {
2124  nBuffer = nBufferCellsNoExtrude;
2125  }
2126 
2127  for (label ibuffer = 0; ibuffer < nBuffer + 1; ibuffer++)
2128  {
2129  labelList tempCounter(nPatchFaceLayers);
2130 
2131  boolList foundNeighbour(pp.nPoints(), false);
2132 
2133  forAll(pp.meshPoints(), patchPointi)
2134  {
2135  forAll(pointFaces[patchPointi], pointFacei)
2136  {
2137  label facei = pointFaces[patchPointi][pointFacei];
2138 
2139  if
2140  (
2141  nPatchFaceLayers[facei] != -1
2142  && maxLevel[facei] > 0
2143  )
2144  {
2145  foundNeighbour[patchPointi] = true;
2146  break;
2147  }
2148  }
2149  }
2150 
2152  (
2153  mesh,
2154  pp.meshPoints(),
2155  foundNeighbour,
2156  orEqOp<bool>(),
2157  false // null value
2158  );
2159 
2160  forAll(pp.meshPoints(), patchPointi)
2161  {
2162  if (foundNeighbour[patchPointi])
2163  {
2164  forAll(pointFaces[patchPointi], pointFacei)
2165  {
2166  label facei = pointFaces[patchPointi][pointFacei];
2167  if
2168  (
2169  nPatchFaceLayers[facei] == -1
2170  && maxLevel[facei] > 0
2171  && ilevel < maxLevel[facei]
2172  )
2173  {
2174  tempCounter[facei] = ilevel;
2175  }
2176  }
2177  }
2178  }
2179  nPatchFaceLayers = tempCounter;
2180  }
2181  }
2182 
2183  forAll(pp.localFaces(), patchFacei)
2184  {
2185  if (nPatchFaceLayers[patchFacei] == -1)
2186  {
2187  nPatchFaceLayers[patchFacei] = maxLevel[patchFacei];
2188  }
2189  }
2190 
2191  forAll(pp.meshPoints(), patchPointi)
2192  {
2193  if (extrudeStatus[patchPointi] != NOEXTRUDE)
2194  {
2195  forAll(pointFaces[patchPointi], pointFacei)
2196  {
2197  label face = pointFaces[patchPointi][pointFacei];
2198  nPatchPointLayers[patchPointi] = max
2199  (
2200  nPatchPointLayers[patchPointi],
2201  nPatchFaceLayers[face]
2202  );
2203  }
2204  }
2205  else
2206  {
2207  nPatchPointLayers[patchPointi] = 0;
2208  }
2209  }
2211  (
2212  mesh,
2213  pp.meshPoints(),
2214  nPatchPointLayers,
2215  maxEqOp<label>(),
2216  label(0) // null value
2217  );
2218  }
2219 }
2220 
2221 
2222 // Does any of the cells use a face from faces?
2223 bool Foam::snappyLayerDriver::cellsUseFace
2224 (
2225  const polyMesh& mesh,
2226  const labelList& cellLabels,
2227  const labelHashSet& faces
2228 )
2229 {
2230  forAll(cellLabels, i)
2231  {
2232  const cell& cFaces = mesh.cells()[cellLabels[i]];
2233 
2234  forAll(cFaces, cFacei)
2235  {
2236  if (faces.found(cFaces[cFacei]))
2237  {
2238  return true;
2239  }
2240  }
2241  }
2242  return false;
2243 }
2244 
2245 
2246 // Checks the newly added cells and locally unmarks points so they
2247 // will not get extruded next time round. Returns global number of unmarked
2248 // points (0 if all was fine)
2249 Foam::label Foam::snappyLayerDriver::checkAndUnmark
2250 (
2251  const addPatchCellLayer& addLayer,
2252  const dictionary& meshQualityDict,
2253  const bool additionalReporting,
2254  const List<labelPair>& baffles,
2255  const indirectPrimitivePatch& pp,
2256  const fvMesh& newMesh,
2257 
2258  pointField& patchDisp,
2259  labelList& patchNLayers,
2260  List<extrudeMode>& extrudeStatus
2261 )
2262 {
2263  // Check the resulting mesh for errors
2264  Info<< nl << "Checking mesh with layer ..." << endl;
2265  faceSet wrongFaces(newMesh, "wrongFaces", newMesh.nFaces()/1000);
2267  (
2268  false,
2269  newMesh,
2270  meshQualityDict,
2271  identity(newMesh.nFaces()),
2272  baffles,
2273  wrongFaces
2274  );
2275  Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
2276  << " illegal faces"
2277  << " (concave, zero area or negative cell pyramid volume)"
2278  << endl;
2279 
2280  // Undo local extrusion if
2281  // - any of the added cells in error
2282 
2283  label nChanged = 0;
2284 
2285  // Get all cells in the layer.
2286  labelListList addedCells
2287  (
2289  (
2290  newMesh,
2291  addLayer.layerFaces()
2292  )
2293  );
2294 
2295  // Check if any of the faces in error uses any face of an added cell
2296  // - if additionalReporting print the few remaining areas for ease of
2297  // finding out where the problems are.
2298 
2299  const label nReportMax = 10;
2300  DynamicField<point> disabledFaceCentres(nReportMax);
2301 
2302  forAll(addedCells, oldPatchFacei)
2303  {
2304  // Get the cells (in newMesh labels) per old patch face (in mesh
2305  // labels)
2306  const labelList& fCells = addedCells[oldPatchFacei];
2307 
2308  if (cellsUseFace(newMesh, fCells, wrongFaces))
2309  {
2310  // Unmark points on old mesh
2311  if
2312  (
2313  unmarkExtrusion
2314  (
2315  pp.localFaces()[oldPatchFacei],
2316  patchDisp,
2317  patchNLayers,
2318  extrudeStatus
2319  )
2320  )
2321  {
2322  if (additionalReporting && (nChanged < nReportMax))
2323  {
2324  disabledFaceCentres.append
2325  (
2326  pp.faceCentres()[oldPatchFacei]
2327  );
2328  }
2329 
2330  nChanged++;
2331  }
2332  }
2333  }
2334 
2335 
2336  label nChangedTotal = returnReduce(nChanged, sumOp<label>());
2337 
2338  if (additionalReporting)
2339  {
2340  // Limit the number of points to be printed so that
2341  // not too many points are reported when running in parallel
2342  // Not accurate, i.e. not always nReportMax points are written,
2343  // but this estimation avoid some communication here.
2344  // The important thing, however, is that when only a few faces
2345  // are disabled, their coordinates are printed, and this should be
2346  // the case
2347  label nReportLocal = nChanged;
2348  if (nChangedTotal > nReportMax)
2349  {
2350  nReportLocal = min
2351  (
2352  max(nChangedTotal / Pstream::nProcs(), 1),
2353  min
2354  (
2355  nChanged,
2356  max(nReportMax / Pstream::nProcs(), 1)
2357  )
2358  );
2359  }
2360 
2361  if (nReportLocal)
2362  {
2363  Pout<< "Checked mesh with layers. Disabled extrusion at " << endl;
2364  for (label i=0; i < nReportLocal; i++)
2365  {
2366  Pout<< " " << disabledFaceCentres[i] << endl;
2367  }
2368  }
2369 
2370  label nReportTotal = returnReduce(nReportLocal, sumOp<label>());
2371 
2372  if (nReportTotal < nChangedTotal)
2373  {
2374  Info<< "Suppressed disabled extrusion message for other "
2375  << nChangedTotal - nReportTotal << " faces." << endl;
2376  }
2377  }
2378 
2379  return nChangedTotal;
2380 }
2381 
2382 
2383 Foam::label Foam::snappyLayerDriver::countExtrusion
2384 (
2385  const indirectPrimitivePatch& pp,
2386  const List<extrudeMode>& extrudeStatus
2387 )
2388 {
2389  // Count number of extruded patch faces
2390  label nExtruded = 0;
2391  {
2392  const faceList& localFaces = pp.localFaces();
2393 
2394  forAll(localFaces, i)
2395  {
2396  const face& localFace = localFaces[i];
2397 
2398  forAll(localFace, fp)
2399  {
2400  if (extrudeStatus[localFace[fp]] != NOEXTRUDE)
2401  {
2402  nExtruded++;
2403  break;
2404  }
2405  }
2406  }
2407  }
2408 
2409  return returnReduce(nExtruded, sumOp<label>());
2410 }
2411 
2412 
2413 // Collect layer faces and layer cells into mesh fields for ease of handling
2414 void Foam::snappyLayerDriver::getLayerCellsFaces
2415 (
2416  const polyMesh& mesh,
2417  const addPatchCellLayer& addLayer,
2418  const scalarField& oldRealThickness,
2419 
2420  labelList& cellNLayers,
2421  scalarField& faceRealThickness
2422 )
2423 {
2424  cellNLayers.setSize(mesh.nCells());
2425  cellNLayers = 0;
2426  faceRealThickness.setSize(mesh.nFaces());
2427  faceRealThickness = 0;
2428 
2429  // Mark all faces in the layer
2430  const labelListList& layerFaces = addLayer.layerFaces();
2431 
2432  // Mark all cells in the layer.
2433  labelListList addedCells(addPatchCellLayer::addedCells(mesh, layerFaces));
2434 
2435  forAll(addedCells, oldPatchFacei)
2436  {
2437  const labelList& added = addedCells[oldPatchFacei];
2438 
2439  const labelList& layer = layerFaces[oldPatchFacei];
2440 
2441  if (layer.size())
2442  {
2443  forAll(added, i)
2444  {
2445  cellNLayers[added[i]] = layer.size()-1;
2446  }
2447  }
2448  }
2449 
2450  forAll(layerFaces, oldPatchFacei)
2451  {
2452  const labelList& layer = layerFaces[oldPatchFacei];
2453  const scalar realThickness = oldRealThickness[oldPatchFacei];
2454 
2455  if (layer.size())
2456  {
2457  // Layer contains both original boundary face and new boundary
2458  // face so is nLayers+1. Leave out old internal face.
2459  for (label i = 1; i < layer.size(); i++)
2460  {
2461  faceRealThickness[layer[i]] = realThickness;
2462  }
2463  }
2464  }
2465 }
2466 
2467 
2468 void Foam::snappyLayerDriver::printLayerData
2469 (
2470  const fvMesh& mesh,
2471  const labelList& patchIDs,
2472  const labelList& cellNLayers,
2473  const scalarField& faceWantedThickness,
2474  const scalarField& faceRealThickness
2475 ) const
2476 {
2477  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
2478 
2479  int oldPrecision = Info().precision();
2480 
2481  // Find maximum length of a patch name, for a nicer output
2482  label maxPatchNameLen = 0;
2483  forAll(patchIDs, i)
2484  {
2485  label patchi = patchIDs[i];
2486  word patchName = pbm[patchi].name();
2487  maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
2488  }
2489 
2490  Info<< nl
2491  << setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
2492  << setw(0) << " faces layers overall thickness" << nl
2493  << setf(ios_base::left) << setw(maxPatchNameLen) << " "
2494  << setw(0) << " [m] [%]" << nl
2495  << setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
2496  << setw(0) << " ----- ------ --- ---" << endl;
2497 
2498 
2499  forAll(patchIDs, i)
2500  {
2501  label patchi = patchIDs[i];
2502  const polyPatch& pp = pbm[patchi];
2503 
2504  label sumSize = pp.size();
2505 
2506  // Number of layers
2507  const labelList& faceCells = pp.faceCells();
2508  label sumNLayers = 0;
2509  forAll(faceCells, i)
2510  {
2511  sumNLayers += cellNLayers[faceCells[i]];
2512  }
2513 
2514  // Thickness
2515  scalarField::subField patchWanted = pbm[patchi].patchSlice
2516  (
2517  faceWantedThickness
2518  );
2519  scalarField::subField patchReal = pbm[patchi].patchSlice
2520  (
2521  faceRealThickness
2522  );
2523 
2524  scalar sumRealThickness = sum(patchReal);
2525  scalar sumFraction = 0;
2526  forAll(patchReal, i)
2527  {
2528  if (patchWanted[i] > vSmall)
2529  {
2530  sumFraction += (patchReal[i]/patchWanted[i]);
2531  }
2532  }
2533 
2534 
2535  reduce(sumSize, sumOp<label>());
2536  reduce(sumNLayers, sumOp<label>());
2537  reduce(sumRealThickness, sumOp<scalar>());
2538  reduce(sumFraction, sumOp<scalar>());
2539 
2540 
2541  scalar avgLayers = 0;
2542  scalar avgReal = 0;
2543  scalar avgFraction = 0;
2544  if (sumSize > 0)
2545  {
2546  avgLayers = scalar(sumNLayers)/sumSize;
2547  avgReal = sumRealThickness/sumSize;
2548  avgFraction = sumFraction/sumSize;
2549  }
2550 
2551  Info<< setf(ios_base::left) << setw(maxPatchNameLen)
2552  << pbm[patchi].name() << setprecision(3)
2553  << " " << setw(8) << sumSize
2554  << " " << setw(8) << avgLayers
2555  << " " << setw(8) << avgReal
2556  << " " << setw(8) << 100*avgFraction
2557  << endl;
2558  }
2559  Info<< setprecision(oldPrecision) << endl;
2560 }
2561 
2562 
2563 bool Foam::snappyLayerDriver::writeLayerData
2564 (
2565  const fvMesh& mesh,
2566  const labelList& patchIDs,
2567  const labelList& cellNLayers,
2568  const scalarField& faceWantedThickness,
2569  const scalarField& faceRealThickness
2570 ) const
2571 {
2572  bool allOk = true;
2573 
2575  {
2576  {
2577  label nAdded = 0;
2578  forAll(cellNLayers, celli)
2579  {
2580  if (cellNLayers[celli] > 0)
2581  {
2582  nAdded++;
2583  }
2584  }
2585  cellSet addedCellSet(mesh, "addedCells", nAdded);
2586  forAll(cellNLayers, celli)
2587  {
2588  if (cellNLayers[celli] > 0)
2589  {
2590  addedCellSet.insert(celli);
2591  }
2592  }
2593  addedCellSet.instance() = meshRefiner_.timeName();
2594  Info<< "Writing "
2595  << returnReduce(addedCellSet.size(), sumOp<label>())
2596  << " added cells to cellSet "
2597  << addedCellSet.name() << endl;
2598  bool ok = addedCellSet.write();
2599  allOk = allOk & ok;
2600  }
2601  {
2602  label nAdded = 0;
2603  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
2604  {
2605  if (faceRealThickness[facei] > 0)
2606  {
2607  nAdded++;
2608  }
2609  }
2610 
2611  faceSet layerFacesSet(mesh, "layerFaces", nAdded);
2612  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
2613  {
2614  if (faceRealThickness[facei] > 0)
2615  {
2616  layerFacesSet.insert(facei);
2617  }
2618  }
2619  layerFacesSet.instance() = meshRefiner_.timeName();
2620  Info<< "Writing "
2621  << returnReduce(layerFacesSet.size(), sumOp<label>())
2622  << " faces inside added layer to faceSet "
2623  << layerFacesSet.name() << endl;
2624  bool ok = layerFacesSet.write();
2625  allOk = allOk & ok;
2626  }
2627  }
2628 
2630  {
2631  Info<< nl << "Writing fields with layer information:" << incrIndent
2632  << endl;
2633  {
2635  (
2636  IOobject
2637  (
2638  "nSurfaceLayers",
2639  mesh.time().timeName(),
2640  mesh,
2643  false
2644  ),
2645  mesh,
2647  fixedValueFvPatchScalarField::typeName
2648  );
2649  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
2650 
2651  volScalarField::Boundary& fldBf =
2652  fld.boundaryFieldRef();
2653 
2654  forAll(patchIDs, i)
2655  {
2656  label patchi = patchIDs[i];
2657  const polyPatch& pp = pbm[patchi];
2658  const labelList& faceCells = pp.faceCells();
2659  scalarField pfld(faceCells.size());
2660  forAll(faceCells, i)
2661  {
2662  pfld[i] = cellNLayers[faceCells[i]];
2663  }
2664  fldBf[patchi] == pfld;
2665  }
2666  Info<< indent << fld.name() << " : actual number of layers"
2667  << endl;
2668  bool ok = fld.write();
2669  allOk = allOk & ok;
2670  }
2671  {
2673  (
2674  IOobject
2675  (
2676  "thickness",
2677  mesh.time().timeName(),
2678  mesh,
2681  false
2682  ),
2683  mesh,
2685  fixedValueFvPatchScalarField::typeName
2686  );
2687 
2688  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
2689 
2690  volScalarField::Boundary& fldBf =
2691  fld.boundaryFieldRef();
2692 
2693  forAll(patchIDs, i)
2694  {
2695  label patchi = patchIDs[i];
2696  fldBf[patchi] == pbm[patchi].patchSlice
2697  (
2698  faceRealThickness
2699  );
2700  }
2701  Info<< indent << fld.name() << " : overall layer thickness"
2702  << endl;
2703  bool ok = fld.write();
2704  allOk = allOk & ok;
2705  }
2706  {
2708  (
2709  IOobject
2710  (
2711  "thicknessFraction",
2712  mesh.time().timeName(),
2713  mesh,
2716  false
2717  ),
2718  mesh,
2720  fixedValueFvPatchScalarField::typeName
2721  );
2722 
2723  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
2724 
2725  volScalarField::Boundary& fldBf =
2726  fld.boundaryFieldRef();
2727 
2728  forAll(patchIDs, i)
2729  {
2730  label patchi = patchIDs[i];
2731 
2732  scalarField::subField patchWanted = pbm[patchi].patchSlice
2733  (
2734  faceWantedThickness
2735  );
2736  scalarField::subField patchReal = pbm[patchi].patchSlice
2737  (
2738  faceRealThickness
2739  );
2740 
2741  // Convert patchReal to relative thickness
2742  scalarField pfld(patchReal.size(), 0.0);
2743  forAll(patchReal, i)
2744  {
2745  if (patchWanted[i] > vSmall)
2746  {
2747  pfld[i] = patchReal[i]/patchWanted[i];
2748  }
2749  }
2750 
2751  fldBf[patchi] == pfld;
2752  }
2753  Info<< indent << fld.name()
2754  << " : overall layer thickness (fraction"
2755  << " of desired thickness)" << endl;
2756  bool ok = fld.write();
2757  allOk = allOk & ok;
2758  }
2759  Info<< decrIndent<< endl;
2760  }
2761 
2762  // if (meshRefinement::outputLevel() & meshRefinement::OUTPUTLAYERINFO)
2763  {
2764  printLayerData
2765  (
2766  mesh,
2767  patchIDs,
2768  cellNLayers,
2769  faceWantedThickness,
2770  faceRealThickness
2771  );
2772  }
2773 
2774  return allOk;
2775 }
2776 
2777 
2778 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2779 
2782  meshRefinement& meshRefiner,
2783  const labelList& globalToMasterPatch,
2784  const labelList& globalToSlavePatch
2785 )
2786 :
2787  meshRefiner_(meshRefiner),
2788  globalToMasterPatch_(globalToMasterPatch),
2789  globalToSlavePatch_(globalToSlavePatch)
2790 {}
2791 
2792 
2793 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2794 
2797  const layerParameters& layerParams,
2798  const dictionary& motionDict
2799 )
2800 {
2801  // Clip to 30 degrees. Not helpful!
2802  // scalar planarAngle = min(30.0, layerParams.featureAngle());
2803  scalar planarAngle = layerParams.featureAngle();
2804  scalar minCos = Foam::cos(degToRad(planarAngle));
2805 
2806  scalar concaveCos = Foam::cos(degToRad(layerParams.concaveAngle()));
2807 
2808  Info<< nl
2809  << "Merging all faces of a cell" << nl
2810  << "---------------------------" << nl
2811  << " - which are on the same patch" << nl
2812  << " - which make an angle < " << planarAngle
2813  << " degrees"
2814  << nl
2815  << " (cos:" << minCos << ')' << nl
2816  << " - as long as the resulting face doesn't become concave"
2817  << " by more than "
2818  << layerParams.concaveAngle() << " degrees" << nl
2819  << " (0=straight, 180=fully concave)" << nl
2820  << endl;
2821 
2822  const fvMesh& mesh = meshRefiner_.mesh();
2823 
2825 
2826  labelList duplicateFace(mesh.nFaces(), -1);
2827  forAll(couples, i)
2828  {
2829  const labelPair& cpl = couples[i];
2830  duplicateFace[cpl[0]] = cpl[1];
2831  duplicateFace[cpl[1]] = cpl[0];
2832  }
2833 
2834  // Get a set of which patches are to have faces merged
2835  labelHashSet patchIDs(meshRefiner_.meshedPatches());
2836  forAll(mesh.boundaryMesh(), patchi)
2837  {
2838  if (layerParams.mergeFaces()[patchi] == layerParameters::mergeFace::no)
2839  {
2840  patchIDs.unset(patchi);
2841  }
2842  if (layerParams.mergeFaces()[patchi] == layerParameters::mergeFace::yes)
2843  {
2844  patchIDs.set(patchi);
2845  }
2846  }
2847 
2848  label nChanged = meshRefiner_.mergePatchFacesUndo
2849  (
2850  minCos,
2851  concaveCos,
2852  patchIDs,
2853  motionDict,
2854  duplicateFace
2855  );
2856 
2857  nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
2858 }
2859 
2860 
2863  const layerParameters& layerParams,
2864  const dictionary& motionDict,
2865  const labelList& patchIDs,
2866  const label nAllowableErrors,
2867  decompositionMethod& decomposer,
2868  fvMeshDistribute& distributor
2869 )
2870 {
2871  fvMesh& mesh = meshRefiner_.mesh();
2872 
2873  // Create baffles (pairs of faces that share the same points)
2874  // Baffles stored as owner and neighbour face that have been created.
2875  List<labelPair> baffles;
2876  meshRefiner_.createZoneBaffles
2877  (
2878  globalToMasterPatch_,
2879  globalToSlavePatch_,
2880  baffles
2881  );
2882 
2883  if (debug&meshRefinement::MESH)
2884  {
2885  const_cast<Time&>(mesh.time())++;
2886  Info<< "Writing baffled mesh to time "
2887  << meshRefiner_.timeName() << endl;
2888  meshRefiner_.write
2889  (
2892  (
2895  ),
2896  mesh.time().path()/meshRefiner_.timeName()
2897  );
2898  }
2899 
2900 
2902  (
2904  (
2905  mesh,
2906  patchIDs
2907  )
2908  );
2909 
2910 
2911  // Global face indices engine
2912  const globalIndex globalFaces(mesh.nFaces());
2913 
2914  // Determine extrudePatch.edgeFaces in global numbering (so across
2915  // coupled patches). This is used only to string up edges between coupled
2916  // faces (all edges between same (global)face indices get extruded).
2917  labelListList edgeGlobalFaces
2918  (
2920  (
2921  mesh,
2922  globalFaces,
2923  pp
2924  )
2925  );
2926 
2927  // Determine patches for extruded boundary edges. Calculates if any
2928  // additional processor patches need to be constructed.
2929 
2930  labelList sidePatchID;
2931  determineSidePatches
2932  (
2933  globalFaces,
2934  edgeGlobalFaces,
2935  pp,
2936 
2937  sidePatchID
2938  );
2939 
2940 
2941  // Point-wise extrusion data
2942  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2943 
2944  // Displacement for all pp.localPoints.
2945  vectorField patchDisp(pp().nPoints(), vector::one);
2946 
2947  // Number of layers for all pp.localPoints. Note: only valid if
2948  // extrudeStatus = EXTRUDE.
2949  labelList patchNLayers(pp().nPoints(), 0);
2950 
2951  // Ideal number of cells added
2952  label nIdealTotAddedCells = 0;
2953 
2954  // Whether to add edge for all pp.localPoints.
2955  List<extrudeMode> extrudeStatus(pp().nPoints(), EXTRUDE);
2956 
2957 
2958  {
2959  // Get number of layer per point from number of layers per patch
2960  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2961 
2962  setNumLayers
2963  (
2964  layerParams.numLayers(), // per patch the num layers
2965  patchIDs, // patches that are being moved
2966  pp, // indirectpatch for all faces moving
2967 
2968  patchDisp,
2969  patchNLayers,
2970  extrudeStatus,
2971  nIdealTotAddedCells
2972  );
2973 
2974  // Precalculate mesh edge labels for patch edges
2975  labelList meshEdges(pp().meshEdges(mesh.edges(), mesh.pointEdges()));
2976 
2977  // Disable extrusion on non-manifold points
2978  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2979 
2980  handleNonManifolds
2981  (
2982  pp,
2983  meshEdges,
2984  edgeGlobalFaces,
2985 
2986  patchDisp,
2987  patchNLayers,
2988  extrudeStatus
2989  );
2990 
2991  // Disable extrusion on feature angles
2992  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2993 
2994  handleFeatureAngle
2995  (
2996  pp,
2997  meshEdges,
2998  degToRad(layerParams.featureAngle()),
2999 
3000  patchDisp,
3001  patchNLayers,
3002  extrudeStatus
3003  );
3004 
3005  // Disable extrusion on warped faces
3006  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3007  // It is hard to calculate some length scale if not in relative
3008  // mode so disable this check.
3009  if (layerParams.relativeSizes())
3010  {
3011  // Undistorted edge length
3012  const scalar edge0Len =
3013  meshRefiner_.meshCutter().level0EdgeLength();
3014  const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
3015 
3016  handleWarpedFaces
3017  (
3018  pp,
3019  layerParams.maxFaceThicknessRatio(),
3020  edge0Len,
3021  cellLevel,
3022 
3023  patchDisp,
3024  patchNLayers,
3025  extrudeStatus
3026  );
3027  }
3028 
3031  //
3032  // handleMultiplePatchFaces
3033  //(
3034  // pp,
3035  //
3036  // patchDisp,
3037  // patchNLayers,
3038  // extrudeStatus
3039  //);
3040 
3041 
3042  // Grow out region of non-extrusion
3043  for (label i = 0; i < layerParams.nGrow(); i++)
3044  {
3045  growNoExtrusion
3046  (
3047  pp,
3048  patchDisp,
3049  patchNLayers,
3050  extrudeStatus
3051  );
3052  }
3053  }
3054 
3055 
3056  // Undistorted edge length
3057  const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
3058  const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
3059 
3060  // Determine (wanted) point-wise overall layer thickness and expansion
3061  // ratio
3062  scalarField thickness(pp().nPoints());
3063  scalarIOField minThickness
3064  (
3065  IOobject
3066  (
3067  "minThickness",
3068  meshRefiner_.timeName(),
3069  mesh,
3071  ),
3072  pp().nPoints()
3073  );
3074  scalarField expansionRatio(pp().nPoints());
3075  calculateLayerThickness
3076  (
3077  pp,
3078  patchIDs,
3079  layerParams,
3080  cellLevel,
3081  patchNLayers,
3082  edge0Len,
3083 
3084  thickness,
3085  minThickness,
3086  expansionRatio
3087  );
3088 
3089 
3090 
3091  // Current set of topology changes. (changing mesh clears out
3092  // polyTopoChange)
3093  polyTopoChange savedMeshMod(mesh.boundaryMesh().size());
3094  // Per cell 0 or number of layers in the cell column it is part of
3095  labelList cellNLayers;
3096  // Per face actual overall layer thickness
3097  scalarField faceRealThickness;
3098  // Per face wanted overall layer thickness
3099  scalarField faceWantedThickness(mesh.nFaces(), 0.0);
3100  {
3101  UIndirectList<scalar>(faceWantedThickness, pp().addressing()) =
3102  avgPointData(pp, thickness);
3103  }
3104 
3105 
3106  {
3107  // Overall displacement field
3108  pointVectorField displacement
3109  (
3110  makeLayerDisplacementField
3111  (
3112  pointMesh::New(mesh),
3113  layerParams.numLayers()
3114  )
3115  );
3116 
3117  // Allocate run-time selectable mesh mover
3118  autoPtr<externalDisplacementMeshMover> medialAxisMoverPtr;
3119  {
3120  // Set up controls for meshMover
3121  dictionary combinedDict(layerParams.dict());
3122  // Add mesh quality constraints
3123  combinedDict.merge(motionDict);
3124  // Where to get minThickness from
3125  combinedDict.add("minThicknessName", minThickness.name());
3126 
3127  // Take over patchDisp as boundary conditions on displacement
3128  // pointVectorField
3129  medialAxisMoverPtr = externalDisplacementMeshMover::New
3130  (
3131  layerParams.meshShrinker(),
3132  combinedDict,
3133  baffles,
3134  displacement
3135  );
3136  }
3137 
3138 
3139  // Saved old points
3140  pointField oldPoints(mesh.points());
3141 
3142  for
3143  (
3144  label iteration = 0;
3145  iteration < layerParams.nLayerIter();
3146  iteration++
3147  )
3148  {
3149  Info<< nl
3150  << "Layer addition iteration " << iteration << nl
3151  << "--------------------------" << endl;
3152 
3153 
3154  // Unset the extrusion at the pp.
3155  const dictionary& meshQualityDict =
3156  (
3157  iteration < layerParams.nRelaxedIter()
3158  ? motionDict
3159  : motionDict.subDict("relaxed")
3160  );
3161 
3162  if (iteration >= layerParams.nRelaxedIter())
3163  {
3164  Info<< "Switched to relaxed meshQuality constraints." << endl;
3165  }
3166 
3167 
3168 
3169  // Make sure displacement is equal on both sides of coupled patches.
3170  syncPatchDisplacement
3171  (
3172  pp,
3173  minThickness,
3174  patchDisp,
3175  patchNLayers,
3176  extrudeStatus
3177  );
3178 
3179  // Displacement acc. to pointnormals
3180  getPatchDisplacement
3181  (
3182  pp,
3183  thickness,
3184  minThickness,
3185  patchDisp,
3186  patchNLayers,
3187  extrudeStatus
3188  );
3189 
3190  // Shrink mesh by displacement value first.
3191  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3192 
3193  {
3194  pointField oldPatchPos(pp().localPoints());
3195 
3196  // Take over patchDisp into pointDisplacement field and
3197  // adjust both for multi-patch constraints
3199  (
3200  patchIDs,
3201  pp,
3202  patchDisp,
3203  displacement
3204  );
3205 
3206 
3207  // Move mesh
3208  // ~~~~~~~~~
3209 
3210  // Set up controls for meshMover
3211  dictionary combinedDict(layerParams.dict());
3212  // Add standard quality constraints
3213  combinedDict.merge(motionDict);
3214  // Add relaxed constraints (overrides standard ones)
3215  combinedDict.merge(meshQualityDict);
3216  // Where to get minThickness from
3217  combinedDict.add("minThicknessName", minThickness.name());
3218 
3219  labelList checkFaces(identity(mesh.nFaces()));
3220  medialAxisMoverPtr().move
3221  (
3222  combinedDict,
3223  nAllowableErrors,
3224  checkFaces
3225  );
3226 
3227  pp().clearGeom();
3228 
3229  // Update patchDisp (since not all might have been honoured)
3230  patchDisp = oldPatchPos - pp().localPoints();
3231  }
3232 
3233  // Truncate displacements that are too small (this will do internal
3234  // ones, coupled ones have already been truncated by
3235  // syncPatchDisplacement)
3236  faceSet dummySet(mesh, "wrongPatchFaces", 0);
3237  truncateDisplacement
3238  (
3239  globalFaces,
3240  edgeGlobalFaces,
3241  pp,
3242  minThickness,
3243  dummySet,
3244  patchDisp,
3245  patchNLayers,
3246  extrudeStatus
3247  );
3248 
3249 
3250  // Dump to .obj file for debugging.
3251  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
3252  {
3253  dumpDisplacement
3254  (
3255  mesh.time().path()/"layer_" + meshRefiner_.timeName(),
3256  pp(),
3257  patchDisp,
3258  extrudeStatus
3259  );
3260 
3261  const_cast<Time&>(mesh.time())++;
3262  Info<< "Writing shrunk mesh to time "
3263  << meshRefiner_.timeName() << endl;
3264 
3265  // See comment in snappySnapDriver why we should not remove
3266  // meshPhi using mesh.clearOut().
3267 
3268  meshRefiner_.write
3269  (
3272  (
3275  ),
3276  mesh.time().path()/meshRefiner_.timeName()
3277  );
3278  }
3279 
3280 
3281  // Mesh topo change engine
3282  polyTopoChange meshMod(mesh);
3283 
3284  // Grow layer of cells on to patch. Handles zero sized displacement.
3285  addPatchCellLayer addLayer(mesh);
3286 
3287  // Determine per point/per face number of layers to extrude. Also
3288  // handles the slow termination of layers when going switching
3289  // layers
3290 
3291  labelList nPatchPointLayers(pp().nPoints(), -1);
3292  labelList nPatchFaceLayers(pp().size(), -1);
3293  setupLayerInfoTruncation
3294  (
3295  pp,
3296  patchNLayers,
3297  extrudeStatus,
3298  layerParams.nBufferCellsNoExtrude(),
3299  nPatchPointLayers,
3300  nPatchFaceLayers
3301  );
3302 
3303  // Calculate displacement for final layer for addPatchLayer.
3304  // (layer of cells next to the original mesh)
3305  vectorField finalDisp(patchNLayers.size(), Zero);
3306 
3307  forAll(nPatchPointLayers, i)
3308  {
3309  scalar ratio = layerParams.finalLayerThicknessRatio
3310  (
3311  nPatchPointLayers[i],
3312  expansionRatio[i]
3313  );
3314  finalDisp[i] = ratio*patchDisp[i];
3315  }
3316 
3317 
3318  const scalarField invExpansionRatio(1.0/expansionRatio);
3319 
3320  // Add topo regardless of whether extrudeStatus is extruderemove.
3321  // Not add layer if patchDisp is zero.
3322  addLayer.setRefinement
3323  (
3324  globalFaces,
3325  edgeGlobalFaces,
3326 
3327  invExpansionRatio,
3328  pp(),
3329  sidePatchID, // boundary patch for extruded bnd edges
3330  labelList(0), // exposed patchIDs (not used adding layers)
3331  nPatchFaceLayers, // layers per face
3332  nPatchPointLayers, // layers per point
3333  finalDisp, // thickness of layer nearest internal mesh
3334  meshMod
3335  );
3336 
3337  if (debug)
3338  {
3339  const_cast<Time&>(mesh.time())++;
3340  }
3341 
3342  // Store mesh changes for if mesh is correct.
3343  savedMeshMod = meshMod;
3344 
3345 
3346  // With the stored topo changes we create a new mesh so we can
3347  // undo if necessary.
3348 
3349  autoPtr<fvMesh> newMeshPtr;
3350  autoPtr<polyTopoChangeMap> map = meshMod.makeMesh
3351  (
3352  newMeshPtr,
3353  IOobject
3354  (
3355  // mesh.name()+"_layer",
3356  mesh.name(),
3357  static_cast<polyMesh&>(mesh).instance(),
3358  mesh.time(), // register with runTime
3360  static_cast<polyMesh&>(mesh).writeOpt()
3361  ), // io params from original mesh but new name
3362  mesh, // original mesh
3363  true // parallel sync
3364  );
3365  fvMesh& newMesh = newMeshPtr();
3366 
3367  //?necessary? Update fields
3368  newMesh.topoChange(map);
3369 
3370  newMesh.setInstance(meshRefiner_.timeName());
3371 
3372  // Update numbering on addLayer:
3373  // - cell/point labels to be newMesh.
3374  // - patchFaces to remain in oldMesh order.
3375  addLayer.topoChange
3376  (
3377  map,
3378  identity(pp().size()),
3379  identity(pp().nPoints())
3380  );
3381 
3382  // Update numbering of baffles
3383  List<labelPair> newMeshBaffles(baffles.size());
3384  forAll(baffles, i)
3385  {
3386  const labelPair& p = baffles[i];
3387  newMeshBaffles[i][0] = map().reverseFaceMap()[p[0]];
3388  newMeshBaffles[i][1] = map().reverseFaceMap()[p[1]];
3389  }
3390 
3391  // Collect layer faces and cells for outside loop.
3392  getLayerCellsFaces
3393  (
3394  newMesh,
3395  addLayer,
3396  avgPointData(pp, mag(patchDisp))(), // current thickness
3397 
3398  cellNLayers,
3399  faceRealThickness
3400  );
3401 
3402 
3403  // Count number of added cells
3404  label nAddedCells = 0;
3405  forAll(cellNLayers, celli)
3406  {
3407  if (cellNLayers[celli] > 0)
3408  {
3409  nAddedCells++;
3410  }
3411  }
3412 
3413 
3414  if (debug&meshRefinement::MESH)
3415  {
3416  Info<< "Writing layer mesh to time " << meshRefiner_.timeName()
3417  << endl;
3418  newMesh.write();
3419 
3420  cellSet addedCellSet(newMesh, "addedCells", nAddedCells);
3421  forAll(cellNLayers, celli)
3422  {
3423  if (cellNLayers[celli] > 0)
3424  {
3425  addedCellSet.insert(celli);
3426  }
3427  }
3428  addedCellSet.instance() = meshRefiner_.timeName();
3429  Info<< "Writing "
3430  << returnReduce(addedCellSet.size(), sumOp<label>())
3431  << " added cells to cellSet " << addedCellSet.name()
3432  << endl;
3433  addedCellSet.write();
3434 
3435  faceSet layerFacesSet
3436  (
3437  newMesh,
3438  "layerFaces",
3439  newMesh.nFaces()/100
3440  );
3441  for (label facei = 0; facei < newMesh.nInternalFaces(); facei++)
3442  {
3443  if (faceRealThickness[facei] > 0)
3444  {
3445  layerFacesSet.insert(facei);
3446  }
3447  }
3448  layerFacesSet.instance() = meshRefiner_.timeName();
3449  Info<< "Writing "
3450  << returnReduce(layerFacesSet.size(), sumOp<label>())
3451  << " faces inside added layer to faceSet "
3452  << layerFacesSet.name() << endl;
3453  layerFacesSet.write();
3454  }
3455 
3456 
3457  label nTotChanged = checkAndUnmark
3458  (
3459  addLayer,
3460  meshQualityDict,
3461  layerParams.additionalReporting(),
3462  newMeshBaffles,
3463  pp(),
3464  newMesh,
3465 
3466  patchDisp,
3467  patchNLayers,
3468  extrudeStatus
3469  );
3470 
3471  label nTotExtruded = countExtrusion(pp, extrudeStatus);
3472  label nTotFaces = returnReduce(pp().size(), sumOp<label>());
3473  label nTotAddedCells = returnReduce(nAddedCells, sumOp<label>());
3474 
3475  Info<< "Extruding " << nTotExtruded
3476  << " out of " << nTotFaces
3477  << " faces (" << 100.0*nTotExtruded/nTotFaces << "%)."
3478  << " Removed extrusion at " << nTotChanged << " faces."
3479  << endl
3480  << "Added " << nTotAddedCells << " out of "
3481  << nIdealTotAddedCells << " cells ("
3482  << 100.0*nTotAddedCells/nIdealTotAddedCells << "%)."
3483  << endl;
3484 
3485  if (nTotChanged == 0)
3486  {
3487  break;
3488  }
3489 
3490  // Reset mesh points and start again
3491  mesh.movePoints(oldPoints);
3492  pp().clearGeom();
3493 
3494  // Grow out region of non-extrusion
3495  for (label i = 0; i < layerParams.nGrow(); i++)
3496  {
3497  growNoExtrusion
3498  (
3499  pp,
3500  patchDisp,
3501  patchNLayers,
3502  extrudeStatus
3503  );
3504  }
3505 
3506  Info<< endl;
3507  }
3508  }
3509 
3510  // At this point we have a (shrunk) mesh and a set of topology changes
3511  // which will make a valid mesh with layer. Apply these changes to the
3512  // current mesh.
3513 
3514  // Apply the stored topo changes to the current mesh.
3515  autoPtr<polyTopoChangeMap> map = savedMeshMod.changeMesh(mesh, false);
3516 
3517  // Hack to remove meshPhi/V0 - mapped incorrectly. TBD.
3518  mesh.clearOut();
3519 
3520  // Update fields
3521  mesh.topoChange(map);
3522 
3523  // Move mesh (since morphing does not do this)
3524  if (map().hasMotionPoints())
3525  {
3526  mesh.movePoints(map().preMotionPoints());
3527  }
3528  else
3529  {
3530  // Delete mesh volumes.
3531  mesh.clearOut();
3532  }
3533 
3534  // Reset the instance for if in overwrite mode
3535  mesh.setInstance(meshRefiner_.timeName());
3536 
3537  meshRefiner_.topoChange(map, labelList(0));
3538 
3539  // Update numbering of faceWantedThickness
3540  meshRefinement::updateList(map().faceMap(), scalar(0), faceWantedThickness);
3541 
3542  // Update numbering on baffles
3543  forAll(baffles, i)
3544  {
3545  labelPair& p = baffles[i];
3546  p[0] = map().reverseFaceMap()[p[0]];
3547  p[1] = map().reverseFaceMap()[p[1]];
3548  }
3549 
3550 
3551  label nBaffles = returnReduce(baffles.size(), sumOp<label>());
3552  if (nBaffles > 0)
3553  {
3554  // Merge any baffles
3555  Info<< "Converting " << nBaffles
3556  << " baffles back into zoned faces ..."
3557  << endl;
3558 
3559  autoPtr<polyTopoChangeMap> map = meshRefiner_.mergeBaffles(baffles);
3560 
3561  inplaceReorder(map().reverseCellMap(), cellNLayers);
3562  inplaceReorder(map().reverseFaceMap(), faceWantedThickness);
3563  inplaceReorder(map().reverseFaceMap(), faceRealThickness);
3564 
3565  Info<< "Converted baffles in = "
3566  << meshRefiner_.mesh().time().cpuTimeIncrement()
3567  << " s\n" << nl << endl;
3568  }
3569 
3570 
3571  // Do final balancing
3572  // ~~~~~~~~~~~~~~~~~~
3573 
3574  if (Pstream::parRun())
3575  {
3576  Info<< nl
3577  << "Doing final balancing" << nl
3578  << "---------------------" << nl
3579  << endl;
3580 
3581  if (debug)
3582  {
3583  const_cast<Time&>(mesh.time())++;
3584  }
3585 
3586  // Hack to remove meshPhi/V0 - mapped incorrectly. TBD.
3587  mesh.clearOut();
3588 
3589  // Balance. No restriction on face zones and baffles.
3590  autoPtr<polyDistributionMap> map = meshRefiner_.balance
3591  (
3592  false,
3593  false,
3594  scalarField(mesh.nCells(), 1.0),
3595  decomposer,
3596  distributor
3597  );
3598 
3599  // Re-distribute flag of layer faces and cells
3600  map().distributeCellData(cellNLayers);
3601  map().distributeFaceData(faceWantedThickness);
3602  map().distributeFaceData(faceRealThickness);
3603  }
3604 
3605 
3606  // Write mesh data
3607  // ~~~~~~~~~~~~~~~
3608 
3609  writeLayerData
3610  (
3611  mesh,
3612  patchIDs,
3613  cellNLayers,
3614  faceWantedThickness,
3615  faceRealThickness
3616  );
3617 }
3618 
3619 
3622  const dictionary& shrinkDict,
3623  const dictionary& motionDict,
3624  const layerParameters& layerParams,
3625  const bool preBalance,
3626  decompositionMethod& decomposer,
3627  fvMeshDistribute& distributor
3628 )
3629 {
3630  const fvMesh& mesh = meshRefiner_.mesh();
3631 
3632  Info<< nl
3633  << "Shrinking and layer addition phase" << nl
3634  << "----------------------------------" << nl
3635  << endl;
3636 
3637  Info<< "Using mesh parameters " << motionDict << nl << endl;
3638 
3639  // Merge coplanar boundary faces
3640  mergePatchFacesUndo(layerParams, motionDict);
3641 
3642  // Per patch the number of layers (-1 or 0 if no layer)
3643  const labelList& numLayers = layerParams.numLayers();
3644 
3645  // Patches that need to get a layer
3646  DynamicList<label> patchIDs(numLayers.size());
3647  label nFacesWithLayers = 0;
3648  forAll(numLayers, patchi)
3649  {
3650  if (numLayers[patchi] > 0)
3651  {
3652  const polyPatch& pp = mesh.boundaryMesh()[patchi];
3653 
3654  if (!pp.coupled())
3655  {
3656  patchIDs.append(patchi);
3657  nFacesWithLayers += mesh.boundaryMesh()[patchi].size();
3658  }
3659  else
3660  {
3662  << "Ignoring layers on coupled patch " << pp.name()
3663  << endl;
3664  }
3665  }
3666  }
3667  patchIDs.shrink();
3668 
3669  if (returnReduce(nFacesWithLayers, sumOp<label>()) == 0)
3670  {
3671  Info<< nl << "No layers to generate ..." << endl;
3672  }
3673  else
3674  {
3675  // Check that outside of mesh is not multiply connected.
3676  checkMeshManifold();
3677 
3678  // Check initial mesh
3679  Info<< "Checking initial mesh ..." << endl;
3680  labelHashSet wrongFaces(mesh.nFaces()/100);
3681  motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
3682  const label nInitErrors = returnReduce
3683  (
3684  wrongFaces.size(),
3685  sumOp<label>()
3686  );
3687 
3688  Info<< "Detected " << nInitErrors << " illegal faces"
3689  << " (concave, zero area or negative cell pyramid volume)"
3690  << endl;
3691 
3692 
3693  // Balance
3694  if (Pstream::parRun() && preBalance)
3695  {
3696  Info<< nl
3697  << "Doing initial balancing" << nl
3698  << "-----------------------" << nl
3699  << endl;
3700 
3701  scalarField cellWeights(mesh.nCells(), 1);
3702  forAll(numLayers, patchi)
3703  {
3704  if (numLayers[patchi] > 0)
3705  {
3706  const polyPatch& pp = mesh.boundaryMesh()[patchi];
3707  forAll(pp.faceCells(), i)
3708  {
3709  cellWeights[pp.faceCells()[i]] += numLayers[patchi];
3710  }
3711  }
3712  }
3713 
3714  // Balance mesh (and meshRefinement). Restrict faceZones to
3715  // be on internal faces only since they will be converted into
3716  // baffles.
3717  autoPtr<polyDistributionMap> map = meshRefiner_.balance
3718  (
3719  true, // false, // keepZoneFaces
3720  false,
3721  cellWeights,
3722  decomposer,
3723  distributor
3724  );
3725  }
3726 
3727 
3728  // Do all topo changes
3729  addLayers
3730  (
3731  layerParams,
3732  motionDict,
3733  patchIDs,
3734  nInitErrors,
3735  decomposer,
3736  distributor
3737  );
3738  }
3739 }
3740 
3741 
3742 // ************************************************************************* //
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.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void clearOut()
Clear all geometry and addressing.
Definition: fvMesh.C:258
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Simple container to keep together layer specific information.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
virtual Ostream & write(const char)=0
Write character.
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
const word & name() const
Return name.
const word & name() const
Return name.
Definition: IOobject.H:315
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
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:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
scalar concaveAngle() const
label nInternalFaces() const
const labelListList & pointEdges() const
label nFaces() 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.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
List< face > faceList
Definition: faceListFwd.H:43
label nCells() const
static tmp< GeometricField< vector, pointPatchField, pointMesh > > New(const word &name, const Internal &, const PtrList< pointPatchField< vector >> &)
Return a temporary field constructed from name,.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const List< mergeFace > & mergeFaces() const
Whether to merge boundary faces of the same layer cell.
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
label nBufferCellsNoExtrude() const
Create buffer region for new layer terminations.
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.
GeometricBoundaryField< scalar, fvPatchField, volMesh > Boundary
Type of the boundary field.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const dimensionSet dimless
const dictionary & dict() const
autoPtr< polyTopoChangeMap > 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.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:372
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
virtual bool write(const bool write=true) const
Write mesh using IO settings from time.
Definition: fvMesh.C:1658
fvMesh & mesh
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:69
const labelList & numLayers() const
How many layers to add:
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
virtual void topoChange(const polyTopoChangeMap &map)
Update mesh corresponding to the given map.
Definition: fvMesh.C:1151
static const word & calculatedType()
Return the type of the calculated for of fvPatchField.
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)
const dimensionSet dimLength
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1211
label nRelaxedIter() const
Number of iterations after which relaxed motion rules.
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:1002
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:58
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
static void updateList(const labelList &newToOld, const T &nullValue, List< T > &elems)
Helper: reorder list according to map.
scalar maxFaceThicknessRatio() const
Stop layer growth on highly warped cells.
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
snappyLayerDriver(meshRefinement &meshRefiner, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch)
Construct from components.
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
scalar finalLayerThicknessRatio(const label nLayers, const scalar expansionRatio) const
Determine ratio of final layer thickness to.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:316
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)
void topoChange(const polyTopoChangeMap &, const labelList &faceMap, const labelList &pointMap)
Update any locally stored mesh information. Gets additional.
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:340
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh edges.
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.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
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
static const zero Zero
Definition: zero.H:97
Abstract base class for decomposition.
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
errorManip< error > abort(error &err)
Definition: errorManip.H:131
scalar featureAngle() const
label nLayerIter() const
Number of overall layer addition iterations.
Istream and Ostream manipulators taking arguments.
static const char nl
Definition: Ostream.H:260
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:235
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
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)
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:77
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:1065
static tmp< pointField > pointNormals(const polyMesh &, const PrimitivePatch< FaceList, PointField > &)
Return parallel consistent point normals for patches using mesh points.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
List< word > wordList
A List of words.
Definition: fileName.H:54
void setSize(const label)
Reset size of List.
Definition: List.C:281
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
bool relativeSizes() const
Are size parameters relative to inner cell size or.
label patchi
label nGrow() const
If points get not extruded do nGrow layers of connected faces.
vector point
Point is a vector.
Definition: point.H:41
const word & meshShrinker() const
Type of mesh shrinker.
#define WarningInFunction
Report a warning using Foam::Warning.
const word & name() const
Return reference to name.
Definition: fvMesh.H:386
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A collection of cell labels.
Definition: cellSet.H:48
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:95
Direct mesh changes based on v1.3 polyTopoChange syntax.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh points.
const Switch & additionalReporting() const
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:100
bool merge(const dictionary &)
Merge entries from the given dictionary.
Definition: dictionary.C:1438
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.
labelListList addedCells() const
Added cells given current mesh & layerfaces.
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:52
fileName & instance() const
Return the instance directory, constant, system, <time> etc.
Definition: IOobject.C:355
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
virtual bool write(const bool write=true) const
Write using setting from DB.
volScalarField & p
A class for managing temporary objects.
Definition: PtrList.H:53
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: TimePaths.H:139
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:228
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.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98
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.
A primitive field of type <T> with automated input and output.
Definition: IOField.H:50
static const label labelMin
Definition: label.H:61
Namespace for OpenFOAM.