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