FreeFOAM The Cross-Platform CFD Toolkit
autoLayerDriverShrink.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) 1991-2010 OpenCFD Ltd.
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  Shrinking mesh (part of adding cell layers)
26 
27 \*----------------------------------------------------------------------------*/
28 
29 #include "autoLayerDriver.H"
30 #include <finiteVolume/fvMesh.H>
31 #include <OpenFOAM/Time.H>
32 #include <OpenFOAM/pointFields.H>
34 #include <autoMesh/pointData.H>
36 
37 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
38 
39 // Calculate inverse sum of edge weights (currently always 1.0)
40 void Foam::autoLayerDriver::sumWeights
41 (
42  const PackedBoolList& isMasterEdge,
43  const labelList& meshEdges,
44  const labelList& meshPoints,
45  const edgeList& edges,
46  scalarField& invSumWeight
47 ) const
48 {
49  invSumWeight = 0;
50 
51  forAll(edges, edgeI)
52  {
53  if (isMasterEdge.get(meshEdges[edgeI]) == 1)
54  {
55  const edge& e = edges[edgeI];
56  //scalar eWeight = edgeWeights[edgeI];
57  scalar eWeight = 1.0;
58 
59  invSumWeight[e[0]] += eWeight;
60  invSumWeight[e[1]] += eWeight;
61  }
62  }
63 
65  (
66  meshRefiner_.mesh(),
67  meshPoints,
68  invSumWeight,
69  plusEqOp<scalar>(),
70  scalar(0.0), // null value
71  false // no separation
72  );
73 
74  forAll(invSumWeight, pointI)
75  {
76  scalar w = invSumWeight[pointI];
77 
78  if (w > 0.0)
79  {
80  invSumWeight[pointI] = 1.0/w;
81  }
82  }
83 }
84 
85 
86 // Smooth field on moving patch
87 void Foam::autoLayerDriver::smoothField
88 (
89  const motionSmoother& meshMover,
90  const PackedBoolList& isMasterEdge,
91  const labelList& meshEdges,
92  const scalarField& fieldMin,
93  const label nSmoothDisp,
94  scalarField& field
95 ) const
96 {
97  const indirectPrimitivePatch& pp = meshMover.patch();
98  const edgeList& edges = pp.edges();
99  const labelList& meshPoints = pp.meshPoints();
100 
101  scalarField invSumWeight(pp.nPoints());
102  sumWeights
103  (
104  isMasterEdge,
105  meshEdges,
106  meshPoints,
107  edges,
108  invSumWeight
109  );
110 
111  // Get smoothly varying patch field.
112  Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
113 
114  for (label iter = 0; iter < nSmoothDisp; iter++)
115  {
116  scalarField average(pp.nPoints());
117  averageNeighbours
118  (
119  meshMover.mesh(),
120  isMasterEdge,
121  meshEdges,
122  meshPoints,
123  pp.edges(),
124  invSumWeight,
125  field,
126  average
127  );
128 
129  // Transfer to field
130  forAll(field, pointI)
131  {
132  //full smoothing neighbours + point value
133  average[pointI] = 0.5*(field[pointI]+average[pointI]);
134 
135  // perform monotonic smoothing
136  if
137  (
138  average[pointI] < field[pointI]
139  && average[pointI] >= fieldMin[pointI]
140  )
141  {
142  field[pointI] = average[pointI];
143  }
144  }
145 
146  // Do residual calculation every so often.
147  if ((iter % 10) == 0)
148  {
149  Info<< " Iteration " << iter << " residual "
150  << gSum(mag(field-average))
151  /returnReduce(average.size(), sumOp<label>())
152  << endl;
153  }
154  }
155 }
156 
157 
158 // Smooth normals on moving patch.
159 void Foam::autoLayerDriver::smoothPatchNormals
160 (
161  const motionSmoother& meshMover,
162  const PackedBoolList& isMasterEdge,
163  const labelList& meshEdges,
164  const label nSmoothDisp,
165  pointField& normals
166 ) const
167 {
168  const indirectPrimitivePatch& pp = meshMover.patch();
169  const edgeList& edges = pp.edges();
170  const labelList& meshPoints = pp.meshPoints();
171 
172  // Calculate inverse sum of weights
173 
174  scalarField invSumWeight(pp.nPoints());
175  sumWeights
176  (
177  isMasterEdge,
178  meshEdges,
179  meshPoints,
180  edges,
181  invSumWeight
182  );
183 
184  // Get smoothly varying internal normals field.
185  Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
186 
187  for (label iter = 0; iter < nSmoothDisp; iter++)
188  {
189  vectorField average(pp.nPoints());
190  averageNeighbours
191  (
192  meshMover.mesh(),
193  isMasterEdge,
194  meshEdges,
195  meshPoints,
196  pp.edges(),
197  invSumWeight,
198  normals,
199  average
200  );
201 
202  // Do residual calculation every so often.
203  if ((iter % 10) == 0)
204  {
205  Info<< " Iteration " << iter << " residual "
206  << gSum(mag(normals-average))
207  /returnReduce(average.size(), sumOp<label>())
208  << endl;
209  }
210 
211  // Transfer to normals vector field
212  forAll(average, pointI)
213  {
214  // full smoothing neighbours + point value
215  average[pointI] = 0.5*(normals[pointI]+average[pointI]);
216  normals[pointI] = average[pointI];
217  normals[pointI] /= mag(normals[pointI]) + VSMALL;
218  }
219  }
220 }
221 
222 
223 // Smooth normals in interior.
224 void Foam::autoLayerDriver::smoothNormals
225 (
226  const label nSmoothDisp,
227  const PackedBoolList& isMasterEdge,
228  const labelList& fixedPoints,
229  pointVectorField& normals
230 ) const
231 {
232  // Get smoothly varying internal normals field.
233  Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
234 
235  const fvMesh& mesh = meshRefiner_.mesh();
236  const edgeList& edges = mesh.edges();
237 
238  // Points that do not change.
239  PackedBoolList isFixedPoint(mesh.nPoints());
240 
241  // Internal points that are fixed
242  forAll(fixedPoints, i)
243  {
244  label meshPointI = fixedPoints[i];
245  isFixedPoint.set(meshPointI, 1);
246  }
247 
248  // Correspondence between local edges/points and mesh edges/points
249  const labelList meshEdges(identity(mesh.nEdges()));
250  const labelList meshPoints(identity(mesh.nPoints()));
251 
252  // Calculate inverse sum of weights
253 
254  scalarField invSumWeight(mesh.nPoints(), 0);
255  sumWeights
256  (
257  isMasterEdge,
258  meshEdges,
259  meshPoints,
260  edges,
261  invSumWeight
262  );
263 
264  Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
265 
266  for (label iter = 0; iter < nSmoothDisp; iter++)
267  {
268  vectorField average(mesh.nPoints());
269  averageNeighbours
270  (
271  mesh,
272  isMasterEdge,
273  meshEdges,
274  meshPoints,
275  edges,
276  invSumWeight,
277  normals,
278  average
279  );
280 
281  // Do residual calculation every so often.
282  if ((iter % 10) == 0)
283  {
284  Info<< " Iteration " << iter << " residual "
285  << gSum(mag(normals-average))
286  /returnReduce(average.size(), sumOp<label>())
287  << endl;
288  }
289 
290 
291  // Transfer to normals vector field
292  forAll(average, pointI)
293  {
294  if (isFixedPoint.get(pointI) == 0)
295  {
296  //full smoothing neighbours + point value
297  average[pointI] = 0.5*(normals[pointI]+average[pointI]);
298  normals[pointI] = average[pointI];
299  normals[pointI] /= mag(normals[pointI]) + VSMALL;
300  }
301  }
302  }
303 }
304 
305 
306 // Tries and find a medial axis point. Done by comparing vectors to nearest
307 // wall point for both vertices of edge.
308 bool Foam::autoLayerDriver::isMaxEdge
309 (
310  const List<pointData>& pointWallDist,
311  const label edgeI,
312  const scalar minCos
313 ) const
314 {
315  const fvMesh& mesh = meshRefiner_.mesh();
316  const pointField& points = mesh.points();
317 
318  // Do not mark edges with one side on moving wall.
319 
320  const edge& e = mesh.edges()[edgeI];
321 
322  vector v0(points[e[0]] - pointWallDist[e[0]].origin());
323  scalar magV0(mag(v0));
324 
325  if (magV0 < SMALL)
326  {
327  return false;
328  }
329 
330  vector v1(points[e[1]] - pointWallDist[e[1]].origin());
331  scalar magV1(mag(v1));
332 
333  if (magV1 < SMALL)
334  {
335  return false;
336  }
337 
338  v0 /= magV0;
339  v1 /= magV1;
340 
341  // Test angle.
342  if ((v0 & v1) < minCos)
343  {
344  return true;
345  }
346  else
347  {
348  return false;
349  }
350 }
351 
352 
353 // Stop layer growth where mesh wraps around edge with a
354 // large feature angle
355 void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
356 (
357  const indirectPrimitivePatch& pp,
358  const scalar minCos,
359  List<extrudeMode>& extrudeStatus,
360  pointField& patchDisp,
361  labelList& patchNLayers,
362  label& nPointCounter
363 ) const
364 {
365  // Mark faces that have all points extruded
366  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
367 
368  boolList extrudedFaces(pp.size(), true);
369 
370  forAll(pp.localFaces(), faceI)
371  {
372  const face& f = pp.localFaces()[faceI];
373 
374  forAll(f, fp)
375  {
376  if (extrudeStatus[f[fp]] == NOEXTRUDE)
377  {
378  extrudedFaces[faceI] = false;
379  break;
380  }
381  }
382  }
383 
384 
385  // Detect situation where two featureedge-neighbouring faces are partly or
386  // not extruded and the edge itself is extruded. In this case unmark the
387  // edge for extrusion.
388 
389  forAll(pp.edgeFaces(), edgeI)
390  {
391  const labelList& eFaces = pp.edgeFaces()[edgeI];
392 
393  if (eFaces.size() == 2)
394  {
395  const edge& e = pp.edges()[edgeI];
396  label v0 = e[0];
397  label v1 = e[1];
398 
399  if
400  (
401  extrudeStatus[v0] != NOEXTRUDE
402  || extrudeStatus[v1] != NOEXTRUDE
403  )
404  {
405  if (!extrudedFaces[eFaces[0]] || !extrudedFaces[eFaces[1]])
406  {
407  const vector& n0 = pp.faceNormals()[eFaces[0]];
408  const vector& n1 = pp.faceNormals()[eFaces[1]];
409 
410  if ((n0 & n1) < minCos)
411  {
412  if
413  (
414  unmarkExtrusion
415  (
416  v0,
417  patchDisp,
418  patchNLayers,
419  extrudeStatus
420  )
421  )
422  {
423  nPointCounter++;
424  }
425  if
426  (
427  unmarkExtrusion
428  (
429  v1,
430  patchDisp,
431  patchNLayers,
432  extrudeStatus
433  )
434  )
435  {
436  nPointCounter++;
437  }
438  }
439  }
440  }
441  }
442  }
443 }
444 
445 
446 // Find isolated islands (points, edges and faces and layer terminations)
447 // in the layer mesh and stop any layer growth at these points.
448 void Foam::autoLayerDriver::findIsolatedRegions
449 (
450  const indirectPrimitivePatch& pp,
451  const PackedBoolList& isMasterEdge,
452  const labelList& meshEdges,
453  const scalar minCosLayerTermination,
454  scalarField& field,
455  List<extrudeMode>& extrudeStatus,
456  pointField& patchDisp,
457  labelList& patchNLayers
458 ) const
459 {
460  const fvMesh& mesh = meshRefiner_.mesh();
461 
462  Info<< "shrinkMeshDistance : Removing isolated regions ..." << endl;
463 
464  // Keep count of number of points unextruded
465  label nPointCounter = 0;
466 
467  while (true)
468  {
469  // Stop layer growth where mesh wraps around edge with a
470  // large feature angle
471  handleFeatureAngleLayerTerminations
472  (
473  pp,
474  minCosLayerTermination,
475 
476  extrudeStatus,
477  patchDisp,
478  patchNLayers,
479  nPointCounter
480  );
481 
482 
483  // Do not extrude from point where all neighbouring
484  // faces are not grown
485  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
486  boolList extrudedFaces(pp.size(), true);
487  forAll(pp.localFaces(), faceI)
488  {
489  const face& f = pp.localFaces()[faceI];
490  forAll(f, fp)
491  {
492  if (extrudeStatus[f[fp]] == NOEXTRUDE)
493  {
494  extrudedFaces[faceI] = false;
495  break;
496  }
497  }
498  }
499 
500  const labelListList& pointFaces = pp.pointFaces();
501 
502  boolList keptPoints(pp.nPoints(), false);
503  forAll(keptPoints, patchPointI)
504  {
505  const labelList& pFaces = pointFaces[patchPointI];
506 
507  forAll(pFaces, i)
508  {
509  label faceI = pFaces[i];
510  if (extrudedFaces[faceI])
511  {
512  keptPoints[patchPointI] = true;
513  break;
514  }
515  }
516  }
517 
519  (
520  mesh,
521  pp.meshPoints(),
522  keptPoints,
523  orEqOp<bool>(),
524  false, // null value
525  false // no separation
526  );
527 
528  label nChanged = 0;
529 
530  forAll(keptPoints, patchPointI)
531  {
532  if (!keptPoints[patchPointI])
533  {
534  if
535  (
536  unmarkExtrusion
537  (
538  patchPointI,
539  patchDisp,
540  patchNLayers,
541  extrudeStatus
542  )
543  )
544  {
545  nPointCounter++;
546  nChanged++;
547  }
548  }
549  }
550 
551 
552  if (returnReduce(nChanged, sumOp<label>()) == 0)
553  {
554  break;
555  }
556  }
557 
558  const edgeList& edges = pp.edges();
559 
560 
561  // Count number of mesh edges using a point
562  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
563 
564  labelList isolatedPoint(pp.nPoints(),0);
565 
566  forAll(edges, edgeI)
567  {
568  if (isMasterEdge.get(meshEdges[edgeI]) == 1)
569  {
570  const edge& e = edges[edgeI];
571 
572  label v0 = e[0];
573  label v1 = e[1];
574 
575  if (extrudeStatus[v1] != NOEXTRUDE)
576  {
577  isolatedPoint[v0] += 1;
578  }
579  if (extrudeStatus[v0] != NOEXTRUDE)
580  {
581  isolatedPoint[v1] += 1;
582  }
583  }
584  }
585 
587  (
588  mesh,
589  pp.meshPoints(),
590  isolatedPoint,
591  plusEqOp<label>(),
592  0, // null value
593  false // no separation
594  );
595 
596  // stop layer growth on isolated faces
597  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
598  forAll(pp, faceI)
599  {
600  const face& f = pp.localFaces()[faceI];
601  bool failed = false;
602  forAll(f, fp)
603  {
604  if (isolatedPoint[f[fp]] > 2)
605  {
606  failed = true;
607  break;
608  }
609  }
610  bool allPointsExtruded = true;
611  if (!failed)
612  {
613  forAll(f, fp)
614  {
615  if (extrudeStatus[f[fp]] == NOEXTRUDE)
616  {
617  allPointsExtruded = false;
618  break;
619  }
620  }
621 
622  if (allPointsExtruded)
623  {
624  forAll(f, fp)
625  {
626  if
627  (
628  unmarkExtrusion
629  (
630  f[fp],
631  patchDisp,
632  patchNLayers,
633  extrudeStatus
634  )
635  )
636  {
637  field[f[fp]] = 0.0;
638  nPointCounter++;
639  }
640  }
641  }
642  }
643  }
644 
645  reduce(nPointCounter, sumOp<label>());
646  Info<< "Number isolated points extrusion stopped : "<< nPointCounter<< endl;
647 }
648 
649 
650 // Calculates medial axis fields:
651 // dispVec : normal of nearest wall. Where this normal changes direction
652 // defines the medial axis
653 // medialDist : distance to medial axis
654 // medialRatio : ratio of medial distance to wall distance.
655 // (1 at wall, 0 at medial axis)
656 void Foam::autoLayerDriver::medialAxisSmoothingInfo
657 (
658  const motionSmoother& meshMover,
659  const label nSmoothNormals,
660  const label nSmoothSurfaceNormals,
661  const scalar minMedianAxisAngleCos,
662 
663  pointVectorField& dispVec,
664  pointScalarField& medialRatio,
665  pointScalarField& medialDist
666 ) const
667 {
668 
669  Info<< "medialAxisSmoothingInfo :"
670  << " Calculate distance to Medial Axis ..." << endl;
671 
672  const polyMesh& mesh = meshMover.mesh();
673  const pointField& points = mesh.points();
674 
675  const indirectPrimitivePatch& pp = meshMover.patch();
676  const vectorField& faceNormals = pp.faceNormals();
677  const labelList& meshPoints = pp.meshPoints();
678 
679  // Predetermine mesh edges
680  // ~~~~~~~~~~~~~~~~~~~~~~~
681 
682  // Precalulate master edge (only relevant for shared edges)
683  PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
684  // Precalculate meshEdge per pp edge
685  labelList meshEdges(pp.nEdges());
686 
687  forAll(meshEdges, patchEdgeI)
688  {
689  const edge& e = pp.edges()[patchEdgeI];
690 
691  label v0 = pp.meshPoints()[e[0]];
692  label v1 = pp.meshPoints()[e[1]];
693  meshEdges[patchEdgeI] = meshTools::findEdge
694  (
695  mesh.edges(),
696  mesh.pointEdges()[v0],
697  v0,
698  v1
699  );
700  }
701 
702 
703  // Determine pointNormal
704  // ~~~~~~~~~~~~~~~~~~~~~
705 
706  pointField pointNormals(pp.nPoints(), vector::zero);
707  {
708  labelList nPointFaces(pp.nPoints(), 0);
709 
710  forAll(faceNormals, faceI)
711  {
712  const face& f = pp.localFaces()[faceI];
713 
714  forAll(f, fp)
715  {
716  pointNormals[f[fp]] += faceNormals[faceI];
717  nPointFaces[f[fp]] ++;
718  }
719  }
720 
722  (
723  mesh,
724  meshPoints,
725  pointNormals,
726  plusEqOp<vector>(),
727  vector::zero, // null value
728  false // no separation
729  );
730 
732  (
733  mesh,
734  meshPoints,
735  nPointFaces,
736  plusEqOp<label>(),
737  0, // null value
738  false // no separation
739  );
740 
741  forAll(pointNormals, i)
742  {
743  pointNormals[i] /= nPointFaces[i];
744  }
745  }
746 
747  // Smooth patch normal vectors
748  smoothPatchNormals
749  (
750  meshMover,
751  isMasterEdge,
752  meshEdges,
753  nSmoothSurfaceNormals,
754  pointNormals
755  );
756 
757 
758  // Calculate distance to pp points
759  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
760 
761  // Distance to wall
762  List<pointData> pointWallDist(mesh.nPoints());
763 
764 
765  // 1. Calculate distance to points where displacement is specified.
766  {
767  // Seed data.
768  List<pointData> wallInfo(meshPoints.size());
769 
770  forAll(meshPoints, patchPointI)
771  {
772  label pointI = meshPoints[patchPointI];
773  wallInfo[patchPointI] = pointData
774  (
775  points[pointI],
776  0.0,
777  pointI, // passive scalar
778  pointNormals[patchPointI] // surface normals
779  );
780  }
781 
782  // Do all calculations
783  List<pointData> edgeWallDist(mesh.nEdges());
784  PointEdgeWave<pointData> wallDistCalc
785  (
786  mesh,
787  meshPoints,
788  wallInfo,
789  pointWallDist,
790  edgeWallDist,
791  mesh.globalData().nTotalPoints() // max iterations
792  );
793  }
794 
795  // 2. Find points with max distance and transport information back to
796  // wall.
797  {
798  List<pointData> pointMedialDist(mesh.nPoints());
799  List<pointData> edgeMedialDist(mesh.nEdges());
800 
801  // Seed point data.
802  DynamicList<pointData> maxInfo(meshPoints.size());
803  DynamicList<label> maxPoints(meshPoints.size());
804 
805  // 1. Medial axis points
806 
807  const edgeList& edges = mesh.edges();
808 
809  forAll(edges, edgeI)
810  {
811  if (isMaxEdge(pointWallDist, edgeI, minMedianAxisAngleCos))
812  {
813  // Both end points of edge have very different nearest wall
814  // point. Mark both points as medial axis points.
815  const edge& e = edges[edgeI];
816 
817  forAll(e, ep)
818  {
819  label pointI = e[ep];
820 
821  if (!pointMedialDist[pointI].valid())
822  {
823  maxPoints.append(pointI);
824  maxInfo.append
825  (
826  pointData
827  (
828  points[pointI],
829  0.0,
830  pointI, // passive data
831  vector::zero // passive data
832  )
833  );
834  pointMedialDist[pointI] = maxInfo[maxInfo.size()-1];
835  }
836  }
837  }
838  }
839 
840 
841  // 2. Seed non-adapt patches
842  const polyBoundaryMesh& patches = mesh.boundaryMesh();
843 
844  labelHashSet adaptPatches(meshMover.adaptPatchIDs());
845 
846  forAll(patches, patchI)
847  {
848  const polyPatch& pp = patches[patchI];
849 
850  if
851  (
852  !pp.coupled()
853  && !isA<emptyPolyPatch>(pp)
854  && !adaptPatches.found(patchI)
855  )
856  {
857  Info<< "Inserting points on patch " << pp.name() << endl;
858 
859  const labelList& meshPoints = pp.meshPoints();
860 
861  forAll(meshPoints, i)
862  {
863  label pointI = meshPoints[i];
864 
865  if (!pointMedialDist[pointI].valid())
866  {
867  maxPoints.append(pointI);
868  maxInfo.append
869  (
870  pointData
871  (
872  points[pointI],
873  0.0,
874  pointI, // passive data
875  vector::zero // passive data
876  )
877  );
878  pointMedialDist[pointI] = maxInfo[maxInfo.size()-1];
879  }
880  }
881  }
882  }
883 
884  maxInfo.shrink();
885  maxPoints.shrink();
886 
887  // Do all calculations
888  PointEdgeWave<pointData> medialDistCalc
889  (
890  mesh,
891  maxPoints,
892  maxInfo,
893 
894  pointMedialDist,
895  edgeMedialDist,
896  mesh.globalData().nTotalPoints() // max iterations
897  );
898 
899  // Extract medial axis distance as pointScalarField
900  forAll(pointMedialDist, pointI)
901  {
902  medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr());
903  }
904  }
905 
906  // Extract transported surface normals as pointVectorField
907  forAll(dispVec, i)
908  {
909  dispVec[i] = pointWallDist[i].v();
910  }
911 
912  // Smooth normal vectors. Do not change normals on pp.meshPoints
913  smoothNormals(nSmoothNormals, isMasterEdge, meshPoints, dispVec);
914 
915  // Calculate ratio point medial distance to point wall distance
916  forAll(medialRatio, pointI)
917  {
918  scalar wDist2 = pointWallDist[pointI].distSqr();
919  scalar mDist = medialDist[pointI];
920 
921  if (wDist2 < sqr(SMALL) && mDist < SMALL)
922  {
923  medialRatio[pointI] = 0.0;
924  }
925  else
926  {
927  medialRatio[pointI] = mDist / (Foam::sqrt(wDist2) + mDist);
928  }
929  }
930 
931  if (debug)
932  {
933  Info<< "medialAxisSmoothingInfo :"
934  << " Writing:" << nl
935  << " " << dispVec.name()
936  << " : normalised direction of nearest displacement" << nl
937  << " " << medialDist.name()
938  << " : distance to medial axis" << nl
939  << " " << medialRatio.name()
940  << " : ratio of medial distance to wall distance" << nl
941  << endl;
942  dispVec.write();
943  medialDist.write();
944  medialRatio.write();
945  }
946 }
947 
948 
949 void Foam::autoLayerDriver::shrinkMeshMedialDistance
950 (
951  motionSmoother& meshMover,
952  const dictionary& meshQualityDict,
953  const label nSmoothThickness,
954  const scalar maxThicknessToMedialRatio,
955  const label nAllowableErrors,
956  const label nSnap,
957  const scalar minCosLayerTermination,
958 
959  const scalarField& layerThickness,
960  const scalarField& minThickness,
961 
962  const pointVectorField& dispVec,
963  const pointScalarField& medialRatio,
964  const pointScalarField& medialDist,
965 
966  List<extrudeMode>& extrudeStatus,
967  pointField& patchDisp,
968  labelList& patchNLayers
969 ) const
970 {
971  Info<< "shrinkMeshMedialDistance : Smoothing using Medial Axis ..." << endl;
972 
973  const polyMesh& mesh = meshMover.mesh();
974 
975  const indirectPrimitivePatch& pp = meshMover.patch();
976  const labelList& meshPoints = pp.meshPoints();
977 
978  // Precalulate master edge (only relevant for shared edges)
979  PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
980  // Precalculate meshEdge per pp edge
981  labelList meshEdges(pp.nEdges());
982 
983  forAll(meshEdges, patchEdgeI)
984  {
985  const edge& e = pp.edges()[patchEdgeI];
986 
987  label v0 = pp.meshPoints()[e[0]];
988  label v1 = pp.meshPoints()[e[1]];
989  meshEdges[patchEdgeI] = meshTools::findEdge
990  (
991  mesh.edges(),
992  mesh.pointEdges()[v0],
993  v0,
994  v1
995  );
996  }
997 
998 
999  scalarField thickness(layerThickness.size());
1000 
1001  thickness = mag(patchDisp);
1002 
1003  forAll(thickness, patchPointI)
1004  {
1005  if (extrudeStatus[patchPointI] == NOEXTRUDE)
1006  {
1007  thickness[patchPointI] = 0.0;
1008  }
1009  }
1010 
1011  label numThicknessRatioExclude = 0;
1012 
1013  // reduce thickness where thickness/medial axis distance large
1014  forAll(meshPoints, patchPointI)
1015  {
1016  if (extrudeStatus[patchPointI] != NOEXTRUDE)
1017  {
1018  label pointI = meshPoints[patchPointI];
1019 
1020  scalar mDist = medialDist[pointI];
1021 
1022  scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
1023 
1024  if (thicknessRatio > maxThicknessToMedialRatio)
1025  {
1026  // Truncate thickness.
1027  thickness[patchPointI] =
1028  0.5*(minThickness[patchPointI]+thickness[patchPointI]);
1029 
1030  patchDisp[patchPointI] =
1031  thickness[patchPointI]
1032  * patchDisp[patchPointI]
1033  / (mag(patchDisp[patchPointI]) + VSMALL);
1034  numThicknessRatioExclude++;
1035  }
1036  }
1037  }
1038 
1039  reduce(numThicknessRatioExclude, sumOp<label>());
1040  Info<< "shrinkMeshMedialDistance : Reduce layer thickness at "
1041  << numThicknessRatioExclude
1042  << " nodes where thickness to medial axis distance is large " << endl;
1043 
1044  // find points where layer growth isolated to a lone point, edge or face
1045 
1046  findIsolatedRegions
1047  (
1048  pp,
1049  isMasterEdge,
1050  meshEdges,
1051  minCosLayerTermination,
1052  thickness,
1053  extrudeStatus,
1054  patchDisp,
1055  patchNLayers
1056  );
1057 
1058  // smooth layer thickness on moving patch
1059  smoothField
1060  (
1061  meshMover,
1062  isMasterEdge,
1063  meshEdges,
1064  minThickness,
1065  nSmoothThickness,
1066  thickness
1067  );
1068 
1069  List<pointData> pointWallDist(mesh.nPoints());
1070 
1071  const pointField& points = mesh.points();
1072  // 1. Calculate distance to points where displacement is specified.
1073  // This wave is used to transport layer thickness
1074  {
1075  // Distance to wall and medial axis on edges.
1076  List<pointData> edgeWallDist(mesh.nEdges());
1077  labelList wallPoints(meshPoints.size());
1078 
1079  // Seed data.
1080  List<pointData> wallInfo(meshPoints.size());
1081 
1082  forAll(meshPoints, patchPointI)
1083  {
1084  label pointI = meshPoints[patchPointI];
1085  wallPoints[patchPointI] = pointI;
1086  wallInfo[patchPointI] = pointData
1087  (
1088  points[pointI],
1089  0.0,
1090  thickness[patchPointI], // transport layer thickness
1091  vector::zero // passive vector
1092  );
1093  }
1094 
1095  // Do all calculations
1096  PointEdgeWave<pointData> wallDistCalc
1097  (
1098  mesh,
1099  wallPoints,
1100  wallInfo,
1101  pointWallDist,
1102  edgeWallDist,
1103  mesh.globalData().nTotalPoints() // max iterations
1104  );
1105  }
1106 
1107  // Calculate scaled displacement vector
1108  pointVectorField& displacement = meshMover.displacement();
1109 
1110  forAll(displacement, pointI)
1111  {
1112  // 1) displacement on nearest wall point, scaled by medialRatio
1113  // (wall distance / medial distance)
1114  // 2) pointWallDist[pointI].s() is layer thickness transported
1115  // from closest wall point.
1116  // 3) shrink in opposite direction of addedPoints
1117  displacement[pointI] =
1118  -medialRatio[pointI]
1119  *pointWallDist[pointI].s()
1120  *dispVec[pointI];
1121  }
1122 
1123  // Current faces to check. Gets modified in meshMover.scaleMesh
1124  labelList checkFaces(identity(mesh.nFaces()));
1125 
1126  Info<< "shrinkMeshMedialDistance : Moving mesh ..." << endl;
1127  scalar oldErrorReduction = -1;
1128 
1129  for (label iter = 0; iter < 2*nSnap ; iter++)
1130  {
1131  Info<< "Iteration " << iter << endl;
1132  if (iter == nSnap)
1133  {
1134  Info<< "Displacement scaling for error reduction set to 0." << endl;
1135  oldErrorReduction = meshMover.setErrorReduction(0.0);
1136  }
1137 
1138  if
1139  (
1140  meshMover.scaleMesh
1141  (
1142  checkFaces,
1143  List<labelPair>(0),
1144  meshMover.paramDict(),
1145  meshQualityDict,
1146  true,
1147  nAllowableErrors
1148  )
1149  )
1150  {
1151  Info<< "shrinkMeshMedialDistance : Successfully moved mesh" << endl;
1152  break;
1153  }
1154  }
1155 
1156  if (oldErrorReduction >= 0)
1157  {
1158  meshMover.setErrorReduction(oldErrorReduction);
1159  }
1160 
1161  Info<< "shrinkMeshMedialDistance : Finished moving mesh ..." << endl;
1162 }
1163 
1164 
1165 // ************************ vim: set sw=4 sts=4 et: ************************ //