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