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