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