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