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) 2011-2013 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  Shrinking mesh (part of adding cell layers)
26 
27 \*----------------------------------------------------------------------------*/
28 
29 #include "autoLayerDriver.H"
30 #include "fvMesh.H"
31 #include "Time.H"
32 #include "pointFields.H"
33 #include "motionSmoother.H"
34 #include "pointData.H"
35 #include "PointEdgeWave.H"
36 #include "OBJstream.H"
37 #include "meshTools.H"
38 #include "PatchTools.H"
39 #include "unitConversion.H"
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 // Calculate inverse sum of edge weights (currently always 1.0)
44 void Foam::autoLayerDriver::sumWeights
45 (
46  const PackedBoolList& isMasterEdge,
47  const labelList& meshEdges,
48  const labelList& meshPoints,
49  const edgeList& edges,
50  scalarField& invSumWeight
51 ) const
52 {
53  const pointField& pts = meshRefiner_.mesh().points();
54 
55  invSumWeight = 0;
56 
57  forAll(edges, edgeI)
58  {
59  if (isMasterEdge.get(meshEdges[edgeI]) == 1)
60  {
61  const edge& e = edges[edgeI];
62  //scalar eWeight = edgeWeights[edgeI];
63  //scalar eWeight = 1.0;
64 
65  scalar eMag = max
66  (
67  VSMALL,
68  mag
69  (
70  pts[meshPoints[e[1]]]
71  - pts[meshPoints[e[0]]]
72  )
73  );
74  scalar eWeight = 1.0/eMag;
75 
76  invSumWeight[e[0]] += eWeight;
77  invSumWeight[e[1]] += eWeight;
78  }
79  }
80 
82  (
83  meshRefiner_.mesh(),
84  meshPoints,
85  invSumWeight,
86  plusEqOp<scalar>(),
87  scalar(0.0) // null value
88  );
89 
90  forAll(invSumWeight, pointI)
91  {
92  scalar w = invSumWeight[pointI];
93 
94  if (w > 0.0)
95  {
96  invSumWeight[pointI] = 1.0/w;
97  }
98  }
99 }
100 
101 
102 // Smooth field on moving patch
103 void Foam::autoLayerDriver::smoothField
104 (
105  const motionSmoother& meshMover,
106  const PackedBoolList& isMasterPoint,
107  const PackedBoolList& isMasterEdge,
108  const labelList& meshEdges,
109  const scalarField& fieldMin,
110  const label nSmoothDisp,
111  scalarField& field
112 ) const
113 {
114  const indirectPrimitivePatch& pp = meshMover.patch();
115  const edgeList& edges = pp.edges();
116  const labelList& meshPoints = pp.meshPoints();
117 
118  scalarField invSumWeight(pp.nPoints());
119  sumWeights
120  (
121  isMasterEdge,
122  meshEdges,
123  meshPoints,
124  edges,
125  invSumWeight
126  );
127 
128  // Get smoothly varying patch field.
129  Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
130 
131  for (label iter = 0; iter < nSmoothDisp; iter++)
132  {
133  scalarField average(pp.nPoints());
134  averageNeighbours
135  (
136  meshMover.mesh(),
137  isMasterEdge,
138  meshEdges,
139  meshPoints,
140  edges,
141  invSumWeight,
142  field,
143  average
144  );
145 
146  // Transfer to field
147  forAll(field, pointI)
148  {
149  //full smoothing neighbours + point value
150  average[pointI] = 0.5*(field[pointI]+average[pointI]);
151 
152  // perform monotonic smoothing
153  if
154  (
155  average[pointI] < field[pointI]
156  && average[pointI] >= fieldMin[pointI]
157  )
158  {
159  field[pointI] = average[pointI];
160  }
161  }
162 
163  // Do residual calculation every so often.
164  if ((iter % 10) == 0)
165  {
166  scalar resid = meshRefinement::gAverage
167  (
168  meshMover.mesh(),
169  isMasterPoint,
170  meshPoints,
171  mag(field-average)()
172  );
173  Info<< " Iteration " << iter << " residual " << resid << endl;
174  }
175  }
176 }
177 //XXXXXXXXX
178 //void Foam::autoLayerDriver::smoothField
179 //(
180 // const motionSmoother& meshMover,
181 // const PackedBoolList& isMasterEdge,
182 // const labelList& meshEdges,
183 // const scalarField& fieldMin,
184 // const label nSmoothDisp,
185 // scalarField& field
186 //) const
187 //{
188 // const indirectPrimitivePatch& pp = meshMover.patch();
189 // const edgeList& edges = pp.edges();
190 // const labelList& meshPoints = pp.meshPoints();
191 //
192 // scalarField invSumWeight(pp.nPoints());
193 // sumWeights
194 // (
195 // isMasterEdge,
196 // meshEdges,
197 // meshPoints,
198 // edges,
199 // invSumWeight
200 // );
201 //
202 // // Get smoothly varying patch field.
203 // Info<< "shrinkMeshDistance : (lambda-mu) Smoothing field ..." << endl;
204 //
205 //
206 // const scalar lambda = 0.33;
207 // const scalar mu = -0.34;
208 //
209 // for (label iter = 0; iter < 90; iter++)
210 // {
211 // scalarField average(pp.nPoints());
212 //
213 // // Calculate average of field
214 // averageNeighbours
215 // (
216 // meshMover.mesh(),
217 // isMasterEdge,
218 // meshEdges,
219 // meshPoints,
220 // pp.edges(),
221 // invSumWeight,
222 // field,
223 // average
224 // );
225 //
226 // forAll(field, i)
227 // {
228 // if (field[i] >= fieldMin[i])
229 // {
230 // field[i] = (1-lambda)*field[i]+lambda*average[i];
231 // }
232 // }
233 //
234 //
235 // // Calculate average of field
236 // averageNeighbours
237 // (
238 // meshMover.mesh(),
239 // isMasterEdge,
240 // meshEdges,
241 // meshPoints,
242 // pp.edges(),
243 // invSumWeight,
244 // field,
245 // average
246 // );
247 //
248 // forAll(field, i)
249 // {
250 // if (field[i] >= fieldMin[i])
251 // {
252 // field[i] = (1-mu)*field[i]+mu*average[i];
253 // }
254 // }
255 //
256 //
257 // // Do residual calculation every so often.
258 // if ((iter % 10) == 0)
259 // {
260 // Info<< " Iteration " << iter << " residual "
261 // << gSum(mag(field-average))
262 // /returnReduce(average.size(), sumOp<label>())
263 // << endl;
264 // }
265 // }
266 //}
267 //XXXXXXXXX
268 
269 // Smooth normals on moving patch.
270 void Foam::autoLayerDriver::smoothPatchNormals
271 (
272  const motionSmoother& meshMover,
273  const PackedBoolList& isMasterPoint,
274  const PackedBoolList& isMasterEdge,
275  const labelList& meshEdges,
276  const label nSmoothDisp,
277  pointField& normals
278 ) const
279 {
280  const indirectPrimitivePatch& pp = meshMover.patch();
281  const edgeList& edges = pp.edges();
282  const labelList& meshPoints = pp.meshPoints();
283 
284  // Calculate inverse sum of weights
285 
286  scalarField invSumWeight(pp.nPoints());
287  sumWeights
288  (
289  isMasterEdge,
290  meshEdges,
291  meshPoints,
292  edges,
293  invSumWeight
294  );
295 
296  // Get smoothly varying internal normals field.
297  Info<< "shrinkMeshDistance : Smoothing normals ..." << endl;
298 
299  for (label iter = 0; iter < nSmoothDisp; iter++)
300  {
301  vectorField average(pp.nPoints());
302  averageNeighbours
303  (
304  meshMover.mesh(),
305  isMasterEdge,
306  meshEdges,
307  meshPoints,
308  edges,
309  invSumWeight,
310  normals,
311  average
312  );
313 
314  // Do residual calculation every so often.
315  if ((iter % 10) == 0)
316  {
317  scalar resid = meshRefinement::gAverage
318  (
319  meshMover.mesh(),
320  isMasterPoint,
321  meshPoints,
322  mag(normals-average)()
323  );
324  Info<< " Iteration " << iter << " residual " << resid << endl;
325  }
326 
327  // Transfer to normals vector field
328  forAll(average, pointI)
329  {
330  // full smoothing neighbours + point value
331  average[pointI] = 0.5*(normals[pointI]+average[pointI]);
332  normals[pointI] = average[pointI];
333  normals[pointI] /= mag(normals[pointI]) + VSMALL;
334  }
335  }
336 }
337 
338 
339 // Smooth normals in interior.
340 void Foam::autoLayerDriver::smoothNormals
341 (
342  const label nSmoothDisp,
343  const PackedBoolList& isMasterPoint,
344  const PackedBoolList& isMasterEdge,
345  const labelList& fixedPoints,
346  pointVectorField& normals
347 ) const
348 {
349  // Get smoothly varying internal normals field.
350  Info<< "shrinkMeshDistance : Smoothing normals in interior ..." << endl;
351 
352  const fvMesh& mesh = meshRefiner_.mesh();
353  const edgeList& edges = mesh.edges();
354 
355  // Points that do not change.
356  PackedBoolList isFixedPoint(mesh.nPoints());
357 
358  // Internal points that are fixed
359  forAll(fixedPoints, i)
360  {
361  label meshPointI = fixedPoints[i];
362  isFixedPoint.set(meshPointI, 1);
363  }
364 
365  // Make sure that points that are coupled to meshPoints but not on a patch
366  // are fixed as well
367  syncTools::syncPointList(mesh, isFixedPoint, maxEqOp<unsigned int>(), 0);
368 
369 
370  // Correspondence between local edges/points and mesh edges/points
371  const labelList meshEdges(identity(mesh.nEdges()));
372  const labelList meshPoints(identity(mesh.nPoints()));
373 
374  // Calculate inverse sum of weights
375 
376  scalarField invSumWeight(mesh.nPoints(), 0);
377  sumWeights
378  (
379  isMasterEdge,
380  meshEdges,
381  meshPoints,
382  edges,
383  invSumWeight
384  );
385 
386  for (label iter = 0; iter < nSmoothDisp; iter++)
387  {
388  vectorField average(mesh.nPoints());
389  averageNeighbours
390  (
391  mesh,
392  isMasterEdge,
393  meshEdges,
394  meshPoints,
395  edges,
396  invSumWeight,
397  normals,
398  average
399  );
400 
401  // Do residual calculation every so often.
402  if ((iter % 10) == 0)
403  {
404  scalar resid = meshRefinement::gAverage
405  (
406  mesh,
407  isMasterPoint,
408  mag(normals-average)()
409  );
410  Info<< " Iteration " << iter << " residual " << resid << endl;
411  }
412 
413 
414  // Transfer to normals vector field
415  forAll(average, pointI)
416  {
417  if (isFixedPoint.get(pointI) == 0)
418  {
419  //full smoothing neighbours + point value
420  average[pointI] = 0.5*(normals[pointI]+average[pointI]);
421  normals[pointI] = average[pointI];
422  normals[pointI] /= mag(normals[pointI]) + VSMALL;
423  }
424  }
425  }
426 }
427 
428 
429 // Tries and find a medial axis point. Done by comparing vectors to nearest
430 // wall point for both vertices of edge.
431 bool Foam::autoLayerDriver::isMaxEdge
432 (
433  const List<pointData>& pointWallDist,
434  const label edgeI,
435  const scalar minCos
436 ) const
437 {
438  const fvMesh& mesh = meshRefiner_.mesh();
439  const pointField& points = mesh.points();
440 
441  // Do not mark edges with one side on moving wall.
442 
443  const edge& e = mesh.edges()[edgeI];
444 
445  vector v0(points[e[0]] - pointWallDist[e[0]].origin());
446  scalar magV0(mag(v0));
447 
448  if (magV0 < SMALL)
449  {
450  return false;
451  }
452 
453  vector v1(points[e[1]] - pointWallDist[e[1]].origin());
454  scalar magV1(mag(v1));
455 
456  if (magV1 < SMALL)
457  {
458  return false;
459  }
460 
461 
462  //- Detect based on vector to nearest point differing for both endpoints
463  //v0 /= magV0;
464  //v1 /= magV1;
465  //
467  //if ((v0 & v1) < minCos)
468  //{
469  // return true;
470  //}
471  //else
472  //{
473  // return false;
474  //}
475 
476  //- Detect based on extrusion vector differing for both endpoints
477  // the idea is that e.g. a sawtooth wall can still be extruded
478  // successfully as long as it is done all to the same direction.
479  if ((pointWallDist[e[0]].v() & pointWallDist[e[1]].v()) < minCos)
480  {
481  return true;
482  }
483  else
484  {
485  return false;
486  }
487 }
488 
489 
490 // Stop layer growth where mesh wraps around edge with a
491 // large feature angle
492 void Foam::autoLayerDriver::handleFeatureAngleLayerTerminations
493 (
494  const scalar minCos,
495  const PackedBoolList& isMasterPoint,
496  const indirectPrimitivePatch& pp,
497  const labelList& meshEdges,
498  List<extrudeMode>& extrudeStatus,
499  pointField& patchDisp,
500  labelList& patchNLayers,
501  label& nPointCounter
502 ) const
503 {
504  const fvMesh& mesh = meshRefiner_.mesh();
505 
506  // Mark faces that have all points extruded
507  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
508 
509  boolList extrudedFaces(pp.size(), true);
510 
511  forAll(pp.localFaces(), faceI)
512  {
513  const face& f = pp.localFaces()[faceI];
514 
515  forAll(f, fp)
516  {
517  if (extrudeStatus[f[fp]] == NOEXTRUDE)
518  {
519  extrudedFaces[faceI] = false;
520  break;
521  }
522  }
523  }
524 
525 
526 
527  //label nOldPointCounter = nPointCounter;
528 
529  // Detect situation where two featureedge-neighbouring faces are partly or
530  // not extruded and the edge itself is extruded. In this case unmark the
531  // edge for extrusion.
532 
533 
534  List<List<point> > edgeFaceNormals(pp.nEdges());
535  List<List<bool> > edgeFaceExtrude(pp.nEdges());
536 
537  const labelListList& edgeFaces = pp.edgeFaces();
538  const vectorField& faceNormals = pp.faceNormals();
539  const labelList& meshPoints = pp.meshPoints();
540 
541  forAll(edgeFaces, edgeI)
542  {
543  const labelList& eFaces = edgeFaces[edgeI];
544 
545  edgeFaceNormals[edgeI].setSize(eFaces.size());
546  edgeFaceExtrude[edgeI].setSize(eFaces.size());
547  forAll(eFaces, i)
548  {
549  label faceI = eFaces[i];
550  edgeFaceNormals[edgeI][i] = faceNormals[faceI];
551  edgeFaceExtrude[edgeI][i] = extrudedFaces[faceI];
552  }
553  }
554 
556  (
557  mesh,
558  meshEdges,
559  edgeFaceNormals,
560  globalMeshData::ListPlusEqOp<List<point> >(), // combine operator
561  List<point>() // null value
562  );
563 
565  (
566  mesh,
567  meshEdges,
568  edgeFaceExtrude,
569  globalMeshData::ListPlusEqOp<List<bool> >(), // combine operator
570  List<bool>() // null value
571  );
572 
573 
574  forAll(edgeFaceNormals, edgeI)
575  {
576  const List<point>& eFaceNormals = edgeFaceNormals[edgeI];
577  const List<bool>& eFaceExtrude = edgeFaceExtrude[edgeI];
578 
579  if (eFaceNormals.size() == 2)
580  {
581  const edge& e = pp.edges()[edgeI];
582  label v0 = e[0];
583  label v1 = e[1];
584 
585  if
586  (
587  extrudeStatus[v0] != NOEXTRUDE
588  || extrudeStatus[v1] != NOEXTRUDE
589  )
590  {
591  if (!eFaceExtrude[0] || !eFaceExtrude[1])
592  {
593  const vector& n0 = eFaceNormals[0];
594  const vector& n1 = eFaceNormals[1];
595 
596  if ((n0 & n1) < minCos)
597  {
598  if
599  (
600  unmarkExtrusion
601  (
602  v0,
603  patchDisp,
604  patchNLayers,
605  extrudeStatus
606  )
607  )
608  {
609  if (isMasterPoint[meshPoints[v0]])
610  {
611  nPointCounter++;
612  }
613  }
614  if
615  (
616  unmarkExtrusion
617  (
618  v1,
619  patchDisp,
620  patchNLayers,
621  extrudeStatus
622  )
623  )
624  {
625  if (isMasterPoint[meshPoints[v1]])
626  {
627  nPointCounter++;
628  }
629  }
630  }
631  }
632  }
633  }
634  }
635 
636  //Info<< "Added "
637  // << returnReduce(nPointCounter-nOldPointCounter, sumOp<label>())
638  // << " point not to extrude." << endl;
639 }
640 
641 
642 // Find isolated islands (points, edges and faces and layer terminations)
643 // in the layer mesh and stop any layer growth at these points.
644 void Foam::autoLayerDriver::findIsolatedRegions
645 (
646  const scalar minCosLayerTermination,
647  const PackedBoolList& isMasterPoint,
648  const PackedBoolList& isMasterEdge,
649  const indirectPrimitivePatch& pp,
650  const labelList& meshEdges,
651  const scalarField& minThickness,
652  List<extrudeMode>& extrudeStatus,
653  pointField& patchDisp,
654  labelList& patchNLayers
655 ) const
656 {
657  const fvMesh& mesh = meshRefiner_.mesh();
658 
659  Info<< "shrinkMeshDistance : Removing isolated regions ..." << endl;
660 
661  // Keep count of number of points unextruded
662  label nPointCounter = 0;
663 
664  while (true)
665  {
666  // Stop layer growth where mesh wraps around edge with a
667  // large feature angle
668  handleFeatureAngleLayerTerminations
669  (
670  minCosLayerTermination,
671  isMasterPoint,
672  pp,
673  meshEdges,
674 
675  extrudeStatus,
676  patchDisp,
677  patchNLayers,
678  nPointCounter
679  );
680 
681  syncPatchDisplacement
682  (
683  pp,
684  minThickness,
685  patchDisp,
686  patchNLayers,
687  extrudeStatus
688  );
689 
690 
691  // Do not extrude from point where all neighbouring
692  // faces are not grown
693  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
694  boolList extrudedFaces(pp.size(), true);
695  forAll(pp.localFaces(), faceI)
696  {
697  const face& f = pp.localFaces()[faceI];
698  forAll(f, fp)
699  {
700  if (extrudeStatus[f[fp]] == NOEXTRUDE)
701  {
702  extrudedFaces[faceI] = false;
703  break;
704  }
705  }
706  }
707 
708  const labelListList& pointFaces = pp.pointFaces();
709 
710  boolList keptPoints(pp.nPoints(), false);
711  forAll(keptPoints, patchPointI)
712  {
713  const labelList& pFaces = pointFaces[patchPointI];
714 
715  forAll(pFaces, i)
716  {
717  label faceI = pFaces[i];
718  if (extrudedFaces[faceI])
719  {
720  keptPoints[patchPointI] = true;
721  break;
722  }
723  }
724  }
725 
727  (
728  mesh,
729  pp.meshPoints(),
730  keptPoints,
731  orEqOp<bool>(),
732  false // null value
733  );
734 
735  label nChanged = 0;
736 
737  forAll(keptPoints, patchPointI)
738  {
739  if (!keptPoints[patchPointI])
740  {
741  if
742  (
743  unmarkExtrusion
744  (
745  patchPointI,
746  patchDisp,
747  patchNLayers,
748  extrudeStatus
749  )
750  )
751  {
752  nPointCounter++;
753  nChanged++;
754  }
755  }
756  }
757 
758 
759  if (returnReduce(nChanged, sumOp<label>()) == 0)
760  {
761  break;
762  }
763  }
764 
765  const edgeList& edges = pp.edges();
766 
767 
768  // Count number of mesh edges using a point
769  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
770 
771  labelList isolatedPoint(pp.nPoints(),0);
772 
773  forAll(edges, edgeI)
774  {
775  if (isMasterEdge.get(meshEdges[edgeI]) == 1)
776  {
777  const edge& e = edges[edgeI];
778 
779  label v0 = e[0];
780  label v1 = e[1];
781 
782  if (extrudeStatus[v1] != NOEXTRUDE)
783  {
784  isolatedPoint[v0] += 1;
785  }
786  if (extrudeStatus[v0] != NOEXTRUDE)
787  {
788  isolatedPoint[v1] += 1;
789  }
790  }
791  }
792 
794  (
795  mesh,
796  pp.meshPoints(),
797  isolatedPoint,
798  plusEqOp<label>(),
799  label(0) // null value
800  );
801 
802  // stop layer growth on isolated faces
803  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
804  forAll(pp, faceI)
805  {
806  const face& f = pp.localFaces()[faceI];
807  bool failed = false;
808  forAll(f, fp)
809  {
810  if (isolatedPoint[f[fp]] > 2)
811  {
812  failed = true;
813  break;
814  }
815  }
816  bool allPointsExtruded = true;
817  if (!failed)
818  {
819  forAll(f, fp)
820  {
821  if (extrudeStatus[f[fp]] == NOEXTRUDE)
822  {
823  allPointsExtruded = false;
824  break;
825  }
826  }
827 
828  if (allPointsExtruded)
829  {
830  forAll(f, fp)
831  {
832  if
833  (
834  unmarkExtrusion
835  (
836  f[fp],
837  patchDisp,
838  patchNLayers,
839  extrudeStatus
840  )
841  )
842  {
843  nPointCounter++;
844  }
845  }
846  }
847  }
848  }
849 
850  reduce(nPointCounter, sumOp<label>());
851  Info<< "Number isolated points extrusion stopped : "<< nPointCounter
852  << endl;
853 }
854 
855 
856 // Calculates medial axis fields:
857 // dispVec : normal of nearest wall. Where this normal changes direction
858 // defines the medial axis
859 // medialDist : distance to medial axis
860 // medialRatio : ratio of medial distance to wall distance.
861 // (1 at wall, 0 at medial axis)
862 void Foam::autoLayerDriver::medialAxisSmoothingInfo
863 (
864  const motionSmoother& meshMover,
865  const label nSmoothNormals,
866  const label nSmoothSurfaceNormals,
867  const scalar minMedialAxisAngleCos,
868  const scalar featureAngle,
869 
870  pointVectorField& dispVec,
871  pointScalarField& medialRatio,
872  pointScalarField& medialDist,
873  pointVectorField& medialVec
874 ) const
875 {
876 
877  Info<< "medialAxisSmoothingInfo :"
878  << " Calculate distance to Medial Axis ..." << endl;
879 
880  const polyMesh& mesh = meshMover.mesh();
881  const pointField& points = mesh.points();
882 
883  const indirectPrimitivePatch& pp = meshMover.patch();
884  const labelList& meshPoints = pp.meshPoints();
885 
886  // Predetermine mesh edges
887  // ~~~~~~~~~~~~~~~~~~~~~~~
888 
889  // Precalulate master point/edge (only relevant for shared points/edges)
890  const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
891  const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
892  // Precalculate meshEdge per pp edge
893  const labelList meshEdges
894  (
895  pp.meshEdges
896  (
897  mesh.edges(),
898  mesh.pointEdges()
899  )
900  );
901 
902 
903  // Determine pointNormal
904  // ~~~~~~~~~~~~~~~~~~~~~
905 
906  pointField pointNormals(PatchTools::pointNormals(mesh, pp));
907 
908  // pointNormals
909  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
910  {
911  pointField meshPointNormals(mesh.nPoints(), point(1, 0, 0));
912  UIndirectList<point>(meshPointNormals, meshPoints) = pointNormals;
914  (
915  "pointNormals",
916  mesh,
917  meshPointNormals
918  );
919  }
920 
921  // Smooth patch normal vectors
922  smoothPatchNormals
923  (
924  meshMover,
925  isMasterPoint,
926  isMasterEdge,
927  meshEdges,
928  nSmoothSurfaceNormals,
929  pointNormals
930  );
931 
932  // smoothed pointNormals
933  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
934  {
935  pointField meshPointNormals(mesh.nPoints(), point(1, 0, 0));
936  UIndirectList<point>(meshPointNormals, meshPoints) = pointNormals;
938  (
939  "smoothed pointNormals",
940  mesh,
941  meshPointNormals
942  );
943  }
944 
945  // Calculate distance to pp points
946  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
947 
948  // Distance to wall
949  List<pointData> pointWallDist(mesh.nPoints());
950 
951  // Dummy additional info for PointEdgeWave
952  int dummyTrackData = 0;
953 
954 
955  // 1. Calculate distance to points where displacement is specified.
956  {
957  // Seed data.
958  List<pointData> wallInfo(meshPoints.size());
959 
960  forAll(meshPoints, patchPointI)
961  {
962  label pointI = meshPoints[patchPointI];
963  wallInfo[patchPointI] = pointData
964  (
965  points[pointI],
966  0.0,
967  pointI, // passive scalar
968  pointNormals[patchPointI] // surface normals
969  );
970  }
971 
972  // Do all calculations
973  List<pointData> edgeWallDist(mesh.nEdges());
974  PointEdgeWave<pointData> wallDistCalc
975  (
976  mesh,
977  meshPoints,
978  wallInfo,
979  pointWallDist,
980  edgeWallDist,
981  mesh.globalData().nTotalPoints(), // max iterations
982  dummyTrackData
983  );
984 
985 
986  label nUnvisit = returnReduce
987  (
988  wallDistCalc.getUnsetPoints(),
989  sumOp<label>()
990  );
991 
992  if (nUnvisit > 0)
993  {
994  WarningIn("autoLayerDriver::medialAxisSmoothingInfo(..)")
995  << "Walking did not visit all points." << nl
996  << " Did not visit " << nUnvisit
997  << " out of " << mesh.globalData().nTotalPoints()
998  << " points. This is not necessarily a problem" << nl
999  << " and might be due to faceZones splitting of part"
1000  << " of the domain." << nl << endl;
1001  }
1002  }
1003 
1004 
1005  // Check sync of wall distance
1006  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
1007  {
1008  pointField origin(pointWallDist.size());
1009  scalarField distSqr(pointWallDist.size());
1010  //NA scalarField passiveS(pointWallDist.size());
1011  pointField passiveV(pointWallDist.size());
1012  forAll(pointWallDist, pointI)
1013  {
1014  origin[pointI] = pointWallDist[pointI].origin();
1015  distSqr[pointI] = pointWallDist[pointI].distSqr();
1016  //passiveS[pointI] = pointWallDist[pointI].s();
1017  passiveV[pointI] = pointWallDist[pointI].v();
1018  }
1019  meshRefinement::testSyncPointList("origin", mesh, origin);
1020  meshRefinement::testSyncPointList("distSqr", mesh, distSqr);
1021  //meshRefinement::testSyncPointList("passiveS", mesh, passiveS);
1022  meshRefinement::testSyncPointList("passiveV", mesh, passiveV);
1023  }
1024 
1025 
1026  // 2. Find points with max distance and transport information back to
1027  // wall.
1028  {
1029  List<pointData> pointMedialDist(mesh.nPoints());
1030  List<pointData> edgeMedialDist(mesh.nEdges());
1031 
1032  // Seed point data.
1033  DynamicList<pointData> maxInfo(meshPoints.size());
1034  DynamicList<label> maxPoints(meshPoints.size());
1035 
1036  // 1. Medial axis points
1037 
1038  const edgeList& edges = mesh.edges();
1039 
1040  forAll(edges, edgeI)
1041  {
1042  const edge& e = edges[edgeI];
1043 
1044  if
1045  (
1046  !pointWallDist[e[0]].valid(dummyTrackData)
1047  || !pointWallDist[e[1]].valid(dummyTrackData)
1048  )
1049  {
1050  // Unvisited point. See above about nUnvisit warning
1051  }
1052  else if (isMaxEdge(pointWallDist, edgeI, minMedialAxisAngleCos))
1053  {
1054  // Both end points of edge have very different nearest wall
1055  // point. Mark both points as medial axis points.
1056 
1057  // Approximate medial axis location on edge.
1058  //const point medialAxisPt = e.centre(points);
1059  vector eVec = e.vec(points);
1060  scalar eMag = mag(eVec);
1061  if (eMag > VSMALL)
1062  {
1063  eVec /= eMag;
1064 
1065  // Calculate distance along edge
1066  const point& p0 = points[e[0]];
1067  const point& p1 = points[e[1]];
1068  scalar dist0 = (p0-pointWallDist[e[0]].origin()) & eVec;
1069  scalar dist1 = (pointWallDist[e[1]].origin()-p1) & eVec;
1070  scalar s = 0.5*(dist1+eMag+dist0);
1071 
1072  point medialAxisPt;
1073  if (s <= dist0)
1074  {
1075  medialAxisPt = p0;
1076  }
1077  else if (s >= dist0+eMag)
1078  {
1079  medialAxisPt = p1;
1080  }
1081  else
1082  {
1083  medialAxisPt = p0+(s-dist0)*eVec;
1084  }
1085 
1086  forAll(e, ep)
1087  {
1088  label pointI = e[ep];
1089 
1090  if (!pointMedialDist[pointI].valid(dummyTrackData))
1091  {
1092  maxPoints.append(pointI);
1093  maxInfo.append
1094  (
1095  pointData
1096  (
1097  medialAxisPt, //points[pointI],
1098  magSqr(points[pointI]-medialAxisPt),//0.0,
1099  pointI, // passive data
1100  vector::zero // passive data
1101  )
1102  );
1103  pointMedialDist[pointI] = maxInfo.last();
1104  }
1105  }
1106  }
1107  }
1108  }
1109 
1110 
1111  // 2. Seed non-adapt patches
1112  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1113 
1114  labelHashSet adaptPatches(meshMover.adaptPatchIDs());
1115 
1116 
1117  forAll(patches, patchI)
1118  {
1119  const polyPatch& pp = patches[patchI];
1120  const pointPatchVectorField& pvf =
1121  meshMover.displacement().boundaryField()[patchI];
1122 
1123  if
1124  (
1125  !pp.coupled()
1126  && !isA<emptyPolyPatch>(pp)
1127  && !adaptPatches.found(patchI)
1128  )
1129  {
1130  const labelList& meshPoints = pp.meshPoints();
1131 
1132  // Check the type of the patchField. The types are
1133  // - fixedValue (0 or more layers) but the >0 layers have
1134  // already been handled in the adaptPatches loop
1135  // - constraint (but not coupled) types, e.g. symmetryPlane,
1136  // slip.
1137  if (pvf.fixesValue())
1138  {
1139  // Disable all movement on fixedValue patchFields
1140  Info<< "Inserting all points on patch " << pp.name()
1141  << endl;
1142 
1143  forAll(meshPoints, i)
1144  {
1145  label pointI = meshPoints[i];
1146  if (!pointMedialDist[pointI].valid(dummyTrackData))
1147  {
1148  maxPoints.append(pointI);
1149  maxInfo.append
1150  (
1151  pointData
1152  (
1153  points[pointI],
1154  0.0,
1155  pointI, // passive data
1156  vector::zero // passive data
1157  )
1158  );
1159  pointMedialDist[pointI] = maxInfo.last();
1160  }
1161  }
1162  }
1163  else
1164  {
1165  // Based on geometry: analyse angle w.r.t. nearest moving
1166  // point. In the pointWallDist we transported the
1167  // normal as the passive vector. Note that this points
1168  // out of the originating wall so inside of the domain
1169  // on this patch.
1170  Info<< "Inserting points on patch " << pp.name()
1171  << " if angle to nearest layer patch > "
1172  << featureAngle << " degrees." << endl;
1173 
1174  scalar featureAngleCos = Foam::cos(degToRad(featureAngle));
1175  pointField pointNormals(PatchTools::pointNormals(mesh, pp));
1176 
1177  forAll(meshPoints, i)
1178  {
1179  label pointI = meshPoints[i];
1180 
1181  if
1182  (
1183  pointWallDist[pointI].valid(dummyTrackData)
1184  && !pointMedialDist[pointI].valid(dummyTrackData)
1185  )
1186  {
1187  // Check if angle not too large.
1188  scalar cosAngle =
1189  (
1190  -pointWallDist[pointI].v()
1191  & pointNormals[i]
1192  );
1193  if (cosAngle > featureAngleCos)
1194  {
1195  // Extrusion direction practically perpendicular
1196  // to the patch. Disable movement at the patch.
1197 
1198  maxPoints.append(pointI);
1199  maxInfo.append
1200  (
1201  pointData
1202  (
1203  points[pointI],
1204  0.0,
1205  pointI, // passive data
1206  vector::zero // passive data
1207  )
1208  );
1209  pointMedialDist[pointI] = maxInfo.last();
1210  }
1211  else
1212  {
1213  // Extrusion direction makes angle with patch
1214  // so allow sliding - don't insert zero points
1215  }
1216  }
1217  }
1218  }
1219  }
1220  }
1221 
1222  maxInfo.shrink();
1223  maxPoints.shrink();
1224 
1225  // Do all calculations
1226  PointEdgeWave<pointData> medialDistCalc
1227  (
1228  mesh,
1229  maxPoints,
1230  maxInfo,
1231 
1232  pointMedialDist,
1233  edgeMedialDist,
1234  mesh.globalData().nTotalPoints(), // max iterations
1235  dummyTrackData
1236  );
1237 
1238  // Extract medial axis distance as pointScalarField
1239  forAll(pointMedialDist, pointI)
1240  {
1241  medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr());
1242  medialVec[pointI] = pointMedialDist[pointI].origin();
1243  }
1244 
1245  // Check
1246  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
1247  {
1248  meshRefinement::testSyncPointList("medialDist", mesh, medialDist);
1249  meshRefinement::testSyncPointList("medialVec", mesh, medialVec);
1250  }
1251  }
1252 
1253  // Extract transported surface normals as pointVectorField
1254  forAll(dispVec, i)
1255  {
1256  if (!pointWallDist[i].valid(dummyTrackData))
1257  {
1258  dispVec[i] = vector(1, 0, 0);
1259  }
1260  else
1261  {
1262  dispVec[i] = pointWallDist[i].v();
1263  }
1264  }
1265 
1266  // Smooth normal vectors. Do not change normals on pp.meshPoints
1267  smoothNormals
1268  (
1269  nSmoothNormals,
1270  isMasterPoint,
1271  isMasterEdge,
1272  meshPoints,
1273  dispVec
1274  );
1275 
1276  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
1277  {
1278  meshRefinement::testSyncPointList("smoothed dispVec", mesh, dispVec);
1279  }
1280 
1281  // Calculate ratio point medial distance to point wall distance
1282  forAll(medialRatio, pointI)
1283  {
1284  if (!pointWallDist[pointI].valid(dummyTrackData))
1285  {
1286  medialRatio[pointI] = 0.0;
1287  }
1288  else
1289  {
1290  scalar wDist2 = pointWallDist[pointI].distSqr();
1291  scalar mDist = medialDist[pointI];
1292 
1293  if (wDist2 < sqr(SMALL) && mDist < SMALL)
1294  //- Note: maybe less strict:
1295  //(
1296  // wDist2 < sqr(meshRefiner_.mergeDistance())
1297  // && mDist < meshRefiner_.mergeDistance()
1298  //)
1299  {
1300  medialRatio[pointI] = 0.0;
1301  }
1302  else
1303  {
1304  medialRatio[pointI] = mDist / (Foam::sqrt(wDist2) + mDist);
1305  }
1306  }
1307  }
1308 
1309 
1310  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
1311  {
1312  // medialRatio
1313  meshRefinement::testSyncPointList("medialRatio", mesh, medialRatio);
1314  }
1315 
1316 
1317  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
1318  {
1319  Info<< "medialAxisSmoothingInfo :"
1320  << " Writing:" << nl
1321  << " " << dispVec.name()
1322  << " : normalised direction of nearest displacement" << nl
1323  << " " << medialDist.name()
1324  << " : distance to medial axis" << nl
1325  << " " << medialVec.name()
1326  << " : nearest point on medial axis" << nl
1327  << " " << medialRatio.name()
1328  << " : ratio of medial distance to wall distance" << nl
1329  << endl;
1330  meshRefiner_.mesh().setInstance(meshRefiner_.timeName());
1331  meshRefiner_.write
1332  (
1335  (
1338  ),
1339  mesh.time().path()/meshRefiner_.timeName()
1340  );
1341  dispVec.write();
1342  medialDist.write();
1343  medialVec.write();
1344  medialRatio.write();
1345  }
1346 }
1347 
1348 
1349 void Foam::autoLayerDriver::shrinkMeshMedialDistance
1350 (
1351  motionSmoother& meshMover,
1352  const dictionary& meshQualityDict,
1353  const List<labelPair>& baffles,
1354  const label nSmoothPatchThickness,
1355  const label nSmoothDisplacement,
1356  const scalar maxThicknessToMedialRatio,
1357  const label nAllowableErrors,
1358  const label nSnap,
1359  const scalar minCosLayerTermination,
1360 
1361  const scalarField& layerThickness,
1362  const scalarField& minThickness,
1363 
1364  const pointVectorField& dispVec,
1365  const pointScalarField& medialRatio,
1366  const pointScalarField& medialDist,
1367  const pointVectorField& medialVec,
1368 
1369  List<extrudeMode>& extrudeStatus,
1370  pointField& patchDisp,
1371  labelList& patchNLayers
1372 ) const
1373 {
1374  Info<< "shrinkMeshMedialDistance : Smoothing using Medial Axis ..." << endl;
1375 
1376  const polyMesh& mesh = meshMover.mesh();
1377 
1378  const indirectPrimitivePatch& pp = meshMover.patch();
1379  const labelList& meshPoints = pp.meshPoints();
1380 
1381  // Precalulate master points/edge (only relevant for shared points/edges)
1382  const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
1383  const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
1384  // Precalculate meshEdge per pp edge
1385  const labelList meshEdges
1386  (
1387  pp.meshEdges
1388  (
1389  mesh.edges(),
1390  mesh.pointEdges()
1391  )
1392  );
1393 
1394 
1395  scalarField thickness(layerThickness.size());
1396 
1397  thickness = mag(patchDisp);
1398 
1399  forAll(thickness, patchPointI)
1400  {
1401  if (extrudeStatus[patchPointI] == NOEXTRUDE)
1402  {
1403  thickness[patchPointI] = 0.0;
1404  }
1405  }
1406 
1407  label numThicknessRatioExclude = 0;
1408 
1409  // reduce thickness where thickness/medial axis distance large
1410  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1411 
1412  autoPtr<OBJstream> str;
1413  if (debug)
1414  {
1415  str.reset
1416  (
1417  new OBJstream
1418  (
1419  mesh.time().path()
1420  / "thicknessRatioExcludePoints_"
1421  + meshRefiner_.timeName()
1422  + ".obj"
1423  )
1424  );
1425  Info<< "Writing points with too large an extrusion distance to "
1426  << str().name() << endl;
1427  }
1428 
1429  autoPtr<OBJstream> medialVecStr;
1430  if (debug)
1431  {
1432  medialVecStr.reset
1433  (
1434  new OBJstream
1435  (
1436  mesh.time().path()
1437  / "thicknessRatioExcludeMedialVec_"
1438  + meshRefiner_.timeName()
1439  + ".obj"
1440  )
1441  );
1442  Info<< "Writing points with too large an extrusion distance to "
1443  << medialVecStr().name() << endl;
1444  }
1445 
1446  forAll(meshPoints, patchPointI)
1447  {
1448  if (extrudeStatus[patchPointI] != NOEXTRUDE)
1449  {
1450  label pointI = meshPoints[patchPointI];
1451 
1452  //- Option 1: look only at extrusion thickness v.s. distance
1453  // to nearest (medial axis or static) point.
1454  scalar mDist = medialDist[pointI];
1455  scalar thicknessRatio = thickness[patchPointI]/(mDist+VSMALL);
1456 
1457  //- Option 2: Look at component in the direction
1458  // of nearest (medial axis or static) point.
1459  vector n =
1460  patchDisp[patchPointI]
1461  / (mag(patchDisp[patchPointI]) + VSMALL);
1462  vector mVec = mesh.points()[pointI]-medialVec[pointI];
1463  mVec /= mag(mVec)+VSMALL;
1464  thicknessRatio *= (n&mVec);
1465 
1466  if (thicknessRatio > maxThicknessToMedialRatio)
1467  {
1468  // Truncate thickness.
1469  if (debug)
1470  {
1471  Pout<< "truncating displacement at "
1472  << mesh.points()[pointI]
1473  << " from " << thickness[patchPointI]
1474  << " to "
1475  << 0.5
1476  *(
1477  minThickness[patchPointI]
1478  +thickness[patchPointI]
1479  )
1480  << " medial direction:" << mVec
1481  << " extrusion direction:" << n
1482  << " with thicknessRatio:" << thicknessRatio
1483  << endl;
1484  }
1485 
1486  thickness[patchPointI] =
1487  0.5*(minThickness[patchPointI]+thickness[patchPointI]);
1488 
1489  patchDisp[patchPointI] = thickness[patchPointI]*n;
1490 
1491  if (isMasterPoint[pointI])
1492  {
1493  numThicknessRatioExclude++;
1494  }
1495 
1496  if (str.valid())
1497  {
1498  const point& pt = mesh.points()[pointI];
1499  str().write(linePointRef(pt, pt+patchDisp[patchPointI]));
1500  }
1501  if (medialVecStr.valid())
1502  {
1503  const point& pt = mesh.points()[pointI];
1504  medialVecStr().write
1505  (
1506  linePointRef
1507  (
1508  pt,
1509  medialVec[pointI]
1510  )
1511  );
1512  }
1513  }
1514  }
1515  }
1516 
1517  reduce(numThicknessRatioExclude, sumOp<label>());
1518  Info<< "shrinkMeshMedialDistance : Reduce layer thickness at "
1519  << numThicknessRatioExclude
1520  << " nodes where thickness to medial axis distance is large " << endl;
1521 
1522 
1523  // find points where layer growth isolated to a lone point, edge or face
1524 
1525  findIsolatedRegions
1526  (
1527  minCosLayerTermination,
1528  isMasterPoint,
1529  isMasterEdge,
1530  pp,
1531  meshEdges,
1532  minThickness,
1533 
1534  extrudeStatus,
1535  patchDisp,
1536  patchNLayers
1537  );
1538 
1539  // Update thickess for changed extrusion
1540  forAll(thickness, patchPointI)
1541  {
1542  if (extrudeStatus[patchPointI] == NOEXTRUDE)
1543  {
1544  thickness[patchPointI] = 0.0;
1545  }
1546  }
1547 
1548 
1549  // smooth layer thickness on moving patch
1550  smoothField
1551  (
1552  meshMover,
1553  isMasterPoint,
1554  isMasterEdge,
1555  meshEdges,
1556  minThickness,
1557  nSmoothPatchThickness,
1558 
1559  thickness
1560  );
1561 
1562  // Dummy additional info for PointEdgeWave
1563  int dummyTrackData = 0;
1564 
1565  List<pointData> pointWallDist(mesh.nPoints());
1566 
1567  const pointField& points = mesh.points();
1568  // 1. Calculate distance to points where displacement is specified.
1569  // This wave is used to transport layer thickness
1570  {
1571  // Distance to wall and medial axis on edges.
1572  List<pointData> edgeWallDist(mesh.nEdges());
1573  labelList wallPoints(meshPoints.size());
1574 
1575  // Seed data.
1576  List<pointData> wallInfo(meshPoints.size());
1577 
1578  forAll(meshPoints, patchPointI)
1579  {
1580  label pointI = meshPoints[patchPointI];
1581  wallPoints[patchPointI] = pointI;
1582  wallInfo[patchPointI] = pointData
1583  (
1584  points[pointI],
1585  0.0,
1586  thickness[patchPointI], // transport layer thickness
1587  vector::zero // passive vector
1588  );
1589  }
1590 
1591  // Do all calculations
1592  PointEdgeWave<pointData> wallDistCalc
1593  (
1594  mesh,
1595  wallPoints,
1596  wallInfo,
1597  pointWallDist,
1598  edgeWallDist,
1599  mesh.globalData().nTotalPoints(), // max iterations
1600  dummyTrackData
1601  );
1602  }
1603 
1604  // Calculate scaled displacement vector
1605  pointVectorField& displacement = meshMover.displacement();
1606 
1607  forAll(displacement, pointI)
1608  {
1609  if (!pointWallDist[pointI].valid(dummyTrackData))
1610  {
1611  displacement[pointI] = vector::zero;
1612  }
1613  else
1614  {
1615  // 1) displacement on nearest wall point, scaled by medialRatio
1616  // (wall distance / medial distance)
1617  // 2) pointWallDist[pointI].s() is layer thickness transported
1618  // from closest wall point.
1619  // 3) shrink in opposite direction of addedPoints
1620  displacement[pointI] =
1621  -medialRatio[pointI]
1622  *pointWallDist[pointI].s()
1623  *dispVec[pointI];
1624  }
1625  }
1626 
1627 
1628 
1629  // Smear displacement away from fixed values (medialRatio=0 or 1)
1630  if (nSmoothDisplacement > 0)
1631  {
1632  scalarField invSumWeight(mesh.nPoints());
1633  sumWeights
1634  (
1635  isMasterEdge,
1636  identity(mesh.nEdges()),
1637  identity(mesh.nPoints()),
1638  mesh.edges(),
1639  invSumWeight
1640  );
1641 
1642 
1643  // Get smoothly varying patch field.
1644  Info<< "shrinkMeshDistance : Smoothing displacement ..." << endl;
1645 
1646  const scalar lambda = 0.33;
1647  const scalar mu = -0.34;
1648 
1649  pointField average(mesh.nPoints());
1650  for (label iter = 0; iter < nSmoothDisplacement; iter++)
1651  {
1652  // Calculate average of field
1653  averageNeighbours
1654  (
1655  mesh,
1656  isMasterEdge,
1657  identity(mesh.nEdges()), //meshEdges,
1658  identity(mesh.nPoints()), //meshPoints,
1659  mesh.edges(), //edges,
1660  invSumWeight,
1661  displacement,
1662  average
1663  );
1664 
1665  forAll(displacement, i)
1666  {
1667  if (medialRatio[i] > SMALL && medialRatio[i] < 1-SMALL)
1668  {
1669  displacement[i] =
1670  (1-lambda)*displacement[i]
1671  +lambda*average[i];
1672  }
1673  }
1674 
1675 
1676  // Calculate average of field
1677  averageNeighbours
1678  (
1679  mesh,
1680  isMasterEdge,
1681  identity(mesh.nEdges()), //meshEdges,
1682  identity(mesh.nPoints()), //meshPoints,
1683  mesh.edges(), //edges,
1684  invSumWeight,
1685  displacement,
1686  average
1687  );
1688 
1689  forAll(displacement, i)
1690  {
1691  if (medialRatio[i] > SMALL && medialRatio[i] < 1-SMALL)
1692  {
1693  displacement[i] = (1-mu)*displacement[i]+mu*average[i];
1694  }
1695  }
1696 
1697 
1698  // Do residual calculation every so often.
1699  if ((iter % 10) == 0)
1700  {
1701  scalar resid = meshRefinement::gAverage
1702  (
1703  mesh,
1705  mag(displacement-average)()
1706  );
1707  Info<< " Iteration " << iter << " residual " << resid
1708  << endl;
1709  }
1710  }
1711  }
1712 
1713  // Make sure displacement boundary conditions is uptodate with
1714  // internal field
1715  meshMover.setDisplacementPatchFields();
1716 
1717 
1718  // Check a bit the sync of displacements
1719  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
1720  {
1721  // initial mesh
1722  meshRefinement::testSyncPointList("mesh.points()", mesh, mesh.points());
1723 
1724  // pointWallDist
1725  scalarField pWallDist(pointWallDist.size());
1726  forAll(pointWallDist, pointI)
1727  {
1728  pWallDist[pointI] = pointWallDist[pointI].s();
1729  }
1730  meshRefinement::testSyncPointList("pointWallDist", mesh, pWallDist);
1731 
1732  // dispVec
1733  meshRefinement::testSyncPointList("dispVec", mesh, dispVec);
1734 
1735  // displacement before and after correction
1737  (
1738  "displacement BEFORE",
1739  mesh,
1740  displacement
1741  );
1742 
1743  meshMover.correctBoundaryConditions(displacement);
1745  (
1746  "displacement AFTER",
1747  mesh,
1748  displacement
1749  );
1750  }
1751 
1752 
1753  if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
1754  {
1755  const_cast<Time&>(mesh.time())++;
1756  Info<< "Writing wanted-displacement mesh (possibly illegal) to "
1757  << meshRefiner_.timeName() << endl;
1758  pointField oldPoints(mesh.points());
1759 
1760  meshRefiner_.mesh().movePoints(meshMover.curPoints());
1761  // Warn meshMover for changed geometry
1762  meshMover.movePoints();
1763 
1764  // Above move will have changed the instance only on the points (which
1765  // is correct).
1766  // However the previous mesh written will be the mesh with layers
1767  // (see autoLayerDriver.C) so we now have to force writing all files
1768  // so we can easily step through time steps. Note that if you
1769  // don't write the mesh with layers this is not necessary.
1770  meshRefiner_.mesh().setInstance(meshRefiner_.timeName());
1771 
1772  meshRefiner_.write
1773  (
1776  (
1779  ),
1780  mesh.time().path()/meshRefiner_.timeName()
1781  );
1782  dispVec.write();
1783  medialDist.write();
1784  medialRatio.write();
1785 
1786  // Move mesh back
1787  meshRefiner_.mesh().movePoints(oldPoints);
1788  // Warn meshMover for changed geometry
1789  meshMover.movePoints();
1790  }
1791 
1792 
1793  // Current faces to check. Gets modified in meshMover.scaleMesh
1794  labelList checkFaces(identity(mesh.nFaces()));
1795 
1796  Info<< "shrinkMeshMedialDistance : Moving mesh ..." << endl;
1797  scalar oldErrorReduction = -1;
1798 
1799  for (label iter = 0; iter < 2*nSnap ; iter++)
1800  {
1801  Info<< "Iteration " << iter << endl;
1802  if (iter == nSnap)
1803  {
1804  Info<< "Displacement scaling for error reduction set to 0." << endl;
1805  oldErrorReduction = meshMover.setErrorReduction(0.0);
1806  }
1807 
1808  if
1809  (
1810  meshMover.scaleMesh
1811  (
1812  checkFaces,
1813  baffles,
1814  meshMover.paramDict(),
1815  meshQualityDict,
1816  true,
1817  nAllowableErrors
1818  )
1819  )
1820  {
1821  Info<< "shrinkMeshMedialDistance : Successfully moved mesh" << endl;
1822  break;
1823  }
1824  }
1825 
1826  if (oldErrorReduction >= 0)
1827  {
1828  meshMover.setErrorReduction(oldErrorReduction);
1829  }
1830 
1831  Info<< "shrinkMeshMedialDistance : Finished moving mesh ..." << endl;
1832 }
1833 
1834 
1835 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
const dimensionedScalar mu
Atomic mass unit.
Definition: createFields.H:13
bool write() const
Write mesh and all data.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
vector point
Point is a vector.
Definition: point.H:41
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: fvMesh.C:728
static T gAverage(const PackedBoolList &isMasterElem, const UList< T > &values)
Helper: calculate average.
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
static PackedBoolList getMasterEdges(const polyMesh &)
Get per edge whether it is uncoupled or a master of a.
Definition: syncTools.C:112
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
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
label n
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
#define WarningIn(functionName)
Report a warning using Foam::Warning.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
static void testSyncPointList(const string &msg, const polyMesh &mesh, const List< scalar > &fld)
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
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
dimensionedScalar cos(const dimensionedScalar &ds)
Unit conversion functions.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
List< edge > edgeList
Definition: edgeList.H:38
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 writeType writeLevel()
Get/set write level.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const Vector zero
Definition: Vector.H:80
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
Definition: pointFields.H:50
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
Definition: pointFields.H:49
dimensionedSymmTensor sqr(const dimensionedVector &dv)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
volScalarField scalarField(fieldObject, mesh)
const fvMesh & mesh() const
Reference to mesh.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
pointPatchField< vector > pointPatchVectorField