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