medialAxisMeshMover.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) 2014-2016 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 \*---------------------------------------------------------------------------*/
25 
26 #include "medialAxisMeshMover.H"
28 #include "pointFields.H"
29 #include "valuePointPatchFields.H"
30 #include "PointEdgeWave.H"
31 #include "meshRefinement.H"
32 #include "unitConversion.H"
33 #include "PatchTools.H"
34 #include "OBJstream.H"
35 #include "pointData.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42  defineTypeNameAndDebug(medialAxisMeshMover, 0);
43 
45  (
46  externalDisplacementMeshMover,
47  medialAxisMeshMover,
48  dictionary
49  );
50 }
51 
52 
53 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
54 
55 Foam::labelList Foam::medialAxisMeshMover::getFixedValueBCs
56 (
57  const pointVectorField& fld
58 )
59 {
60  DynamicList<label> adaptPatchIDs;
61  forAll(fld.boundaryField(), patchi)
62  {
63  const pointPatchField<vector>& patchFld =
64  fld.boundaryField()[patchi];
65 
66  if (isA<valuePointPatchField<vector>>(patchFld))
67  {
68  if (isA<zeroFixedValuePointPatchField<vector>>(patchFld))
69  {
70  // Special condition of fixed boundary condition. Does not
71  // get adapted
72  }
73  else
74  {
75  adaptPatchIDs.append(patchi);
76  }
77  }
78  }
79  return adaptPatchIDs;
80 }
81 
82 
84 Foam::medialAxisMeshMover::getPatch
85 (
86  const polyMesh& mesh,
87  const labelList& patchIDs
88 )
89 {
90  const polyBoundaryMesh& patches = mesh.boundaryMesh();
91 
92  // Count faces.
93  label nFaces = 0;
94 
95  forAll(patchIDs, i)
96  {
97  const polyPatch& pp = patches[patchIDs[i]];
98 
99  nFaces += pp.size();
100  }
101 
102  // Collect faces.
103  labelList addressing(nFaces);
104  nFaces = 0;
105 
106  forAll(patchIDs, i)
107  {
108  const polyPatch& pp = patches[patchIDs[i]];
109 
110  label meshFacei = pp.start();
111 
112  forAll(pp, i)
113  {
114  addressing[nFaces++] = meshFacei++;
115  }
116  }
117 
118  return autoPtr<indirectPrimitivePatch>
119  (
121  (
122  IndirectList<face>(mesh.faces(), addressing),
123  mesh.points()
124  )
125  );
126 }
127 
128 
129 void Foam::medialAxisMeshMover::smoothPatchNormals
130 (
131  const label nSmoothDisp,
132  const PackedBoolList& isPatchMasterPoint,
133  const PackedBoolList& isPatchMasterEdge,
134  pointField& normals
135 ) const
136 {
137  const indirectPrimitivePatch& pp = adaptPatchPtr_();
138  const edgeList& edges = pp.edges();
139  const labelList& meshPoints = pp.meshPoints();
140 
141  // Get smoothly varying internal normals field.
142  Info<< typeName << " : Smoothing normals ..." << endl;
143 
144  scalarField edgeWeights(edges.size());
145  scalarField invSumWeight(meshPoints.size());
147  (
148  mesh(),
149  isPatchMasterEdge,
150  meshPoints,
151  edges,
152  edgeWeights,
153  invSumWeight
154  );
155 
156 
158  for (label iter = 0; iter < nSmoothDisp; iter++)
159  {
161  (
162  mesh(),
163  isPatchMasterEdge,
164  meshPoints,
165  edges,
166  edgeWeights,
167  normals,
168  average
169  );
170  average *= invSumWeight;
171 
172  // Do residual calculation every so often.
173  if ((iter % 10) == 0)
174  {
175  scalar resid = meshRefinement::gAverage
176  (
177  isPatchMasterPoint,
178  mag(normals-average)()
179  );
180  Info<< " Iteration " << iter << " residual " << resid << endl;
181  }
182 
183  // Transfer to normals vector field
184  forAll(average, pointi)
185  {
186  // full smoothing neighbours + point value
187  average[pointi] = 0.5*(normals[pointi]+average[pointi]);
188  normals[pointi] = average[pointi];
189  normals[pointi] /= mag(normals[pointi]) + VSMALL;
190  }
191  }
192 }
193 
194 
195 // Smooth normals in interior.
196 void Foam::medialAxisMeshMover::smoothNormals
197 (
198  const label nSmoothDisp,
199  const PackedBoolList& isMeshMasterPoint,
200  const PackedBoolList& isMeshMasterEdge,
201  const labelList& fixedPoints,
202  pointVectorField& normals
203 ) const
204 {
205  // Get smoothly varying internal normals field.
206  Info<< typeName
207  << " : Smoothing normals in interior ..." << endl;
208 
209  const edgeList& edges = mesh().edges();
210 
211  // Points that do not change.
212  PackedBoolList isFixedPoint(mesh().nPoints());
213 
214  // Internal points that are fixed
215  forAll(fixedPoints, i)
216  {
217  label meshPointi = fixedPoints[i];
218  isFixedPoint.set(meshPointi, 1);
219  }
220 
221  // Make sure that points that are coupled to meshPoints but not on a patch
222  // are fixed as well
223  syncTools::syncPointList(mesh(), isFixedPoint, maxEqOp<unsigned int>(), 0);
224 
225 
226  // Correspondence between local edges/points and mesh edges/points
227  const labelList meshPoints(identity(mesh().nPoints()));
228 
229  // Calculate inverse sum of weights
230 
231  scalarField edgeWeights(mesh().nEdges());
232  scalarField invSumWeight(meshPoints.size());
234  (
235  mesh(),
236  isMeshMasterEdge,
237  meshPoints,
238  edges,
239  edgeWeights,
240  invSumWeight
241  );
242 
244  for (label iter = 0; iter < nSmoothDisp; iter++)
245  {
247  (
248  mesh(),
249  isMeshMasterEdge,
250  meshPoints,
251  edges,
252  edgeWeights,
253  normals,
254  average
255  );
256  average *= invSumWeight;
257 
258  // Do residual calculation every so often.
259  if ((iter % 10) == 0)
260  {
261  scalar resid = meshRefinement::gAverage
262  (
263  isMeshMasterPoint,
264  mag(normals-average)()
265  );
266  Info<< " Iteration " << iter << " residual " << resid << endl;
267  }
268 
269 
270  // Transfer to normals vector field
271  forAll(average, pointi)
272  {
273  if (isFixedPoint.get(pointi) == 0)
274  {
275  //full smoothing neighbours + point value
276  average[pointi] = 0.5*(normals[pointi]+average[pointi]);
277  normals[pointi] = average[pointi];
278  normals[pointi] /= mag(normals[pointi]) + VSMALL;
279  }
280  }
281  }
282 }
283 
284 
285 // Tries and find a medial axis point. Done by comparing vectors to nearest
286 // wall point for both vertices of edge.
287 bool Foam::medialAxisMeshMover::isMaxEdge
288 (
289  const List<pointData>& pointWallDist,
290  const label edgeI,
291  const scalar minCos
292 ) const
293 {
294  const pointField& points = mesh().points();
295 
296  // Do not mark edges with one side on moving wall.
297 
298  const edge& e = mesh().edges()[edgeI];
299 
300  vector v0(points[e[0]] - pointWallDist[e[0]].origin());
301  scalar magV0(mag(v0));
302 
303  if (magV0 < SMALL)
304  {
305  return false;
306  }
307 
308  vector v1(points[e[1]] - pointWallDist[e[1]].origin());
309  scalar magV1(mag(v1));
310 
311  if (magV1 < SMALL)
312  {
313  return false;
314  }
315 
316 
317  //- Detect based on vector to nearest point differing for both endpoints
318  //v0 /= magV0;
319  //v1 /= magV1;
320  //
322  //if ((v0 & v1) < minCos)
323  //{
324  // return true;
325  //}
326  //else
327  //{
328  // return false;
329  //}
330 
331  //- Detect based on extrusion vector differing for both endpoints
332  // the idea is that e.g. a sawtooth wall can still be extruded
333  // successfully as long as it is done all to the same direction.
334  if ((pointWallDist[e[0]].v() & pointWallDist[e[1]].v()) < minCos)
335  {
336  return true;
337  }
338  else
339  {
340  return false;
341  }
342 }
343 
344 
345 void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
346 {
347  Info<< typeName
348  << " : Calculating distance to Medial Axis ..." << endl;
349 
350  const pointField& points = mesh().points();
351 
352  const indirectPrimitivePatch& pp = adaptPatchPtr_();
353  const labelList& meshPoints = pp.meshPoints();
354 
355 
356  // Read a few parameters
357  // ~~~~~~~~~~~~~~~~~~~~~
358 
359  //- Smooth surface normals
360  const label nSmoothSurfaceNormals = readLabel
361  (
362  coeffDict.lookup("nSmoothSurfaceNormals")
363  );
364 
365  //- When is medial axis
366  word angleKey = "minMedialAxisAngle";
367  if (!coeffDict.found(angleKey))
368  {
369  // Backwards compatibility
370  angleKey = "minMedianAxisAngle";
371  }
372  scalar minMedialAxisAngleCos = Foam::cos
373  (
374  degToRad(readScalar(coeffDict.lookup(angleKey)))
375  );
376 
377  //- Feature angle when to stop adding layers
378  const scalar featureAngle = readScalar(coeffDict.lookup("featureAngle"));
379 
380  //- When to slip along wall
381  const scalar slipFeatureAngle =
382  (
383  coeffDict.found("slipFeatureAngle")
384  ? readScalar(coeffDict.lookup("slipFeatureAngle"))
385  : 0.5*featureAngle
386  );
387 
388  //- Smooth internal normals
389  const label nSmoothNormals = readLabel
390  (
391  coeffDict.lookup("nSmoothNormals")
392  );
393 
394  //- Number of edges walking out
395  const label nMedialAxisIter = coeffDict.lookupOrDefault<label>
396  (
397  "nMedialAxisIter",
399  );
400 
401 
402  // Predetermine mesh edges
403  // ~~~~~~~~~~~~~~~~~~~~~~~
404 
405  // Precalulate (mesh) master point/edge (only relevant for shared pts/edges)
406  const PackedBoolList isMeshMasterPoint(syncTools::getMasterPoints(mesh()));
407  const PackedBoolList isMeshMasterEdge(syncTools::getMasterEdges(mesh()));
408  // Precalculate meshEdge per pp edge
409  const labelList meshEdges
410  (
411  pp.meshEdges
412  (
413  mesh().edges(),
414  mesh().pointEdges()
415  )
416  );
417 
418  // Precalulate (patch) master point/edge
419  const PackedBoolList isPatchMasterPoint
420  (
422  (
423  mesh(),
424  meshPoints
425  )
426  );
427  const PackedBoolList isPatchMasterEdge
428  (
430  (
431  mesh(),
432  meshEdges
433  )
434  );
435 
436  // Determine pointNormal
437  // ~~~~~~~~~~~~~~~~~~~~~
438 
439  pointField pointNormals(PatchTools::pointNormals(mesh(), pp));
440 
441  // Smooth patch normal vectors
442  smoothPatchNormals
443  (
444  nSmoothSurfaceNormals,
445  isPatchMasterPoint,
446  isPatchMasterEdge,
447  pointNormals
448  );
449 
450 
451  // Calculate distance to pp points
452  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
453 
454  // Distance to wall
455  List<pointData> pointWallDist(mesh().nPoints());
456 
457  // Dummy additional info for PointEdgeWave
458  int dummyTrackData = 0;
459 
460 
461  // 1. Calculate distance to points where displacement is specified.
462  {
463  // Seed data.
464  List<pointData> wallInfo(meshPoints.size());
465 
466  forAll(meshPoints, patchPointi)
467  {
468  label pointi = meshPoints[patchPointi];
469  wallInfo[patchPointi] = pointData
470  (
471  points[pointi],
472  0.0,
473  pointi, // passive scalar
474  pointNormals[patchPointi] // surface normals
475  );
476  }
477 
478  // Do all calculations
479  List<pointData> edgeWallDist(mesh().nEdges());
480  PointEdgeWave<pointData> wallDistCalc
481  (
482  mesh(),
483  meshPoints,
484  wallInfo,
485  pointWallDist,
486  edgeWallDist,
487  0, // max iterations
488  dummyTrackData
489  );
490  wallDistCalc.iterate(nMedialAxisIter);
491 
492  label nUnvisit = returnReduce
493  (
494  wallDistCalc.getUnsetPoints(),
495  sumOp<label>()
496  );
497 
498  if (nUnvisit > 0)
499  {
500  if (nMedialAxisIter > 0)
501  {
502  Info<< typeName
503  << " : Limited walk to " << nMedialAxisIter
504  << " steps. Not visited " << nUnvisit
505  << " out of " << mesh().globalData().nTotalPoints()
506  << " points" << endl;
507  }
508  else
509  {
511  << "Walking did not visit all points." << nl
512  << " Did not visit " << nUnvisit
513  << " out of " << mesh().globalData().nTotalPoints()
514  << " points. This is not necessarily a problem" << nl
515  << " and might be due to faceZones splitting of part"
516  << " of the domain." << nl << endl;
517  }
518  }
519  }
520 
521 
522  // 2. Find points with max distance and transport information back to
523  // wall.
524  {
525  List<pointData> pointMedialDist(mesh().nPoints());
526  List<pointData> edgeMedialDist(mesh().nEdges());
527 
528  // Seed point data.
529  DynamicList<pointData> maxInfo(meshPoints.size());
530  DynamicList<label> maxPoints(meshPoints.size());
531 
532  // 1. Medial axis points
533 
534  const edgeList& edges = mesh().edges();
535 
536  forAll(edges, edgeI)
537  {
538  const edge& e = edges[edgeI];
539 
540  if
541  (
542  !pointWallDist[e[0]].valid(dummyTrackData)
543  || !pointWallDist[e[1]].valid(dummyTrackData)
544  )
545  {
546  // Unvisited point. See above about nUnvisit warning
547  forAll(e, ep)
548  {
549  label pointi = e[ep];
550 
551  if (!pointMedialDist[pointi].valid(dummyTrackData))
552  {
553  maxPoints.append(pointi);
554  maxInfo.append
555  (
556  pointData
557  (
558  points[pointi],
559  0.0,
560  pointi, // passive data
561  Zero // passive data
562  )
563  );
564  pointMedialDist[pointi] = maxInfo.last();
565  }
566  }
567 
568  }
569  else if (isMaxEdge(pointWallDist, edgeI, minMedialAxisAngleCos))
570  {
571  // Both end points of edge have very different nearest wall
572  // point. Mark both points as medial axis points.
573 
574  // Approximate medial axis location on edge.
575  //const point medialAxisPt = e.centre(points);
576  vector eVec = e.vec(points);
577  scalar eMag = mag(eVec);
578  if (eMag > VSMALL)
579  {
580  eVec /= eMag;
581 
582  // Calculate distance along edge
583  const point& p0 = points[e[0]];
584  const point& p1 = points[e[1]];
585  scalar dist0 = (p0-pointWallDist[e[0]].origin()) & eVec;
586  scalar dist1 = (pointWallDist[e[1]].origin()-p1) & eVec;
587  scalar s = 0.5*(dist1+eMag+dist0);
588 
589  point medialAxisPt;
590  if (s <= dist0)
591  {
592  medialAxisPt = p0;
593  }
594  else if (s >= dist0+eMag)
595  {
596  medialAxisPt = p1;
597  }
598  else
599  {
600  medialAxisPt = p0+(s-dist0)*eVec;
601  }
602 
603  forAll(e, ep)
604  {
605  label pointi = e[ep];
606 
607  if (!pointMedialDist[pointi].valid(dummyTrackData))
608  {
609  maxPoints.append(pointi);
610  maxInfo.append
611  (
612  pointData
613  (
614  medialAxisPt, //points[pointi],
615  magSqr(points[pointi]-medialAxisPt),//0.0,
616  pointi, // passive data
617  Zero // passive data
618  )
619  );
620  pointMedialDist[pointi] = maxInfo.last();
621  }
622  }
623  }
624  }
625  }
626 
627 
628  // 2. Seed non-adapt patches
629  const polyBoundaryMesh& patches = mesh().boundaryMesh();
630 
631  labelHashSet adaptPatches(adaptPatchIDs_);
632 
633 
634  forAll(patches, patchi)
635  {
636  const polyPatch& pp = patches[patchi];
637  const pointPatchVectorField& pvf =
638  pointDisplacement().boundaryField()[patchi];
639 
640  if
641  (
642  !pp.coupled()
643  && !isA<emptyPolyPatch>(pp)
644  && !adaptPatches.found(patchi)
645  )
646  {
647  const labelList& meshPoints = pp.meshPoints();
648 
649  // Check the type of the patchField. The types are
650  // - fixedValue (0 or more layers) but the >0 layers have
651  // already been handled in the adaptPatches loop
652  // - constraint (but not coupled) types, e.g. symmetryPlane,
653  // slip.
654  if (pvf.fixesValue())
655  {
656  // Disable all movement on fixedValue patchFields
657  Info<< typeName
658  << " : Inserting all points on patch " << pp.name()
659  << endl;
660 
661  forAll(meshPoints, i)
662  {
663  label pointi = meshPoints[i];
664  if (!pointMedialDist[pointi].valid(dummyTrackData))
665  {
666  maxPoints.append(pointi);
667  maxInfo.append
668  (
669  pointData
670  (
671  points[pointi],
672  0.0,
673  pointi, // passive data
674  Zero // passive data
675  )
676  );
677  pointMedialDist[pointi] = maxInfo.last();
678  }
679  }
680  }
681  else
682  {
683  // Based on geometry: analyse angle w.r.t. nearest moving
684  // point. In the pointWallDist we transported the
685  // normal as the passive vector. Note that this points
686  // out of the originating wall so inside of the domain
687  // on this patch.
688  Info<< typeName
689  << " : Inserting points on patch " << pp.name()
690  << " if angle to nearest layer patch > "
691  << slipFeatureAngle << " degrees." << endl;
692 
693  scalar slipFeatureAngleCos = Foam::cos
694  (
695  degToRad(slipFeatureAngle)
696  );
697  pointField pointNormals
698  (
700  );
701 
702  forAll(meshPoints, i)
703  {
704  label pointi = meshPoints[i];
705 
706  if
707  (
708  pointWallDist[pointi].valid(dummyTrackData)
709  && !pointMedialDist[pointi].valid(dummyTrackData)
710  )
711  {
712  // Check if angle not too large.
713  scalar cosAngle =
714  (
715  -pointWallDist[pointi].v()
716  & pointNormals[i]
717  );
718  if (cosAngle > slipFeatureAngleCos)
719  {
720  // Extrusion direction practically perpendicular
721  // to the patch. Disable movement at the patch.
722 
723  maxPoints.append(pointi);
724  maxInfo.append
725  (
726  pointData
727  (
728  points[pointi],
729  0.0,
730  pointi, // passive data
731  Zero // passive data
732  )
733  );
734  pointMedialDist[pointi] = maxInfo.last();
735  }
736  else
737  {
738  // Extrusion direction makes angle with patch
739  // so allow sliding - don't insert zero points
740  }
741  }
742  }
743  }
744  }
745  }
746 
747  maxInfo.shrink();
748  maxPoints.shrink();
749 
750  // Do all calculations
751  PointEdgeWave<pointData> medialDistCalc
752  (
753  mesh(),
754  maxPoints,
755  maxInfo,
756 
757  pointMedialDist,
758  edgeMedialDist,
759  0, // max iterations
760  dummyTrackData
761  );
762  medialDistCalc.iterate(2*nMedialAxisIter);
763 
764 
765  // Extract medial axis distance as pointScalarField
766  forAll(pointMedialDist, pointi)
767  {
768  if (pointMedialDist[pointi].valid(dummyTrackData))
769  {
770  medialDist_[pointi] = Foam::sqrt
771  (
772  pointMedialDist[pointi].distSqr()
773  );
774  medialVec_[pointi] = pointMedialDist[pointi].origin();
775  }
776  else
777  {
778  // Unvisited. Do as if on medial axis so unmoving
779  medialDist_[pointi] = 0.0;
780  medialVec_[pointi] = point(1, 0, 0);
781  }
782  }
783  }
784 
785  // Extract transported surface normals as pointVectorField
786  forAll(dispVec_, i)
787  {
788  if (!pointWallDist[i].valid(dummyTrackData))
789  {
790  dispVec_[i] = vector(1, 0, 0);
791  }
792  else
793  {
794  dispVec_[i] = pointWallDist[i].v();
795  }
796  }
797 
798  // Smooth normal vectors. Do not change normals on pp.meshPoints
799  smoothNormals
800  (
801  nSmoothNormals,
802  isMeshMasterPoint,
803  isMeshMasterEdge,
804  meshPoints,
805  dispVec_
806  );
807 
808  // Calculate ratio point medial distance to point wall distance
809  forAll(medialRatio_, pointi)
810  {
811  if (!pointWallDist[pointi].valid(dummyTrackData))
812  {
813  medialRatio_[pointi] = 0.0;
814  }
815  else
816  {
817  scalar wDist2 = pointWallDist[pointi].distSqr();
818  scalar mDist = medialDist_[pointi];
819 
820  if (wDist2 < sqr(SMALL) && mDist < SMALL)
821  //- Note: maybe less strict:
822  //(
823  // wDist2 < sqr(meshRefiner_.mergeDistance())
824  // && mDist < meshRefiner_.mergeDistance()
825  //)
826  {
827  medialRatio_[pointi] = 0.0;
828  }
829  else
830  {
831  medialRatio_[pointi] = mDist / (Foam::sqrt(wDist2) + mDist);
832  }
833  }
834  }
835 
836 
837  if (debug)
838  {
839  Info<< typeName
840  << " : Writing medial axis fields:" << nl
841  << incrIndent
842  << "ratio of medial distance to wall distance : "
843  << medialRatio_.name() << nl
844  << "distance to nearest medial axis : "
845  << medialDist_.name() << nl
846  << "nearest medial axis location : "
847  << medialVec_.name() << nl
848  << "normal at nearest wall : "
849  << dispVec_.name() << nl
850  << decrIndent << nl
851  << endl;
852 
853  dispVec_.write();
854  medialRatio_.write();
855  medialDist_.write();
856  medialVec_.write();
857  }
858 }
859 
860 
861 bool Foam::medialAxisMeshMover::unmarkExtrusion
862 (
863  const label patchPointi,
864  pointField& patchDisp,
865  List<snappyLayerDriver::extrudeMode>& extrudeStatus
866 )
867 {
868  if (extrudeStatus[patchPointi] == snappyLayerDriver::EXTRUDE)
869  {
870  extrudeStatus[patchPointi] = snappyLayerDriver::NOEXTRUDE;
871  patchDisp[patchPointi] = Zero;
872  return true;
873  }
874  else if (extrudeStatus[patchPointi] == snappyLayerDriver::EXTRUDEREMOVE)
875  {
876  extrudeStatus[patchPointi] = snappyLayerDriver::NOEXTRUDE;
877  patchDisp[patchPointi] = Zero;
878  return true;
879  }
880  else
881  {
882  return false;
883  }
884 }
885 
886 
887 void Foam::medialAxisMeshMover::syncPatchDisplacement
888 (
889  const scalarField& minThickness,
890  pointField& patchDisp,
891  List<snappyLayerDriver::extrudeMode>& extrudeStatus
892 ) const
893 {
894  const indirectPrimitivePatch& pp = adaptPatchPtr_();
895  const labelList& meshPoints = pp.meshPoints();
896 
897  label nChangedTotal = 0;
898 
899  while (true)
900  {
901  label nChanged = 0;
902 
903  // Sync displacement (by taking min)
905  (
906  mesh(),
907  meshPoints,
908  patchDisp,
909  minMagSqrEqOp<vector>(),
910  point::rootMax // null value
911  );
912 
913  // Unmark if displacement too small
914  forAll(patchDisp, i)
915  {
916  if (mag(patchDisp[i]) < minThickness[i])
917  {
918  if (unmarkExtrusion(i, patchDisp, extrudeStatus))
919  {
920  nChanged++;
921  }
922  }
923  }
924 
925  //labelList syncPatchNLayers(patchNLayers);
926  //
927  //syncTools::syncPointList
928  //(
929  // mesh(),
930  // meshPoints,
931  // syncPatchNLayers,
932  // minEqOp<label>(),
933  // labelMax // null value
934  //);
935  //
938  //forAll(syncPatchNLayers, i)
939  //{
940  // if (syncPatchNLayers[i] != patchNLayers[i])
941  // {
942  // if
943  // (
944  // unmarkExtrusion
945  // (
946  // i,
947  // patchDisp,
948  // patchNLayers,
949  // extrudeStatus
950  // )
951  // )
952  // {
953  // nChanged++;
954  // }
955  // }
956  //}
957  //
958  //syncTools::syncPointList
959  //(
960  // mesh(),
961  // meshPoints,
962  // syncPatchNLayers,
963  // maxEqOp<label>(),
964  // labelMin // null value
965  //);
966  //
969  //forAll(syncPatchNLayers, i)
970  //{
971  // if (syncPatchNLayers[i] != patchNLayers[i])
972  // {
973  // if
974  // (
975  // unmarkExtrusion
976  // (
977  // i,
978  // patchDisp,
979  // patchNLayers,
980  // extrudeStatus
981  // )
982  // )
983  // {
984  // nChanged++;
985  // }
986  // }
987  //}
988 
989  nChangedTotal += nChanged;
990 
991  if (!returnReduce(nChanged, sumOp<label>()))
992  {
993  break;
994  }
995  }
996 
997  //Info<< "Prevented extrusion on "
998  // << returnReduce(nChangedTotal, sumOp<label>())
999  // << " coupled patch points during syncPatchDisplacement." << endl;
1000 }
1001 
1002 
1003 void Foam::medialAxisMeshMover::minSmoothField
1004 (
1005  const label nSmoothDisp,
1006  const PackedBoolList& isPatchMasterPoint,
1007  const PackedBoolList& isPatchMasterEdge,
1008  const scalarField& fieldMin,
1009  scalarField& field
1010 ) const
1011 {
1012  const indirectPrimitivePatch& pp = adaptPatchPtr_();
1013  const edgeList& edges = pp.edges();
1014  const labelList& meshPoints = pp.meshPoints();
1015 
1016  scalarField edgeWeights(edges.size());
1017  scalarField invSumWeight(meshPoints.size());
1019  (
1020  mesh(),
1021  isPatchMasterEdge,
1022  meshPoints,
1023  edges,
1024  edgeWeights,
1025  invSumWeight
1026  );
1027 
1028  // Get smoothly varying patch field.
1029  Info<< typeName << " : Smoothing field ..." << endl;
1030 
1031  for (label iter = 0; iter < nSmoothDisp; iter++)
1032  {
1033  scalarField average(pp.nPoints());
1035  (
1036  mesh(),
1037  isPatchMasterEdge,
1038  meshPoints,
1039  edges,
1040  edgeWeights,
1041  field,
1042  average
1043  );
1044  average *= invSumWeight;
1045 
1046  // Transfer to field
1047  forAll(field, pointi)
1048  {
1049  //full smoothing neighbours + point value
1050  average[pointi] = 0.5*(field[pointi]+average[pointi]);
1051 
1052  // perform monotonic smoothing
1053  if
1054  (
1055  average[pointi] < field[pointi]
1056  && average[pointi] >= fieldMin[pointi]
1057  )
1058  {
1059  field[pointi] = average[pointi];
1060  }
1061  }
1062 
1063  // Do residual calculation every so often.
1064  if ((iter % 10) == 0)
1065  {
1066  scalar resid = meshRefinement::gAverage
1067  (
1068  isPatchMasterPoint,
1069  mag(field-average)()
1070  );
1071  Info<< " Iteration " << iter << " residual " << resid << endl;
1072  }
1073  }
1074 }
1075 
1076 // Stop layer growth where mesh wraps around edge with a
1077 // large feature angle
1078 void Foam::medialAxisMeshMover::
1079 handleFeatureAngleLayerTerminations
1080 (
1081  const scalar minCos,
1082  const PackedBoolList& isPatchMasterPoint,
1083  const labelList& meshEdges,
1084  List<snappyLayerDriver::extrudeMode>& extrudeStatus,
1085  pointField& patchDisp,
1086  label& nPointCounter
1087 ) const
1088 {
1089  const indirectPrimitivePatch& pp = adaptPatchPtr_();
1090 
1091  // Mark faces that have all points extruded
1092  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1093 
1094  boolList extrudedFaces(pp.size(), true);
1095 
1096  forAll(pp.localFaces(), facei)
1097  {
1098  const face& f = pp.localFaces()[facei];
1099 
1100  forAll(f, fp)
1101  {
1102  if (extrudeStatus[f[fp]] == snappyLayerDriver::NOEXTRUDE)
1103  {
1104  extrudedFaces[facei] = false;
1105  break;
1106  }
1107  }
1108  }
1109 
1110 
1111 
1112  //label nOldPointCounter = nPointCounter;
1113 
1114  // Detect situation where two featureedge-neighbouring faces are partly or
1115  // not extruded and the edge itself is extruded. In this case unmark the
1116  // edge for extrusion.
1117 
1118 
1119  List<List<point>> edgeFaceNormals(pp.nEdges());
1120  List<List<bool>> edgeFaceExtrude(pp.nEdges());
1121 
1122  const labelListList& edgeFaces = pp.edgeFaces();
1123  const vectorField& faceNormals = pp.faceNormals();
1124 
1125  forAll(edgeFaces, edgeI)
1126  {
1127  const labelList& eFaces = edgeFaces[edgeI];
1128 
1129  edgeFaceNormals[edgeI].setSize(eFaces.size());
1130  edgeFaceExtrude[edgeI].setSize(eFaces.size());
1131  forAll(eFaces, i)
1132  {
1133  label facei = eFaces[i];
1134  edgeFaceNormals[edgeI][i] = faceNormals[facei];
1135  edgeFaceExtrude[edgeI][i] = extrudedFaces[facei];
1136  }
1137  }
1138 
1140  (
1141  mesh(),
1142  meshEdges,
1143  edgeFaceNormals,
1144  globalMeshData::ListPlusEqOp<List<point>>(), // combine operator
1145  List<point>() // null value
1146  );
1147 
1149  (
1150  mesh(),
1151  meshEdges,
1152  edgeFaceExtrude,
1153  globalMeshData::ListPlusEqOp<List<bool>>(), // combine operator
1154  List<bool>() // null value
1155  );
1156 
1157 
1158  forAll(edgeFaceNormals, edgeI)
1159  {
1160  const List<point>& eFaceNormals = edgeFaceNormals[edgeI];
1161  const List<bool>& eFaceExtrude = edgeFaceExtrude[edgeI];
1162 
1163  if (eFaceNormals.size() == 2)
1164  {
1165  const edge& e = pp.edges()[edgeI];
1166  label v0 = e[0];
1167  label v1 = e[1];
1168 
1169  if
1170  (
1171  extrudeStatus[v0] != snappyLayerDriver::NOEXTRUDE
1172  || extrudeStatus[v1] != snappyLayerDriver::NOEXTRUDE
1173  )
1174  {
1175  if (!eFaceExtrude[0] || !eFaceExtrude[1])
1176  {
1177  const vector& n0 = eFaceNormals[0];
1178  const vector& n1 = eFaceNormals[1];
1179 
1180  if ((n0 & n1) < minCos)
1181  {
1182  if (unmarkExtrusion(v0, patchDisp, extrudeStatus))
1183  {
1184  if (isPatchMasterPoint[v0])
1185  {
1186  nPointCounter++;
1187  }
1188  }
1189  if (unmarkExtrusion(v1, patchDisp, extrudeStatus))
1190  {
1191  if (isPatchMasterPoint[v1])
1192  {
1193  nPointCounter++;
1194  }
1195  }
1196  }
1197  }
1198  }
1199  }
1200  }
1201 
1202  //Info<< "Added "
1203  // << returnReduce(nPointCounter-nOldPointCounter, sumOp<label>())
1204  // << " point not to extrude." << endl;
1205 }
1206 
1207 
1208 // Find isolated islands (points, edges and faces and layer terminations)
1209 // in the layer mesh and stop any layer growth at these points.
1210 void Foam::medialAxisMeshMover::findIsolatedRegions
1211 (
1212  const scalar minCosLayerTermination,
1213  const bool detectExtrusionIsland,
1214  const PackedBoolList& isPatchMasterPoint,
1215  const PackedBoolList& isPatchMasterEdge,
1216  const labelList& meshEdges,
1217  const scalarField& minThickness,
1218  List<snappyLayerDriver::extrudeMode>& extrudeStatus,
1219  pointField& patchDisp
1220 ) const
1221 {
1222  const indirectPrimitivePatch& pp = adaptPatchPtr_();
1223  const labelListList& pointFaces = pp.pointFaces();
1224  const labelList& meshPoints = pp.meshPoints();
1225 
1226  Info<< typeName << " : Removing isolated regions ..." << endl;
1227 
1228  // Keep count of number of points unextruded
1229  label nPointCounter = 0;
1230 
1231 
1232  autoPtr<OBJstream> str;
1233  if (debug)
1234  {
1235  str.reset
1236  (
1237  new OBJstream
1238  (
1239  mesh().time().path()
1240  / "islandExcludePoints_"
1241  + mesh().time().timeName()
1242  + ".obj"
1243  )
1244  );
1245  Info<< typeName
1246  << " : Writing points surrounded by non-extruded points to "
1247  << str().name() << endl;
1248  }
1249 
1250  while (true)
1251  {
1252  // Stop layer growth where mesh wraps around edge with a
1253  // large feature angle
1254  handleFeatureAngleLayerTerminations
1255  (
1256  minCosLayerTermination,
1257  isPatchMasterPoint,
1258  meshEdges,
1259 
1260  extrudeStatus,
1261  patchDisp,
1262  nPointCounter
1263  );
1264 
1265  syncPatchDisplacement(minThickness, patchDisp, extrudeStatus);
1266 
1267 
1268 
1269  // Detect either:
1270  // - point where all surrounding points are not extruded
1271  // (detectExtrusionIsland)
1272  // or
1273  // - point where all the faces surrounding it are not fully
1274  // extruded
1275 
1276  boolList keptPoints(pp.nPoints(), false);
1277 
1278  if (detectExtrusionIsland)
1279  {
1280  // Do not extrude from point where all neighbouring
1281  // points are not grown
1282  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1283 
1284  labelList islandPoint(pp.size(), -1);
1285  forAll(pp, facei)
1286  {
1287  const face& f = pp.localFaces()[facei];
1288 
1289  forAll(f, fp)
1290  {
1291  if (extrudeStatus[f[fp]] != snappyLayerDriver::NOEXTRUDE)
1292  {
1293  if (islandPoint[facei] == -1)
1294  {
1295  // First point to extrude
1296  islandPoint[facei] = f[fp];
1297  }
1298  else if (islandPoint[facei] != -2)
1299  {
1300  // Second or more point to extrude
1301  islandPoint[facei] = -2;
1302  }
1303  }
1304  }
1305  }
1306 
1307  // islandPoint:
1308  // -1 : no point extruded on face
1309  // -2 : >= 2 points extruded on face
1310  // >=0: label of point extruded
1311 
1312  // Check all surrounding faces that I am the islandPoint
1313  forAll(pointFaces, patchPointi)
1314  {
1315  if (extrudeStatus[patchPointi] != snappyLayerDriver::NOEXTRUDE)
1316  {
1317  const labelList& pFaces = pointFaces[patchPointi];
1318 
1319  forAll(pFaces, i)
1320  {
1321  label facei = pFaces[i];
1322  if (islandPoint[facei] != patchPointi)
1323  {
1324  keptPoints[patchPointi] = true;
1325  break;
1326  }
1327  }
1328  }
1329  }
1330  }
1331  else
1332  {
1333  // Do not extrude from point where all neighbouring
1334  // faces are not grown
1335  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1336 
1337  boolList extrudedFaces(pp.size(), true);
1338  forAll(pp.localFaces(), facei)
1339  {
1340  const face& f = pp.localFaces()[facei];
1341  forAll(f, fp)
1342  {
1343  if (extrudeStatus[f[fp]] == snappyLayerDriver::NOEXTRUDE)
1344  {
1345  extrudedFaces[facei] = false;
1346  break;
1347  }
1348  }
1349  }
1350 
1351  const labelListList& pointFaces = pp.pointFaces();
1352 
1353  forAll(keptPoints, patchPointi)
1354  {
1355  const labelList& pFaces = pointFaces[patchPointi];
1356 
1357  forAll(pFaces, i)
1358  {
1359  label facei = pFaces[i];
1360  if (extrudedFaces[facei])
1361  {
1362  keptPoints[patchPointi] = true;
1363  break;
1364  }
1365  }
1366  }
1367  }
1368 
1369 
1371  (
1372  mesh(),
1373  meshPoints,
1374  keptPoints,
1375  orEqOp<bool>(),
1376  false // null value
1377  );
1378 
1379  label nChanged = 0;
1380 
1381  forAll(keptPoints, patchPointi)
1382  {
1383  if (!keptPoints[patchPointi])
1384  {
1385  if (unmarkExtrusion(patchPointi, patchDisp, extrudeStatus))
1386  {
1387  nPointCounter++;
1388  nChanged++;
1389 
1390  if (str.valid())
1391  {
1392  str().write(pp.points()[meshPoints[patchPointi]]);
1393  }
1394  }
1395  }
1396  }
1397 
1398 
1399  if (returnReduce(nChanged, sumOp<label>()) == 0)
1400  {
1401  break;
1402  }
1403  }
1404 
1405  const edgeList& edges = pp.edges();
1406 
1407 
1408  // Count number of mesh edges using a point
1409  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1410 
1411  labelList isolatedPoint(pp.nPoints(),0);
1412 
1413  forAll(edges, edgeI)
1414  {
1415  if (isPatchMasterEdge[edgeI])
1416  {
1417  const edge& e = edges[edgeI];
1418 
1419  label v0 = e[0];
1420  label v1 = e[1];
1421 
1422  if (extrudeStatus[v1] != snappyLayerDriver::NOEXTRUDE)
1423  {
1424  isolatedPoint[v0] += 1;
1425  }
1426  if (extrudeStatus[v0] != snappyLayerDriver::NOEXTRUDE)
1427  {
1428  isolatedPoint[v1] += 1;
1429  }
1430  }
1431  }
1432 
1434  (
1435  mesh(),
1436  meshPoints,
1437  isolatedPoint,
1438  plusEqOp<label>(),
1439  label(0) // null value
1440  );
1441 
1442  // stop layer growth on isolated faces
1443  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1444  forAll(pp, facei)
1445  {
1446  const face& f = pp.localFaces()[facei];
1447  bool failed = false;
1448  forAll(f, fp)
1449  {
1450  if (isolatedPoint[f[fp]] > 2)
1451  {
1452  failed = true;
1453  break;
1454  }
1455  }
1456  bool allPointsExtruded = true;
1457  if (!failed)
1458  {
1459  forAll(f, fp)
1460  {
1461  if (extrudeStatus[f[fp]] == snappyLayerDriver::NOEXTRUDE)
1462  {
1463  allPointsExtruded = false;
1464  break;
1465  }
1466  }
1467 
1468  if (allPointsExtruded)
1469  {
1470  forAll(f, fp)
1471  {
1472  if
1473  (
1474  unmarkExtrusion
1475  (
1476  f[fp],
1477  patchDisp,
1478  extrudeStatus
1479  )
1480  )
1481  {
1482  nPointCounter++;
1483 
1484  if (str.valid())
1485  {
1486  str().write(pp.points()[meshPoints[f[fp]]]);
1487  }
1488  }
1489  }
1490  }
1491  }
1492  }
1493 
1494  reduce(nPointCounter, sumOp<label>());
1495  Info<< typeName
1496  << " : Number of isolated points extrusion stopped : "<< nPointCounter
1497  << endl;
1498 }
1499 
1500 
1501 void Foam::medialAxisMeshMover::smoothLambdaMuDisplacement
1502 (
1503  const label nSmoothDisp,
1504  const PackedBoolList& isMeshMasterPoint,
1505  const PackedBoolList& isMeshMasterEdge,
1506  vectorField& displacement
1507 ) const
1508 {
1509  const edgeList& edges = mesh().edges();
1510 
1511  // Correspondence between local edges/points and mesh edges/points
1512  const labelList meshPoints(identity(mesh().nPoints()));
1513 
1514  // Calculate inverse sum of weights
1515  scalarField edgeWeights(mesh().nEdges());
1516  scalarField invSumWeight(meshPoints.size());
1518  (
1519  mesh(),
1520  isMeshMasterEdge,
1521  meshPoints,
1522  edges,
1523  edgeWeights,
1524  invSumWeight
1525  );
1526 
1527  // Get smoothly varying patch field.
1528  Info<< typeName << " : Smoothing displacement ..." << endl;
1529 
1530  const scalar lambda = 0.33;
1531  const scalar mu = -0.34;
1532 
1534 
1535  for (label iter = 0; iter < nSmoothDisp; iter++)
1536  {
1538  (
1539  mesh(),
1540  isMeshMasterEdge,
1541  meshPoints,
1542  edges,
1543  edgeWeights,
1544  displacement,
1545  average
1546  );
1547  average *= invSumWeight;
1548 
1549  forAll(displacement, i)
1550  {
1551  if (medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL)
1552  {
1553  displacement[i] = (1-lambda)*displacement[i]+lambda*average[i];
1554  }
1555  }
1556 
1558  (
1559  mesh(),
1560  isMeshMasterEdge,
1561  meshPoints,
1562  edges,
1563  edgeWeights,
1564  displacement,
1565  average
1566  );
1567  average *= invSumWeight;
1568 
1569 
1570  forAll(displacement, i)
1571  {
1572  if (medialRatio_[i] > SMALL && medialRatio_[i] < 1-SMALL)
1573  {
1574  displacement[i] = (1-mu)*displacement[i]+mu*average[i];
1575  }
1576  }
1577 
1578 
1579  // Do residual calculation every so often.
1580  if ((iter % 10) == 0)
1581  {
1582  scalar resid = meshRefinement::gAverage
1583  (
1584  isMeshMasterPoint,
1585  mag(displacement-average)()
1586  );
1587  Info<< " Iteration " << iter << " residual " << resid << endl;
1588  }
1589  }
1590 }
1591 
1592 
1593 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1594 
1595 Foam::medialAxisMeshMover::medialAxisMeshMover
1597  const dictionary& dict,
1598  const List<labelPair>& baffles,
1599  pointVectorField& pointDisplacement
1600 )
1601 :
1602  externalDisplacementMeshMover(dict, baffles, pointDisplacement),
1603  adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
1604  adaptPatchPtr_(getPatch(mesh(), adaptPatchIDs_)),
1605  scale_
1606  (
1607  IOobject
1608  (
1609  "scale",
1610  pointDisplacement.time().timeName(),
1611  pointDisplacement.db(),
1614  ),
1615  pMesh(),
1616  dimensionedScalar("scale", dimless, 1.0)
1617  ),
1618  oldPoints_(mesh().points()),
1619  meshMover_
1620  (
1621  const_cast<polyMesh&>(mesh()),
1622  const_cast<pointMesh&>(pMesh()),
1623  adaptPatchPtr_(),
1624  pointDisplacement,
1625  scale_,
1626  oldPoints_,
1627  adaptPatchIDs_,
1628  dict
1629  ),
1630  dispVec_
1631  (
1632  IOobject
1633  (
1634  "dispVec",
1635  pointDisplacement.time().timeName(),
1636  pointDisplacement.db(),
1639  false
1640  ),
1641  pMesh(),
1642  dimensionedVector("dispVec", dimLength, Zero)
1643  ),
1644  medialRatio_
1645  (
1646  IOobject
1647  (
1648  "medialRatio",
1649  pointDisplacement.time().timeName(),
1650  pointDisplacement.db(),
1653  false
1654  ),
1655  pMesh(),
1656  dimensionedScalar("medialRatio", dimless, 0.0)
1657  ),
1658  medialDist_
1659  (
1660  IOobject
1661  (
1662  "pointMedialDist",
1663  pointDisplacement.time().timeName(),
1664  pointDisplacement.db(),
1667  false
1668  ),
1669  pMesh(),
1670  dimensionedScalar("pointMedialDist", dimLength, 0.0)
1671  ),
1672  medialVec_
1673  (
1674  IOobject
1675  (
1676  "medialVec",
1677  pointDisplacement.time().timeName(),
1678  pointDisplacement.db(),
1681  false
1682  ),
1683  pMesh(),
1684  dimensionedVector("medialVec", dimLength, Zero)
1685  )
1686 {
1687  update(dict);
1688 }
1689 
1690 
1691 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1692 
1694 {}
1695 
1696 
1697 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1698 
1699 void Foam::medialAxisMeshMover::calculateDisplacement
1700 (
1701  const dictionary& coeffDict,
1702  const scalarField& minThickness,
1703  List<snappyLayerDriver::extrudeMode>& extrudeStatus,
1704  pointField& patchDisp
1705 )
1706 {
1707  Info<< typeName << " : Smoothing using Medial Axis ..." << endl;
1708 
1709  const indirectPrimitivePatch& pp = adaptPatchPtr_;
1710  const labelList& meshPoints = pp.meshPoints();
1711 
1712 
1713  // Read settings
1714  // ~~~~~~~~~~~~~
1715 
1716  //- (lambda-mu) smoothing of internal displacement
1717  const label nSmoothDisplacement = coeffDict.lookupOrDefault
1718  (
1719  "nSmoothDisplacement",
1720  0
1721  );
1722 
1723  //- Layer thickness too big
1724  const scalar maxThicknessToMedialRatio = readScalar
1725  (
1726  coeffDict.lookup("maxThicknessToMedialRatio")
1727  );
1728 
1729  //- Feature angle when to stop adding layers
1730  const scalar featureAngle = readScalar(coeffDict.lookup("featureAngle"));
1731 
1732  //- Stop layer growth where mesh wraps around sharp edge
1733  const scalar minCosLayerTermination = Foam::cos
1734  (
1735  degToRad(0.5*featureAngle)
1736  );
1737 
1738  //- Smoothing wanted patch thickness
1739  const label nSmoothPatchThickness = readLabel
1740  (
1741  coeffDict.lookup("nSmoothThickness")
1742  );
1743 
1744  //- Number of edges walking out
1745  const label nMedialAxisIter = coeffDict.lookupOrDefault<label>
1746  (
1747  "nMedialAxisIter",
1749  );
1750 
1751  //- Use strict extrusionIsland detection
1752  const Switch detectExtrusionIsland = coeffDict.lookupOrDefault<Switch>
1753  (
1754  "detectExtrusionIsland",
1755  false
1756  );
1757 
1758 
1759  // Precalulate master points/edge (only relevant for shared points/edges)
1760  const PackedBoolList isMeshMasterPoint(syncTools::getMasterPoints(mesh()));
1761  const PackedBoolList isMeshMasterEdge(syncTools::getMasterEdges(mesh()));
1762  // Precalculate meshEdge per pp edge
1763  const labelList meshEdges
1764  (
1765  pp.meshEdges
1766  (
1767  mesh().edges(),
1768  mesh().pointEdges()
1769  )
1770  );
1771 
1772  // Precalulate (patch) master point/edge
1773  const PackedBoolList isPatchMasterPoint
1774  (
1776  (
1777  mesh(),
1778  meshPoints
1779  )
1780  );
1781  const PackedBoolList isPatchMasterEdge
1782  (
1784  (
1785  mesh(),
1786  meshEdges
1787  )
1788  );
1789 
1790 
1791  scalarField thickness(patchDisp.size());
1792 
1793  thickness = mag(patchDisp);
1794 
1795  forAll(thickness, patchPointi)
1796  {
1797  if (extrudeStatus[patchPointi] == snappyLayerDriver::NOEXTRUDE)
1798  {
1799  thickness[patchPointi] = 0.0;
1800  }
1801  }
1802 
1803  label numThicknessRatioExclude = 0;
1804 
1805  // reduce thickness where thickness/medial axis distance large
1806  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1807 
1808  autoPtr<OBJstream> str;
1809  if (debug)
1810  {
1811  str.reset
1812  (
1813  new OBJstream
1814  (
1815  mesh().time().path()
1816  / "thicknessRatioExcludePoints_"
1817  + mesh().time().timeName()
1818  + ".obj"
1819  )
1820  );
1821  Info<< typeName
1822  << " : Writing points with too large an extrusion distance to "
1823  << str().name() << endl;
1824  }
1825 
1826  autoPtr<OBJstream> medialVecStr;
1827  if (debug)
1828  {
1829  medialVecStr.reset
1830  (
1831  new OBJstream
1832  (
1833  mesh().time().path()
1834  / "thicknessRatioExcludeMedialVec_"
1835  + mesh().time().timeName()
1836  + ".obj"
1837  )
1838  );
1839  Info<< typeName
1840  << " : Writing medial axis vectors on points with too large"
1841  << " an extrusion distance to " << medialVecStr().name() << endl;
1842  }
1843 
1844  forAll(meshPoints, patchPointi)
1845  {
1846  if (extrudeStatus[patchPointi] != snappyLayerDriver::NOEXTRUDE)
1847  {
1848  label pointi = meshPoints[patchPointi];
1849 
1850  //- Option 1: look only at extrusion thickness v.s. distance
1851  // to nearest (medial axis or static) point.
1852  scalar mDist = medialDist_[pointi];
1853  scalar thicknessRatio = thickness[patchPointi]/(mDist+VSMALL);
1854 
1855  //- Option 2: Look at component in the direction
1856  // of nearest (medial axis or static) point.
1857  vector n =
1858  patchDisp[patchPointi]
1859  / (mag(patchDisp[patchPointi]) + VSMALL);
1860  vector mVec = mesh().points()[pointi]-medialVec_[pointi];
1861  mVec /= mag(mVec)+VSMALL;
1862  thicknessRatio *= (n&mVec);
1863 
1864  if (thicknessRatio > maxThicknessToMedialRatio)
1865  {
1866  // Truncate thickness.
1867  if (debug&2)
1868  {
1869  Pout<< "truncating displacement at "
1870  << mesh().points()[pointi]
1871  << " from " << thickness[patchPointi]
1872  << " to "
1873  << 0.5
1874  *(
1875  minThickness[patchPointi]
1876  +thickness[patchPointi]
1877  )
1878  << " medial direction:" << mVec
1879  << " extrusion direction:" << n
1880  << " with thicknessRatio:" << thicknessRatio
1881  << endl;
1882  }
1883 
1884  thickness[patchPointi] =
1885  0.5*(minThickness[patchPointi]+thickness[patchPointi]);
1886 
1887  patchDisp[patchPointi] = thickness[patchPointi]*n;
1888 
1889  if (isPatchMasterPoint[patchPointi])
1890  {
1891  numThicknessRatioExclude++;
1892  }
1893 
1894  if (str.valid())
1895  {
1896  const point& pt = mesh().points()[pointi];
1897  str().write(linePointRef(pt, pt+patchDisp[patchPointi]));
1898  }
1899  if (medialVecStr.valid())
1900  {
1901  const point& pt = mesh().points()[pointi];
1902  medialVecStr().write
1903  (
1904  linePointRef
1905  (
1906  pt,
1907  medialVec_[pointi]
1908  )
1909  );
1910  }
1911  }
1912  }
1913  }
1914 
1915  reduce(numThicknessRatioExclude, sumOp<label>());
1916  Info<< typeName << " : Reducing layer thickness at "
1917  << numThicknessRatioExclude
1918  << " nodes where thickness to medial axis distance is large " << endl;
1919 
1920 
1921  // find points where layer growth isolated to a lone point, edge or face
1922 
1923  findIsolatedRegions
1924  (
1925  minCosLayerTermination,
1926  detectExtrusionIsland,
1927 
1928  isPatchMasterPoint,
1929  isPatchMasterEdge,
1930  meshEdges,
1931  minThickness,
1932 
1933  extrudeStatus,
1934  patchDisp
1935  );
1936 
1937  // Update thickess for changed extrusion
1938  forAll(thickness, patchPointi)
1939  {
1940  if (extrudeStatus[patchPointi] == snappyLayerDriver::NOEXTRUDE)
1941  {
1942  thickness[patchPointi] = 0.0;
1943  }
1944  }
1945 
1946 
1947  // smooth layer thickness on moving patch
1948  minSmoothField
1949  (
1950  nSmoothPatchThickness,
1951  isPatchMasterPoint,
1952  isPatchMasterEdge,
1953  minThickness,
1954 
1955  thickness
1956  );
1957 
1958 
1959  // Dummy additional info for PointEdgeWave
1960  int dummyTrackData = 0;
1961 
1962  List<pointData> pointWallDist(mesh().nPoints());
1963 
1964  const pointField& points = mesh().points();
1965  // 1. Calculate distance to points where displacement is specified.
1966  // This wave is used to transport layer thickness
1967  {
1968  // Distance to wall and medial axis on edges.
1969  List<pointData> edgeWallDist(mesh().nEdges());
1970  labelList wallPoints(meshPoints.size());
1971 
1972  // Seed data.
1973  List<pointData> wallInfo(meshPoints.size());
1974 
1975  forAll(meshPoints, patchPointi)
1976  {
1977  label pointi = meshPoints[patchPointi];
1978  wallPoints[patchPointi] = pointi;
1979  wallInfo[patchPointi] = pointData
1980  (
1981  points[pointi],
1982  0.0,
1983  thickness[patchPointi], // transport layer thickness
1984  Zero // passive vector
1985  );
1986  }
1987 
1988  // Do all calculations
1989  PointEdgeWave<pointData> wallDistCalc
1990  (
1991  mesh(),
1992  wallPoints,
1993  wallInfo,
1994  pointWallDist,
1995  edgeWallDist,
1996  0, // max iterations
1997  dummyTrackData
1998  );
1999  wallDistCalc.iterate(nMedialAxisIter);
2000  }
2001 
2002 
2003  // Calculate scaled displacement vector
2004  pointField& displacement = pointDisplacement_;
2005 
2006  forAll(displacement, pointi)
2007  {
2008  if (!pointWallDist[pointi].valid(dummyTrackData))
2009  {
2010  displacement[pointi] = Zero;
2011  }
2012  else
2013  {
2014  // 1) displacement on nearest wall point, scaled by medialRatio
2015  // (wall distance / medial distance)
2016  // 2) pointWallDist[pointi].s() is layer thickness transported
2017  // from closest wall point.
2018  // 3) shrink in opposite direction of addedPoints
2019  displacement[pointi] =
2020  -medialRatio_[pointi]
2021  *pointWallDist[pointi].s()
2022  *dispVec_[pointi];
2023  }
2024  }
2025 
2026 
2027  // Smear displacement away from fixed values (medialRatio=0 or 1)
2028  if (nSmoothDisplacement > 0)
2029  {
2030  smoothLambdaMuDisplacement
2031  (
2032  nSmoothDisplacement,
2033  isMeshMasterPoint,
2034  isMeshMasterEdge,
2035  displacement
2036  );
2037  }
2038 }
2039 
2040 
2041 bool Foam::medialAxisMeshMover::shrinkMesh
2042 (
2043  const dictionary& meshQualityDict,
2044  const label nAllowableErrors,
2045  labelList& checkFaces
2046 )
2047 {
2048  //- Number of attempts shrinking the mesh
2049  const label nSnap = readLabel(meshQualityDict.lookup("nRelaxIter"));
2050 
2051 
2052 
2053 
2054  // Make sure displacement boundary conditions is uptodate with
2055  // internal field
2056  meshMover_.setDisplacementPatchFields();
2057 
2058  Info<< typeName << " : Moving mesh ..." << endl;
2059  scalar oldErrorReduction = -1;
2060 
2061  bool meshOk = false;
2062 
2063  for (label iter = 0; iter < 2*nSnap ; iter++)
2064  {
2065  Info<< typeName
2066  << " : Iteration " << iter << endl;
2067  if (iter == nSnap)
2068  {
2069  Info<< typeName
2070  << " : Displacement scaling for error reduction set to 0."
2071  << endl;
2072  oldErrorReduction = meshMover_.setErrorReduction(0.0);
2073  }
2074 
2075  if
2076  (
2077  meshMover_.scaleMesh
2078  (
2079  checkFaces,
2080  baffles_,
2081  meshMover_.paramDict(),
2082  meshQualityDict,
2083  true,
2084  nAllowableErrors
2085  )
2086  )
2087  {
2088  Info<< typeName << " : Successfully moved mesh" << endl;
2089  meshOk = true;
2090  break;
2091  }
2092  }
2093 
2094  if (oldErrorReduction >= 0)
2095  {
2096  meshMover_.setErrorReduction(oldErrorReduction);
2097  }
2098 
2099  Info<< typeName << " : Finished moving mesh ..." << endl;
2100 
2101  return meshOk;
2102 }
2103 
2104 
2107  const dictionary& moveDict,
2108  const label nAllowableErrors,
2109  labelList& checkFaces
2110 )
2111 {
2112  // Read a few settings
2113  // ~~~~~~~~~~~~~~~~~~~
2114 
2115  //- Name of field specifying min thickness
2116  const word minThicknessName = word(moveDict.lookup("minThicknessName"));
2117 
2118 
2119  // The points have moved so before calculation update
2120  // the mesh and motionSolver accordingly
2121  movePoints(mesh().points());
2122  //
2124  //pointDisplacement_.boundaryField().updateCoeffs();
2125 
2126 
2127  // Extract out patch-wise displacement
2128  const indirectPrimitivePatch& pp = adaptPatchPtr_();
2129 
2130  scalarField zeroMinThickness;
2131  if (minThicknessName == "none")
2132  {
2133  zeroMinThickness = scalarField(pp.nPoints(), 0.0);
2134  }
2135  const scalarField& minThickness =
2136  (
2137  (minThicknessName == "none")
2138  ? zeroMinThickness
2139  : mesh().lookupObject<scalarField>(minThicknessName)
2140  );
2141 
2142 
2143  pointField patchDisp
2144  (
2145  pointDisplacement_.primitiveField(),
2146  pp.meshPoints()
2147  );
2148 
2150  (
2151  pp.nPoints(),
2153  );
2154  forAll(extrudeStatus, pointi)
2155  {
2156  if (mag(patchDisp[pointi]) <= minThickness[pointi]+SMALL)
2157  {
2158  extrudeStatus[pointi] = snappyLayerDriver::NOEXTRUDE;
2159  }
2160  }
2161 
2162 
2163  // Solve displacement
2164  calculateDisplacement(moveDict, minThickness, extrudeStatus, patchDisp);
2165 
2166  //- Move mesh according to calculated displacement
2167  return shrinkMesh
2168  (
2169  moveDict, // meshQualityDict,
2170  nAllowableErrors, // nAllowableErrors
2171  checkFaces
2172  );
2173 }
2174 
2175 
2177 {
2179 
2180  // Update local data for new geometry
2181  adaptPatchPtr_().movePoints(p);
2182 
2183  // Update motionSmoother for new geometry
2184  meshMover_.movePoints();
2185 
2186  // Assume corrent mesh location is correct
2187  meshMover_.correct();
2188 }
2189 
2190 
2191 // ************************************************************************* //
Variant of pointEdgePoint with some transported additional data. WIP - should be templated on data li...
Definition: pointData.H:61
bool isA(const Type &t)
Check if a dynamic_cast to typeid is possible.
Definition: typeInfo.H:134
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
label nPoints() const
Return number of points supporting patch faces.
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
virtual void movePoints(const pointField &)
Update local data for geometry changes.
A line primitive.
Definition: line.H:56
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:221
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
label iterate(const label maxIter)
Iterate until no changes or maxIter reached. Returns actual.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:715
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
Macros for easy insertion into run-time selection tables.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:48
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
static PackedBoolList getMasterEdges(const polyMesh &)
Get per edge whether it is uncoupled or a master of a.
Definition: syncTools.C:109
A list of faces which address into the list of points.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
void reset(T *=0)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
dynamicFvMesh & mesh
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
const pointField & points
List< edge > edgeList
Definition: edgeList.H:38
const Boundary & boundaryField() const
Return const-reference to the boundary field.
static void weightedSum(const polyMesh &mesh, const PackedBoolList &isMasterEdge, const labelList &meshPoints, const edgeList &edges, const scalarField &edgeWeights, const Field< Type > &data, Field< Type > &sum)
Helper: weighted sum (over all subset of mesh points) by.
static void calculateEdgeWeights(const polyMesh &mesh, const PackedBoolList &isMasterEdge, const labelList &meshPoints, const edgeList &edges, scalarField &edgeWeights, scalarField &invSumWeight)
Helper: calculate edge weights (1/length)
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
virtual bool move(const dictionary &, const label nAllowableErrors, labelList &checkFaces)
Move mesh using current pointDisplacement boundary values.
A class for handling words, derived from string.
Definition: word.H:59
label nPoints
Do not extrude. No layers added.
static PackedBoolList getMasterEdges(const polyMesh &mesh, const labelList &meshEdges)
Determine master edge for subset of edges. If coupled.
const Time & time() const
Return time.
Definition: IOobject.C:227
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual void movePoints(const pointField &)
Update local data for geometry changes.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){pointMap[i]=i;}for(label i=0;i< nPoints;i++){if(f[i] > 0.0){hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei){if(edges[ei].mag(points)< SMALL){label start=pointMap[edges[ei].start()];while(start!=pointMap[start]){start=pointMap[start];}label end=pointMap[edges[ei].end()];while(end!=pointMap[end]){end=pointMap[end];}label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;}}cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){cellShape &cs=cellShapes[celli];forAll(cs, i){cs[i]=pointMap[cs[i]];}cs.collapse();}label bcIDs[11]={-1, 0, 2, 4,-1, 5,-1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&symmetryPolyPatch::typeName,&wedgePolyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&symmetryPolyPatch::typeName,&oldCyclicPolyPatch::typeName};enum patchTypeNames{PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={"piston","valve","liner","cylinderHead","axis","wedge","inflow","outflow","presin","presout","symmetryPlane","cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
label nTotalPoints() const
Return total number of points in decomposed mesh. Not.
word timeName
Definition: getTimeIndex.H:3
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:91
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
label readLabel(Istream &is)
Definition: label.H:64
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set).
Definition: autoPtrI.H:83
static tmp< pointField > pointNormals(const polyMesh &, const PrimitivePatch< Face, FaceList, PointField, PointType > &)
Return parallel consistent point normals for patches using mesh points.
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1140
dimensioned< scalar > magSqr(const dimensioned< Type > &)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
virtual const fileName & name() const
Return the name of the stream.
Definition: IOstream.H:297
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
OFstream which keeps track of vertices.
Definition: OBJstream.H:53
static const char nl
Definition: Ostream.H:262
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
defineTypeNameAndDebug(combustionModel, 0)
static PackedBoolList getMasterPoints(const polyMesh &)
Get per point whether it is uncoupled or a master of a.
Definition: syncTools.C:65
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:237
pointPatchField< vector > pointPatchVectorField
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
const dimensionedScalar mu
Atomic mass unit.
Virtual base class for mesh movers with externally provided displacement field giving the boundary co...
void setSize(const label)
Reset size of List.
Definition: List.C:295
A bit-packed bool list.
static T gAverage(const PackedBoolList &isMasterElem, const UList< T > &values)
Helper: calculate average.
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedScalar lambda(laminarTransport.lookup("lambda"))
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
GeometricField< vector, pointPatchField, pointMesh > pointVectorField
Definition: pointFields.H:50
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Field< vector > vectorField
Specialisation of Field<T> for vector.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
virtual Ostream & write(const token &)=0
Write next token to stream.
volScalarField & p
static PackedBoolList getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
Wave propagation of information through grid. Every iteration information goes through one layer of e...
Definition: PointEdgeWave.H:85
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:230
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451