snappySnapDriver.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) 2011-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 Description
25  All to do with snapping to the surface
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "snappySnapDriver.H"
30 #include "motionSmoother.H"
31 #include "polyTopoChange.H"
32 #include "syncTools.H"
33 #include "fvMesh.H"
34 #include "Time.H"
35 #include "OFstream.H"
36 #include "OBJstream.H"
37 #include "mapPolyMesh.H"
38 #include "pointEdgePoint.H"
39 #include "PointEdgeWave.H"
40 #include "mergePoints.H"
41 #include "snapParameters.H"
42 #include "refinementSurfaces.H"
43 #include "unitConversion.H"
44 #include "localPointRegion.H"
45 #include "PatchTools.H"
46 
47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
48 
49 namespace Foam
50 {
51  defineTypeNameAndDebug(snappySnapDriver, 0);
52 }
53 
54 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
55 
56 Foam::label Foam::snappySnapDriver::getCollocatedPoints
57 (
58  const scalar tol,
59  const pointField& points,
60  PackedBoolList& isCollocatedPoint
61 )
62 {
63  labelList pointMap;
64  label nUnique = mergePoints
65  (
66  points, // points
67  tol, // mergeTol
68  false, // verbose
69  pointMap
70  );
71  bool hasMerged = (nUnique < points.size());
72 
73  if (!returnReduce(hasMerged, orOp<bool>()))
74  {
75  return 0;
76  }
77 
78  // Determine which merged points are referenced more than once
79  label nCollocated = 0;
80 
81  // Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
82  // twice)
83  labelList firstOldPoint(nUnique, -1);
84  forAll(pointMap, oldPointi)
85  {
86  label newPointi = pointMap[oldPointi];
87 
88  if (firstOldPoint[newPointi] == -1)
89  {
90  // First use of oldPointi. Store.
91  firstOldPoint[newPointi] = oldPointi;
92  }
93  else if (firstOldPoint[newPointi] == -2)
94  {
95  // Third or more reference of oldPointi -> non-manifold
96  isCollocatedPoint.set(oldPointi, 1u);
97  nCollocated++;
98  }
99  else
100  {
101  // Second reference of oldPointi -> non-manifold
102  isCollocatedPoint.set(firstOldPoint[newPointi], 1u);
103  nCollocated++;
104 
105  isCollocatedPoint.set(oldPointi, 1u);
106  nCollocated++;
107 
108  // Mark with special value to save checking next time round
109  firstOldPoint[newPointi] = -2;
110  }
111  }
112  return returnReduce(nCollocated, sumOp<label>());
113 }
114 
115 
116 Foam::pointField Foam::snappySnapDriver::smoothPatchDisplacement
117 (
118  const motionSmoother& meshMover,
119  const List<labelPair>& baffles
120 )
121 {
122  const indirectPrimitivePatch& pp = meshMover.patch();
123 
124  // Calculate geometrically non-manifold points on the patch to be moved.
125  PackedBoolList nonManifoldPoint(pp.nPoints());
126  label nNonManifoldPoints = getCollocatedPoints
127  (
128  small,
129  pp.localPoints(),
130  nonManifoldPoint
131  );
132  Info<< "Found " << nNonManifoldPoints << " non-manifold point(s)."
133  << endl;
134 
135 
136  // Average points
137  // ~~~~~~~~~~~~~~
138 
139  // We determine three points:
140  // - average of (centres of) connected patch faces
141  // - average of (centres of) connected internal mesh faces
142  // - as fallback: centre of any connected cell
143  // so we can do something moderately sensible for non/manifold points.
144 
145  // Note: the averages are calculated properly parallel. This is
146  // necessary to get the points shared by processors correct.
147 
148 
149  const labelListList& pointFaces = pp.pointFaces();
150  const labelList& meshPoints = pp.meshPoints();
151  const pointField& points = pp.points();
152  const polyMesh& mesh = meshMover.mesh();
153 
154  // Get labels of faces to count (master of coupled faces and baffle pairs)
155  PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
156 
157  {
158  forAll(baffles, i)
159  {
160  label f0 = baffles[i].first();
161  label f1 = baffles[i].second();
162 
163  if (isMasterFace.get(f0))
164  {
165  // Make f1 a slave
166  isMasterFace.unset(f1);
167  }
168  else if (isMasterFace.get(f1))
169  {
170  isMasterFace.unset(f0);
171  }
172  else
173  {
175  << "Both sides of baffle consisting of faces " << f0
176  << " and " << f1 << " are already slave faces."
177  << abort(FatalError);
178  }
179  }
180  }
181 
182 
183  // Get average position of boundary face centres
184  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
185 
186  vectorField avgBoundary(pointFaces.size(), Zero);
187  labelList nBoundary(pointFaces.size(), 0);
188 
189  forAll(pointFaces, patchPointi)
190  {
191  const labelList& pFaces = pointFaces[patchPointi];
192 
193  forAll(pFaces, pfI)
194  {
195  label facei = pFaces[pfI];
196 
197  if (isMasterFace.get(pp.addressing()[facei]))
198  {
199  avgBoundary[patchPointi] +=
200  pp[facei].centre(points) - pp.localPoints()[patchPointi];
201  nBoundary[patchPointi]++;
202  }
203  }
204  }
205 
207  (
208  mesh,
209  pp.meshPoints(),
210  avgBoundary,
211  plusEqOp<point>(), // combine op
212  vector::zero // null value
213  );
215  (
216  mesh,
217  pp.meshPoints(),
218  nBoundary,
219  plusEqOp<label>(), // combine op
220  label(0) // null value
221  );
222 
223  forAll(avgBoundary, i)
224  {
225  avgBoundary[i] = avgBoundary[i]/nBoundary[i] + pp.localPoints()[i];
226  }
227 
228 
229  // Get average position of internal face centres
230  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
231 
232  vectorField avgInternal;
233  labelList nInternal;
234  {
235  vectorField globalSum(mesh.nPoints(), Zero);
236  labelList globalNum(mesh.nPoints(), 0);
237 
238  // Note: no use of pointFaces
239  const faceList& faces = mesh.faces();
240 
241  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
242  {
243  const face& f = faces[facei];
244  const point& fc = mesh.faceCentres()[facei];
245 
246  forAll(f, fp)
247  {
248  globalSum[f[fp]] += fc - mesh.points()[f[fp]];
249  globalNum[f[fp]]++;
250  }
251  }
252 
253  // Count coupled faces as internal ones (but only once)
254  const polyBoundaryMesh& patches = mesh.boundaryMesh();
255 
256  forAll(patches, patchi)
257  {
258  if
259  (
260  patches[patchi].coupled()
261  && refCast<const coupledPolyPatch>(patches[patchi]).owner()
262  )
263  {
264  const coupledPolyPatch& pp =
265  refCast<const coupledPolyPatch>(patches[patchi]);
266 
267  const vectorField::subField faceCentres = pp.faceCentres();
268 
269  forAll(pp, i)
270  {
271  const face& f = pp[i];
272  const point& fc = faceCentres[i];
273 
274  forAll(f, fp)
275  {
276  globalSum[f[fp]] += fc - mesh.points()[f[fp]];
277  globalNum[f[fp]]++;
278  }
279  }
280  }
281  }
282 
284  (
285  mesh,
286  globalSum,
287  plusEqOp<vector>(), // combine op
288  vector::zero // null value
289  );
291  (
292  mesh,
293  globalNum,
294  plusEqOp<label>(), // combine op
295  label(0) // null value
296  );
297 
298  avgInternal.setSize(meshPoints.size());
299  nInternal.setSize(meshPoints.size());
300 
301  forAll(avgInternal, patchPointi)
302  {
303  label meshPointi = meshPoints[patchPointi];
304 
305  nInternal[patchPointi] = globalNum[meshPointi];
306 
307  if (nInternal[patchPointi] == 0)
308  {
309  avgInternal[patchPointi] = globalSum[meshPointi];
310  }
311  else
312  {
313  avgInternal[patchPointi] =
314  globalSum[meshPointi]/nInternal[patchPointi]
315  + mesh.points()[meshPointi];
316  }
317  }
318  }
319 
320 
321  // Precalculate any cell using mesh point (replacement of pointCells()[])
322  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
323 
324  labelList anyCell(mesh.nPoints(), -1);
325  forAll(mesh.faceNeighbour(), facei)
326  {
327  label own = mesh.faceOwner()[facei];
328  const face& f = mesh.faces()[facei];
329 
330  forAll(f, fp)
331  {
332  anyCell[f[fp]] = own;
333  }
334  }
335  for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
336  {
337  label own = mesh.faceOwner()[facei];
338 
339  const face& f = mesh.faces()[facei];
340 
341  forAll(f, fp)
342  {
343  anyCell[f[fp]] = own;
344  }
345  }
346 
347 
348  // Displacement to calculate.
349  pointField patchDisp(meshPoints.size(), Zero);
350 
351  forAll(pointFaces, i)
352  {
353  label meshPointi = meshPoints[i];
354  const point& currentPos = pp.points()[meshPointi];
355 
356  // Now we have the two average points: avgBoundary and avgInternal
357  // and how many boundary/internal faces connect to the point
358  // (nBoundary, nInternal)
359  // Do some blending between the two.
360  // Note: the following section has some reasoning behind it but the
361  // blending factors can be experimented with.
362 
363  point newPos;
364 
365  if (!nonManifoldPoint.get(i))
366  {
367  // Points that are manifold. Weight the internal and boundary
368  // by their number of faces and blend with
369  scalar internalBlend = 0.1;
370  scalar blend = 0.1;
371 
372  point avgPos =
373  (
374  internalBlend*nInternal[i]*avgInternal[i]
375  +(1-internalBlend)*nBoundary[i]*avgBoundary[i]
376  )
377  / (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
378 
379  newPos = (1-blend)*avgPos + blend*currentPos;
380  }
381  else if (nInternal[i] == 0)
382  {
383  // Non-manifold without internal faces. Use any connected cell
384  // as internal point instead. Use precalculated any cell to avoid
385  // e.g. pointCells()[meshPointi][0]
386 
387  const point& cc = mesh.cellCentres()[anyCell[meshPointi]];
388 
389  scalar cellCBlend = 0.8;
390  scalar blend = 0.1;
391 
392  point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
393 
394  newPos = (1-blend)*avgPos + blend*currentPos;
395  }
396  else
397  {
398  // Non-manifold point with internal faces connected to them
399  scalar internalBlend = 0.9;
400  scalar blend = 0.1;
401 
402  point avgPos =
403  internalBlend*avgInternal[i]
404  + (1-internalBlend)*avgBoundary[i];
405 
406  newPos = (1-blend)*avgPos + blend*currentPos;
407  }
408 
409  patchDisp[i] = newPos - currentPos;
410  }
411 
412  return patchDisp;
413 }
414 
415 
416 Foam::tmp<Foam::scalarField> Foam::snappySnapDriver::edgePatchDist
417 (
418  const pointMesh& pMesh,
419  const indirectPrimitivePatch& pp
420 )
421 {
422  const polyMesh& mesh = pMesh();
423 
424  // Set initial changed points to all the patch points
425  List<pointEdgePoint> wallInfo(pp.nPoints());
426 
427  forAll(pp.localPoints(), ppI)
428  {
429  wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
430  }
431 
432  // Current info on points
433  List<pointEdgePoint> allPointInfo(mesh.nPoints());
434 
435  // Current info on edges
436  List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
437 
438  PointEdgeWave<pointEdgePoint> wallCalc
439  (
440  mesh,
441  pp.meshPoints(),
442  wallInfo,
443 
444  allPointInfo,
445  allEdgeInfo,
446  mesh.globalData().nTotalPoints() // max iterations
447  );
448 
449  // Copy edge values into scalarField
450  tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
451  scalarField& edgeDist = tedgeDist.ref();
452 
453  forAll(allEdgeInfo, edgeI)
454  {
455  edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
456  }
457 
458 
459  //{
460  // // For debugging: dump to file
461  // pointScalarField pointDist
462  // (
463  // IOobject
464  // (
465  // "pointDist",
466  // meshRefiner_.timeName(),
467  // mesh.DB(),
468  // IOobject::NO_READ,
469  // IOobject::AUTO_WRITE
470  // ),
471  // pMesh,
472  // dimensionedScalar(dimless, 0)
473  // );
474  //
475  // forAll(allEdgeInfo, edgeI)
476  // {
477  // scalar d = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
478  //
479  // const edge& e = mesh.edges()[edgeI];
480  //
481  // pointDist[e[0]] += d;
482  // pointDist[e[1]] += d;
483  // }
484  // forAll(pointDist, pointi)
485  // {
486  // pointDist[pointi] /= mesh.pointEdges()[pointi].size();
487  // }
488  // Info<< "Writing patch distance to " << pointDist.name()
489  // << " at time " << meshRefiner_.timeName() << endl;
490  //
491  // pointDist.write();
492  //}
493 
494  return tedgeDist;
495 }
496 
497 
498 void Foam::snappySnapDriver::dumpMove
499 (
500  const fileName& fName,
501  const pointField& meshPts,
502  const pointField& surfPts
503 )
504 {
505  // Dump direction of growth into file
506  Info<< "Dumping move direction to " << fName << endl;
507 
508  OFstream nearestStream(fName);
509 
510  label vertI = 0;
511 
512  forAll(meshPts, ptI)
513  {
514  meshTools::writeOBJ(nearestStream, meshPts[ptI]);
515  vertI++;
516 
517  meshTools::writeOBJ(nearestStream, surfPts[ptI]);
518  vertI++;
519 
520  nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
521  }
522 }
523 
524 
525 bool Foam::snappySnapDriver::outwardsDisplacement
526 (
527  const indirectPrimitivePatch& pp,
528  const vectorField& patchDisp
529 )
530 {
531  const vectorField& faceNormals = pp.faceNormals();
532  const labelListList& pointFaces = pp.pointFaces();
533 
534  forAll(pointFaces, pointi)
535  {
536  const labelList& pFaces = pointFaces[pointi];
537 
538  vector disp(patchDisp[pointi]);
539 
540  scalar magDisp = mag(disp);
541 
542  if (magDisp > small)
543  {
544  disp /= magDisp;
545 
546  bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
547 
548  if (!outwards)
549  {
550  Warning<< "Displacement " << patchDisp[pointi]
551  << " at mesh point " << pp.meshPoints()[pointi]
552  << " coord " << pp.points()[pp.meshPoints()[pointi]]
553  << " points through the surrounding patch faces" << endl;
554  return false;
555  }
556  }
557  else
558  {
559  //? Displacement small but in wrong direction. Would probably be ok.
560  }
561  }
562  return true;
563 }
564 
565 
566 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
567 
569 (
570  meshRefinement& meshRefiner,
571  const labelList& globalToMasterPatch,
572  const labelList& globalToSlavePatch
573 )
574 :
575  meshRefiner_(meshRefiner),
576  globalToMasterPatch_(globalToMasterPatch),
577  globalToSlavePatch_(globalToSlavePatch)
578 {}
579 
580 
581 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
582 
584 (
585  const List<labelPair>& baffles
586 )
587 {
588  labelList zonedSurfaces =
589  surfaceZonesInfo::getNamedSurfaces(meshRefiner_.surfaces().surfZones());
590 
592 
593  // No need to sync; all processors will have all same zonedSurfaces.
594  label nBaffles = returnReduce(baffles.size(), sumOp<label>());
595  if (zonedSurfaces.size() && nBaffles > 0)
596  {
597  // Merge any baffles
598  Info<< "Converting " << nBaffles << " baffles back into zoned faces ..."
599  << endl;
600 
601  map = meshRefiner_.mergeBaffles(baffles);
602 
603  Info<< "Converted baffles in = "
604  << meshRefiner_.mesh().time().cpuTimeIncrement()
605  << " s\n" << nl << endl;
606  }
607 
608  return map;
609 }
610 
611 
613 (
614  const fvMesh& mesh,
615  const snapParameters& snapParams,
616  const indirectPrimitivePatch& pp
617 )
618 {
619  const edgeList& edges = pp.edges();
620  const labelListList& pointEdges = pp.pointEdges();
621  const pointField& localPoints = pp.localPoints();
622 
623  scalarField maxEdgeLen(localPoints.size(), -great);
624 
625  forAll(pointEdges, pointi)
626  {
627  const labelList& pEdges = pointEdges[pointi];
628 
629  forAll(pEdges, pEdgeI)
630  {
631  const edge& e = edges[pEdges[pEdgeI]];
632 
633  scalar len = e.mag(localPoints);
634 
635  maxEdgeLen[pointi] = max(maxEdgeLen[pointi], len);
636  }
637  }
638 
640  (
641  mesh,
642  pp.meshPoints(),
643  maxEdgeLen,
644  maxEqOp<scalar>(), // combine op
645  -great // null value
646  );
647 
648  return scalarField(snapParams.snapTol()*maxEdgeLen);
649 }
650 
651 
653 (
654  const meshRefinement& meshRefiner,
655  const snapParameters& snapParams,
656  const label nInitErrors,
657  const List<labelPair>& baffles,
658  motionSmoother& meshMover
659 )
660 {
661  const fvMesh& mesh = meshRefiner.mesh();
662 
663  labelList checkFaces;
664 
665  Info<< "Smoothing patch points ..." << endl;
666  for
667  (
668  label smoothIter = 0;
669  smoothIter < snapParams.nSmoothPatch();
670  smoothIter++
671  )
672  {
673  Info<< "Smoothing iteration " << smoothIter << endl;
674  checkFaces.setSize(mesh.nFaces());
675  forAll(checkFaces, facei)
676  {
677  checkFaces[facei] = facei;
678  }
679 
680  pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
681  // pointField patchDisp
682  //(
683  // smoothLambdaMuPatchDisplacement(meshMover, baffles)
684  //);
685 
686  // The current mesh is the starting mesh to smooth from.
687  meshMover.setDisplacement(patchDisp);
688 
689  meshMover.correct();
690 
691  scalar oldErrorReduction = -1;
692 
693  for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
694  {
695  Info<< nl << "Scaling iteration " << snapIter << endl;
696 
697  if (snapIter == snapParams.nSnap())
698  {
699  Info<< "Displacement scaling for error reduction set to 0."
700  << endl;
701  oldErrorReduction = meshMover.setErrorReduction(0.0);
702  }
703 
704  // Try to adapt mesh to obtain displacement by smoothly
705  // decreasing displacement at error locations.
706  if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
707  {
708  Info<< "Successfully moved mesh" << endl;
709  break;
710  }
711  }
712 
713  if (oldErrorReduction >= 0)
714  {
715  meshMover.setErrorReduction(oldErrorReduction);
716  }
717  Info<< endl;
718  }
719 
720 
721  // The current mesh is the starting mesh to smooth from.
722  meshMover.correct();
723 
724  if (debug&meshRefinement::MESH)
725  {
726  const_cast<Time&>(mesh.time())++;
727  Info<< "Writing patch smoothed mesh to time "
728  << meshRefiner.timeName() << '.' << endl;
729  meshRefiner.write
730  (
733  (
736  ),
737  mesh.time().path()/meshRefiner.timeName()
738  );
739  Info<< "Dumped mesh in = "
740  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
741  }
742 
743  Info<< "Patch points smoothed in = "
744  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
745 }
746 
747 
749 (
750  const fvMesh& mesh,
751  const indirectPrimitivePatch& pp,
752  const word& zoneName
753 )
754 {
755  label zoneI = mesh.faceZones().findZoneID(zoneName);
756 
757  if (zoneI == -1)
758  {
760  << "Cannot find zone " << zoneName
761  << exit(FatalError);
762  }
763 
764  const faceZone& fZone = mesh.faceZones()[zoneI];
765 
766 
767  // Could use PrimitivePatch & localFaces to extract points but might just
768  // as well do it ourselves.
769 
770  boolList pointOnZone(pp.nPoints(), false);
771 
772  forAll(fZone, i)
773  {
774  const face& f = mesh.faces()[fZone[i]];
775 
776  forAll(f, fp)
777  {
778  label meshPointi = f[fp];
779 
781  pp.meshPointMap().find(meshPointi);
782 
783  if (iter != pp.meshPointMap().end())
784  {
785  label pointi = iter();
786  pointOnZone[pointi] = true;
787  }
788  }
789  }
790 
791  return findIndices(pointOnZone, true);
792 }
793 
794 
796 (
797  const fvMesh& mesh,
798  const indirectPrimitivePatch& pp
799 )
800 {
801  const labelListList& pointFaces = pp.pointFaces();
802 
803  tmp<pointField> tavgBoundary
804  (
805  new pointField(pointFaces.size(), Zero)
806  );
807  pointField& avgBoundary = tavgBoundary.ref();
808  labelList nBoundary(pointFaces.size(), 0);
809 
810  forAll(pointFaces, pointi)
811  {
812  const labelList& pFaces = pointFaces[pointi];
813 
814  forAll(pFaces, pfI)
815  {
816  label facei = pFaces[pfI];
817  label meshFacei = pp.addressing()[facei];
818 
819  label own = mesh.faceOwner()[meshFacei];
820  avgBoundary[pointi] +=
821  mesh.cellCentres()[own] - pp.localPoints()[pointi];
822  nBoundary[pointi]++;
823  }
824  }
825 
827  (
828  mesh,
829  pp.meshPoints(),
830  avgBoundary,
831  plusEqOp<point>(), // combine op
832  vector::zero // null value
833  );
835  (
836  mesh,
837  pp.meshPoints(),
838  nBoundary,
839  plusEqOp<label>(), // combine op
840  label(0) // null value
841  );
842 
843  forAll(pointFaces, pointi)
844  {
845  avgBoundary[pointi] =
846  avgBoundary[pointi]/nBoundary[pointi] + pp.localPoints()[pointi];
847  }
848 
849  return tavgBoundary;
850 }
851 
852 
854 (
855  const scalar planarCos,
856  const indirectPrimitivePatch& pp,
857  const pointField& nearestPoint,
858  const vectorField& nearestNormal,
859 
860  vectorField& disp
861 ) const
862 {
863  Info<< "Detecting near surfaces ..." << endl;
864 
865  const pointField& localPoints = pp.localPoints();
866  const labelList& meshPoints = pp.meshPoints();
867  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
868  const fvMesh& mesh = meshRefiner_.mesh();
869 
871  // const scalarField edgeLen(calcEdgeLen(pp));
872  //
875  //
876  //{
877  // const scalar cos45 = Foam::cos(degToRad(45));
878  // vector n(cos45, cos45, cos45);
879  // n /= mag(n);
880  //
881  // pointField start(14*pp.nPoints());
882  // pointField end(start.size());
883  //
884  // label rayI = 0;
885  // forAll(localPoints, pointi)
886  // {
887  // const point& pt = localPoints[pointi];
888  //
889  // // Along coordinate axes
890  //
891  // {
892  // start[rayI] = pt;
893  // point& endPt = end[rayI++];
894  // endPt = pt;
895  // endPt.x() -= edgeLen[pointi];
896  // }
897  // {
898  // start[rayI] = pt;
899  // point& endPt = end[rayI++];
900  // endPt = pt;
901  // endPt.x() += edgeLen[pointi];
902  // }
903  // {
904  // start[rayI] = pt;
905  // point& endPt = end[rayI++];
906  // endPt = pt;
907  // endPt.y() -= edgeLen[pointi];
908  // }
909  // {
910  // start[rayI] = pt;
911  // point& endPt = end[rayI++];
912  // endPt = pt;
913  // endPt.y() += edgeLen[pointi];
914  // }
915  // {
916  // start[rayI] = pt;
917  // point& endPt = end[rayI++];
918  // endPt = pt;
919  // endPt.z() -= edgeLen[pointi];
920  // }
921  // {
922  // start[rayI] = pt;
923  // point& endPt = end[rayI++];
924  // endPt = pt;
925  // endPt.z() += edgeLen[pointi];
926  // }
927  //
928  // // At 45 degrees
929  //
930  // const vector vec(edgeLen[pointi]*n);
931  //
932  // {
933  // start[rayI] = pt;
934  // point& endPt = end[rayI++];
935  // endPt = pt;
936  // endPt.x() += vec.x();
937  // endPt.y() += vec.y();
938  // endPt.z() += vec.z();
939  // }
940  // {
941  // start[rayI] = pt;
942  // point& endPt = end[rayI++];
943  // endPt = pt;
944  // endPt.x() -= vec.x();
945  // endPt.y() += vec.y();
946  // endPt.z() += vec.z();
947  // }
948  // {
949  // start[rayI] = pt;
950  // point& endPt = end[rayI++];
951  // endPt = pt;
952  // endPt.x() += vec.x();
953  // endPt.y() -= vec.y();
954  // endPt.z() += vec.z();
955  // }
956  // {
957  // start[rayI] = pt;
958  // point& endPt = end[rayI++];
959  // endPt = pt;
960  // endPt.x() -= vec.x();
961  // endPt.y() -= vec.y();
962  // endPt.z() += vec.z();
963  // }
964  // {
965  // start[rayI] = pt;
966  // point& endPt = end[rayI++];
967  // endPt = pt;
968  // endPt.x() += vec.x();
969  // endPt.y() += vec.y();
970  // endPt.z() -= vec.z();
971  // }
972  // {
973  // start[rayI] = pt;
974  // point& endPt = end[rayI++];
975  // endPt = pt;
976  // endPt.x() -= vec.x();
977  // endPt.y() += vec.y();
978  // endPt.z() -= vec.z();
979  // }
980  // {
981  // start[rayI] = pt;
982  // point& endPt = end[rayI++];
983  // endPt = pt;
984  // endPt.x() += vec.x();
985  // endPt.y() -= vec.y();
986  // endPt.z() -= vec.z();
987  // }
988  // {
989  // start[rayI] = pt;
990  // point& endPt = end[rayI++];
991  // endPt = pt;
992  // endPt.x() -= vec.x();
993  // endPt.y() -= vec.y();
994  // endPt.z() -= vec.z();
995  // }
996  // }
997  //
998  // labelList surface1;
999  // List<pointIndexHit> hit1;
1000  // labelList region1;
1001  // vectorField normal1;
1002  //
1003  // labelList surface2;
1004  // List<pointIndexHit> hit2;
1005  // labelList region2;
1006  // vectorField normal2;
1007  // surfaces.findNearestIntersection
1008  // (
1009  // unzonedSurfaces, // surfacesToTest,
1010  // start,
1011  // end,
1012  //
1013  // surface1,
1014  // hit1,
1015  // region1,
1016  // normal1,
1017  //
1018  // surface2,
1019  // hit2,
1020  // region2,
1021  // normal2
1022  // );
1023  //
1024  // // All intersections
1025  // {
1026  // OBJstream str
1027  // (
1028  // mesh.time().path()
1029  // / "surfaceHits_" + meshRefiner_.timeName() + ".obj"
1030  // );
1031  //
1032  // Info<< "Dumping intersections with rays to " << str.name()
1033  // << endl;
1034  //
1035  // forAll(hit1, i)
1036  // {
1037  // if (hit1[i].hit())
1038  // {
1039  // str.write(linePointRef(start[i], hit1[i].hitPoint()));
1040  // }
1041  // if (hit2[i].hit())
1042  // {
1043  // str.write(linePointRef(start[i], hit2[i].hitPoint()));
1044  // }
1045  // }
1046  // }
1047  //
1048  // // Co-planar intersections
1049  // {
1050  // OBJstream str
1051  // (
1052  // mesh.time().path()
1053  // / "coplanarHits_" + meshRefiner_.timeName() + ".obj"
1054  // );
1055  //
1056  // Info<< "Dumping intersections with co-planar surfaces to "
1057  // << str.name() << endl;
1058  //
1059  // forAll(localPoints, pointi)
1060  // {
1061  // bool hasNormal = false;
1062  // point surfPointA;
1063  // vector surfNormalA;
1064  // point surfPointB;
1065  // vector surfNormalB;
1066  //
1067  // bool isCoplanar = false;
1068  //
1069  // label rayI = 14*pointi;
1070  // for (label i = 0; i < 14; i++)
1071  // {
1072  // if (hit1[rayI].hit())
1073  // {
1074  // const point& pt = hit1[rayI].hitPoint();
1075  // const vector& n = normal1[rayI];
1076  //
1077  // if (!hasNormal)
1078  // {
1079  // hasNormal = true;
1080  // surfPointA = pt;
1081  // surfNormalA = n;
1082  // }
1083  // else
1084  // {
1085  // if
1086  // (
1087  // meshRefiner_.isGap
1088  // (
1089  // planarCos,
1090  // surfPointA,
1091  // surfNormalA,
1092  // pt,
1093  // n
1094  // )
1095  // )
1096  // {
1097  // isCoplanar = true;
1098  // surfPointB = pt;
1099  // surfNormalB = n;
1100  // break;
1101  // }
1102  // }
1103  // }
1104  // if (hit2[rayI].hit())
1105  // {
1106  // const point& pt = hit2[rayI].hitPoint();
1107  // const vector& n = normal2[rayI];
1108  //
1109  // if (!hasNormal)
1110  // {
1111  // hasNormal = true;
1112  // surfPointA = pt;
1113  // surfNormalA = n;
1114  // }
1115  // else
1116  // {
1117  // if
1118  // (
1119  // meshRefiner_.isGap
1120  // (
1121  // planarCos,
1122  // surfPointA,
1123  // surfNormalA,
1124  // pt,
1125  // n
1126  // )
1127  // )
1128  // {
1129  // isCoplanar = true;
1130  // surfPointB = pt;
1131  // surfNormalB = n;
1132  // break;
1133  // }
1134  // }
1135  // }
1136  //
1137  // rayI++;
1138  // }
1139  //
1140  // if (isCoplanar)
1141  // {
1142  // str.write(linePointRef(surfPointA, surfPointB));
1143  // }
1144  // }
1145  // }
1146  //}
1147 
1148 
1149  const pointField avgCc(avgCellCentres(mesh, pp));
1150 
1151  // Construct rays through localPoints to beyond cell centre
1152  pointField start(pp.nPoints());
1153  pointField end(pp.nPoints());
1154  forAll(localPoints, pointi)
1155  {
1156  const point& pt = localPoints[pointi];
1157  const vector d = 2*(avgCc[pointi]-pt);
1158  start[pointi] = pt - d;
1159  end[pointi] = pt + d;
1160  }
1161 
1162 
1163  autoPtr<OBJstream> gapStr;
1164  if (debug&meshRefinement::ATTRACTION)
1165  {
1166  gapStr.reset
1167  (
1168  new OBJstream
1169  (
1170  mesh.time().path()
1171  / "detectNearSurfaces_" + meshRefiner_.timeName() + ".obj"
1172  )
1173  );
1174  }
1175 
1176 
1177  const PackedBoolList isPatchMasterPoint
1178  (
1180  (
1181  mesh,
1182  meshPoints
1183  )
1184  );
1185 
1186  label nOverride = 0;
1187 
1188  // 1. All points to non-interface surfaces
1189  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1190  {
1191  const labelList unzonedSurfaces =
1193  (
1194  meshRefiner_.surfaces().surfZones()
1195  );
1196 
1197  // Do intersection test
1198  labelList surface1;
1199  List<pointIndexHit> hit1;
1200  labelList region1;
1201  vectorField normal1;
1202 
1203  labelList surface2;
1204  List<pointIndexHit> hit2;
1205  labelList region2;
1206  vectorField normal2;
1207  surfaces.findNearestIntersection
1208  (
1209  unzonedSurfaces,
1210  start,
1211  end,
1212 
1213  surface1,
1214  hit1,
1215  region1,
1216  normal1,
1217 
1218  surface2,
1219  hit2,
1220  region2,
1221  normal2
1222  );
1223 
1224 
1225  forAll(localPoints, pointi)
1226  {
1227  // Current location
1228  const point& pt = localPoints[pointi];
1229 
1230  bool override = false;
1231 
1232  // if (hit1[pointi].hit())
1233  //{
1234  // if
1235  // (
1236  // meshRefiner_.isGap
1237  // (
1238  // planarCos,
1239  // nearestPoint[pointi],
1240  // nearestNormal[pointi],
1241  // hit1[pointi].hitPoint(),
1242  // normal1[pointi]
1243  // )
1244  // )
1245  // {
1246  // disp[pointi] = hit1[pointi].hitPoint()-pt;
1247  // override = true;
1248  // }
1249  //}
1250  // if (hit2[pointi].hit())
1251  //{
1252  // if
1253  // (
1254  // meshRefiner_.isGap
1255  // (
1256  // planarCos,
1257  // nearestPoint[pointi],
1258  // nearestNormal[pointi],
1259  // hit2[pointi].hitPoint(),
1260  // normal2[pointi]
1261  // )
1262  // )
1263  // {
1264  // disp[pointi] = hit2[pointi].hitPoint()-pt;
1265  // override = true;
1266  // }
1267  //}
1268 
1269  if (hit1[pointi].hit() && hit2[pointi].hit())
1270  {
1271  if
1272  (
1273  meshRefiner_.isGap
1274  (
1275  planarCos,
1276  hit1[pointi].hitPoint(),
1277  normal1[pointi],
1278  hit2[pointi].hitPoint(),
1279  normal2[pointi]
1280  )
1281  )
1282  {
1283  // TBD: check if the attraction (to nearest) would attract
1284  // good enough and not override attraction
1285 
1286  if (gapStr.valid())
1287  {
1288  const point& intPt = hit2[pointi].hitPoint();
1289  gapStr().write(linePointRef(pt, intPt));
1290  }
1291 
1292  // Choose hit2 : nearest to end point (so inside the domain)
1293  disp[pointi] = hit2[pointi].hitPoint()-pt;
1294  override = true;
1295  }
1296  }
1297 
1298  if (override && isPatchMasterPoint[pointi])
1299  {
1300  nOverride++;
1301  }
1302  }
1303  }
1304 
1305 
1306  // 2. All points on zones to their respective surface
1307  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1308 
1309  {
1310  // Surfaces with zone information
1311  const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
1312 
1313  const labelList zonedSurfaces = surfaceZonesInfo::getNamedSurfaces
1314  (
1315  surfZones
1316  );
1317 
1318  forAll(zonedSurfaces, i)
1319  {
1320  label zoneSurfI = zonedSurfaces[i];
1321 
1322  const word& faceZoneName = surfZones[zoneSurfI].faceZoneName();
1323 
1324  const labelList surfacesToTest(1, zoneSurfI);
1325 
1326  // Get indices of points both on faceZone and on pp.
1327  labelList zonePointIndices
1328  (
1329  getZoneSurfacePoints
1330  (
1331  mesh,
1332  pp,
1333  faceZoneName
1334  )
1335  );
1336 
1337  // Do intersection test
1338  labelList surface1;
1339  List<pointIndexHit> hit1;
1340  labelList region1;
1341  vectorField normal1;
1342 
1343  labelList surface2;
1344  List<pointIndexHit> hit2;
1345  labelList region2;
1346  vectorField normal2;
1347  surfaces.findNearestIntersection
1348  (
1349  surfacesToTest,
1350  pointField(start, zonePointIndices),
1351  pointField(end, zonePointIndices),
1352 
1353  surface1,
1354  hit1,
1355  region1,
1356  normal1,
1357 
1358  surface2,
1359  hit2,
1360  region2,
1361  normal2
1362  );
1363 
1364 
1365  forAll(hit1, i)
1366  {
1367  label pointi = zonePointIndices[i];
1368 
1369  // Current location
1370  const point& pt = localPoints[pointi];
1371 
1372  bool override = false;
1373 
1374  // if (hit1[i].hit())
1375  //{
1376  // if
1377  // (
1378  // meshRefiner_.isGap
1379  // (
1380  // planarCos,
1381  // nearestPoint[pointi],
1382  // nearestNormal[pointi],
1383  // hit1[i].hitPoint(),
1384  // normal1[i]
1385  // )
1386  // )
1387  // {
1388  // disp[pointi] = hit1[i].hitPoint()-pt;
1389  // override = true;
1390  // }
1391  //}
1392  // if (hit2[i].hit())
1393  //{
1394  // if
1395  // (
1396  // meshRefiner_.isGap
1397  // (
1398  // planarCos,
1399  // nearestPoint[pointi],
1400  // nearestNormal[pointi],
1401  // hit2[i].hitPoint(),
1402  // normal2[i]
1403  // )
1404  // )
1405  // {
1406  // disp[pointi] = hit2[i].hitPoint()-pt;
1407  // override = true;
1408  // }
1409  //}
1410 
1411  if (hit1[i].hit() && hit2[i].hit())
1412  {
1413  if
1414  (
1415  meshRefiner_.isGap
1416  (
1417  planarCos,
1418  hit1[i].hitPoint(),
1419  normal1[i],
1420  hit2[i].hitPoint(),
1421  normal2[i]
1422  )
1423  )
1424  {
1425  if (gapStr.valid())
1426  {
1427  const point& intPt = hit2[i].hitPoint();
1428  gapStr().write(linePointRef(pt, intPt));
1429  }
1430 
1431  disp[pointi] = hit2[i].hitPoint()-pt;
1432  override = true;
1433  }
1434  }
1435 
1436  if (override && isPatchMasterPoint[pointi])
1437  {
1438  nOverride++;
1439  }
1440  }
1441  }
1442  }
1443 
1444  Info<< "Overriding nearest with intersection of close gaps at "
1445  << returnReduce(nOverride, sumOp<label>())
1446  << " out of " << returnReduce(pp.nPoints(), sumOp<label>())
1447  << " points." << endl;
1448 }
1449 
1450 
1453  const meshRefinement& meshRefiner,
1454  const scalarField& snapDist,
1455  const indirectPrimitivePatch& pp,
1456  pointField& nearestPoint,
1457  vectorField& nearestNormal
1458 )
1459 {
1460  Info<< "Calculating patchDisplacement as distance to nearest surface"
1461  << " point ..." << endl;
1462 
1463  const pointField& localPoints = pp.localPoints();
1464  const refinementSurfaces& surfaces = meshRefiner.surfaces();
1465  const fvMesh& mesh = meshRefiner.mesh();
1466 
1467  // Displacement per patch point
1468  vectorField patchDisp(localPoints.size(), Zero);
1469 
1470  if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
1471  {
1472  // Current surface snapped to
1473  labelList snapSurf(localPoints.size(), -1);
1474 
1475  // Divide surfaces into zoned and unzoned
1476  const labelList zonedSurfaces =
1478  (
1479  meshRefiner.surfaces().surfZones()
1480  );
1481  const labelList unzonedSurfaces =
1483  (
1484  meshRefiner.surfaces().surfZones()
1485  );
1486 
1487 
1488  // 1. All points to non-interface surfaces
1489  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1490 
1491  {
1492  List<pointIndexHit> hitInfo;
1493  labelList hitSurface;
1494 
1495  if (nearestNormal.size() == localPoints.size())
1496  {
1497  labelList hitRegion;
1498  vectorField hitNormal;
1499  surfaces.findNearestRegion
1500  (
1501  unzonedSurfaces,
1502  localPoints,
1503  sqr(snapDist),
1504  hitSurface,
1505  hitInfo,
1506  hitRegion,
1507  hitNormal
1508  );
1509 
1510  forAll(hitInfo, pointi)
1511  {
1512  if (hitInfo[pointi].hit())
1513  {
1514  nearestPoint[pointi] = hitInfo[pointi].hitPoint();
1515  nearestNormal[pointi] = hitNormal[pointi];
1516  }
1517  }
1518  }
1519  else
1520  {
1521  surfaces.findNearest
1522  (
1523  unzonedSurfaces,
1524  localPoints,
1525  sqr(snapDist), // sqr of attract distance
1526  hitSurface,
1527  hitInfo
1528  );
1529  }
1530 
1531  forAll(hitInfo, pointi)
1532  {
1533  if (hitInfo[pointi].hit())
1534  {
1535  patchDisp[pointi] =
1536  hitInfo[pointi].hitPoint()
1537  - localPoints[pointi];
1538 
1539  snapSurf[pointi] = hitSurface[pointi];
1540  }
1541  }
1542  }
1543 
1544 
1545 
1546  // 2. All points on zones to their respective surface
1547  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1548 
1549  // Surfaces with zone information
1550  const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
1551 
1552  // Current best snap distance
1553  scalarField minSnapDist(snapDist);
1554 
1555  forAll(zonedSurfaces, i)
1556  {
1557  label zoneSurfI = zonedSurfaces[i];
1558 
1559  const word& faceZoneName = surfZones[zoneSurfI].faceZoneName();
1560 
1561  const labelList surfacesToTest(1, zoneSurfI);
1562 
1563  // Get indices of points both on faceZone and on pp.
1564  labelList zonePointIndices
1565  (
1566  getZoneSurfacePoints
1567  (
1568  mesh,
1569  pp,
1570  faceZoneName
1571  )
1572  );
1573 
1574  // Find nearest for points both on faceZone and pp.
1575  List<pointIndexHit> hitInfo;
1576  labelList hitSurface;
1577 
1578  if (nearestNormal.size() == localPoints.size())
1579  {
1580  labelList hitRegion;
1581  vectorField hitNormal;
1582  surfaces.findNearestRegion
1583  (
1584  surfacesToTest,
1585  pointField(localPoints, zonePointIndices),
1586  sqr(scalarField(minSnapDist, zonePointIndices)),
1587  hitSurface,
1588  hitInfo,
1589  hitRegion,
1590  hitNormal
1591  );
1592 
1593  forAll(hitInfo, i)
1594  {
1595  if (hitInfo[i].hit())
1596  {
1597  label pointi = zonePointIndices[i];
1598  nearestPoint[pointi] = hitInfo[i].hitPoint();
1599  nearestNormal[pointi] = hitNormal[i];
1600  }
1601  }
1602  }
1603  else
1604  {
1605  surfaces.findNearest
1606  (
1607  surfacesToTest,
1608  pointField(localPoints, zonePointIndices),
1609  sqr(scalarField(minSnapDist, zonePointIndices)),
1610  hitSurface,
1611  hitInfo
1612  );
1613  }
1614 
1615  forAll(hitInfo, i)
1616  {
1617  label pointi = zonePointIndices[i];
1618 
1619  if (hitInfo[i].hit())
1620  {
1621  patchDisp[pointi] =
1622  hitInfo[i].hitPoint()
1623  - localPoints[pointi];
1624 
1625  minSnapDist[pointi] = min
1626  (
1627  minSnapDist[pointi],
1628  mag(patchDisp[pointi])
1629  );
1630 
1631  snapSurf[pointi] = zoneSurfI;
1632  }
1633  }
1634  }
1635 
1636  // Check if all points are being snapped
1637  forAll(snapSurf, pointi)
1638  {
1639  if (snapSurf[pointi] == -1)
1640  {
1642  << "For point:" << pointi
1643  << " coordinate:" << localPoints[pointi]
1644  << " did not find any surface within:"
1645  << minSnapDist[pointi]
1646  << " metre." << endl;
1647  }
1648  }
1649 
1650  {
1651  const PackedBoolList isPatchMasterPoint
1652  (
1654  (
1655  mesh,
1656  pp.meshPoints()
1657  )
1658  );
1659 
1660  scalarField magDisp(mag(patchDisp));
1661 
1662  Info<< "Wanted displacement : average:"
1663  << meshRefinement::gAverage(isPatchMasterPoint, magDisp)
1664  << " min:" << gMin(magDisp)
1665  << " max:" << gMax(magDisp) << endl;
1666  }
1667  }
1668 
1669  Info<< "Calculated surface displacement in = "
1670  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1671 
1672 
1673  // Limit amount of movement.
1674  forAll(patchDisp, patchPointi)
1675  {
1676  scalar magDisp = mag(patchDisp[patchPointi]);
1677 
1678  if (magDisp > snapDist[patchPointi])
1679  {
1680  patchDisp[patchPointi] *= snapDist[patchPointi] / magDisp;
1681 
1682  Pout<< "Limiting displacement for " << patchPointi
1683  << " from " << magDisp << " to " << snapDist[patchPointi]
1684  << endl;
1685  }
1686  }
1687 
1688  // Points on zones in one domain but only present as point on other
1689  // will not do condition 2 on all. Sync explicitly.
1691  (
1692  mesh,
1693  pp.meshPoints(),
1694  patchDisp,
1695  minMagSqrEqOp<point>(), // combine op
1696  vector(great, great, great) // null value (note: cannot use vGreat)
1697  );
1698 
1699  return patchDisp;
1700 }
1701 
1702 
1705  const snapParameters& snapParams,
1706  motionSmoother& meshMover
1707 ) const
1708 {
1709  const fvMesh& mesh = meshRefiner_.mesh();
1710  const indirectPrimitivePatch& pp = meshMover.patch();
1711 
1712  Info<< "Smoothing displacement ..." << endl;
1713 
1714  // Set edge diffusivity as inverse of distance to patch
1715  scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + small));
1716  // scalarField edgeGamma(mesh.nEdges(), 1.0);
1717  // scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
1718 
1719  // Get displacement field
1720  pointVectorField& disp = meshMover.displacement();
1721 
1722  for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
1723  {
1724  if ((iter % 10) == 0)
1725  {
1726  Info<< "Iteration " << iter << endl;
1727  }
1728  pointVectorField oldDisp(disp);
1729  meshMover.smooth(oldDisp, edgeGamma, disp);
1730  }
1731  Info<< "Displacement smoothed in = "
1732  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1733 
1734  if (debug&meshRefinement::MESH)
1735  {
1736  const_cast<Time&>(mesh.time())++;
1737  Info<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
1738  << endl;
1739 
1740  // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
1741  // but this will also delete all pointMesh but not pointFields which
1742  // gives an illegal situation.
1743 
1744  meshRefiner_.write
1745  (
1748  (
1751  ),
1752  mesh.time().path()/meshRefiner_.timeName()
1753  );
1754  Info<< "Writing displacement field ..." << endl;
1755  disp.write();
1756  tmp<pointScalarField> magDisp(mag(disp));
1757  magDisp().write();
1758 
1759  Info<< "Writing actual patch displacement ..." << endl;
1760  vectorField actualPatchDisp(disp, pp.meshPoints());
1761  dumpMove
1762  (
1763  mesh.time().path()
1764  / "actualPatchDisplacement_" + meshRefiner_.timeName() + ".obj",
1765  pp.localPoints(),
1766  pp.localPoints() + actualPatchDisp
1767  );
1768  }
1769 }
1770 
1771 
1774  const snapParameters& snapParams,
1775  const label nInitErrors,
1776  const List<labelPair>& baffles,
1777  motionSmoother& meshMover
1778 )
1779 {
1780  const fvMesh& mesh = meshRefiner_.mesh();
1781 
1782  // Relax displacement until correct mesh
1783  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1784  labelList checkFaces(identity(mesh.nFaces()));
1785 
1786  scalar oldErrorReduction = -1;
1787 
1788  bool meshOk = false;
1789 
1790  Info<< "Moving mesh ..." << endl;
1791  for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
1792  {
1793  Info<< nl << "Iteration " << iter << endl;
1794 
1795  if (iter == snapParams.nSnap())
1796  {
1797  Info<< "Displacement scaling for error reduction set to 0." << endl;
1798  oldErrorReduction = meshMover.setErrorReduction(0.0);
1799  }
1800 
1801  meshOk = meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors);
1802 
1803  if (meshOk)
1804  {
1805  Info<< "Successfully moved mesh" << endl;
1806  break;
1807  }
1808  if (debug&meshRefinement::MESH)
1809  {
1810  const_cast<Time&>(mesh.time())++;
1811  Info<< "Writing scaled mesh to time " << meshRefiner_.timeName()
1812  << endl;
1813  mesh.write();
1814 
1815  Info<< "Writing displacement field ..." << endl;
1816  meshMover.displacement().write();
1817  tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
1818  magDisp().write();
1819  }
1820  }
1821 
1822  if (oldErrorReduction >= 0)
1823  {
1824  meshMover.setErrorReduction(oldErrorReduction);
1825  }
1826  Info<< "Moved mesh in = "
1827  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1828 
1829  return meshOk;
1830 }
1831 
1832 
1835  const snapParameters& snapParams,
1836  const labelList& adaptPatchIDs,
1837  const labelList& preserveFaces
1838 )
1839 {
1840  const fvMesh& mesh = meshRefiner_.mesh();
1841  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
1842 
1843  Info<< "Repatching faces according to nearest surface ..." << endl;
1844 
1845  // Get the labels of added patches.
1847  (
1849  (
1850  mesh,
1851  adaptPatchIDs
1852  )
1853  );
1854  indirectPrimitivePatch& pp = ppPtr();
1855 
1856  // Divide surfaces into zoned and unzoned
1857  labelList zonedSurfaces =
1859  labelList unzonedSurfaces =
1861 
1862 
1863  // Faces that do not move
1864  PackedBoolList isZonedFace(mesh.nFaces());
1865  {
1866  // 1. Preserve faces in preserveFaces list
1867  forAll(preserveFaces, facei)
1868  {
1869  if (preserveFaces[facei] != -1)
1870  {
1871  isZonedFace.set(facei, 1);
1872  }
1873  }
1874 
1875  // 2. All faces on zoned surfaces
1876  const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
1877  const meshFaceZones& fZones = mesh.faceZones();
1878 
1879  forAll(zonedSurfaces, i)
1880  {
1881  const label zoneSurfI = zonedSurfaces[i];
1882  const faceZone& fZone = fZones[surfZones[zoneSurfI].faceZoneName()];
1883 
1884  forAll(fZone, i)
1885  {
1886  isZonedFace.set(fZone[i], 1);
1887  }
1888  }
1889  }
1890 
1891 
1892  // Determine per pp face which patch it should be in
1893  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1894 
1895  // Patch that face should be in
1896  labelList closestPatch(pp.size(), -1);
1897  {
1898  // face snap distance as max of point snap distance
1899  scalarField faceSnapDist(pp.size(), -great);
1900  {
1901  // Distance to attract to nearest feature on surface
1902  const scalarField snapDist
1903  (
1904  calcSnapDistance
1905  (
1906  mesh,
1907  snapParams,
1908  pp
1909  )
1910  );
1911 
1912  const faceList& localFaces = pp.localFaces();
1913 
1914  forAll(localFaces, facei)
1915  {
1916  const face& f = localFaces[facei];
1917 
1918  forAll(f, fp)
1919  {
1920  faceSnapDist[facei] = max
1921  (
1922  faceSnapDist[facei],
1923  snapDist[f[fp]]
1924  );
1925  }
1926  }
1927  }
1928 
1929  pointField localFaceCentres(mesh.faceCentres(), pp.addressing());
1930 
1931  // Get nearest surface and region
1932  labelList hitSurface;
1933  labelList hitRegion;
1934  surfaces.findNearestRegion
1935  (
1936  unzonedSurfaces,
1937  localFaceCentres,
1938  sqr(faceSnapDist), // sqr of attract distance
1939  hitSurface,
1940  hitRegion
1941  );
1942 
1943  // Get patch
1944  forAll(pp, i)
1945  {
1946  label facei = pp.addressing()[i];
1947 
1948  if (hitSurface[i] != -1 && !isZonedFace.get(facei))
1949  {
1950  closestPatch[i] = globalToMasterPatch_
1951  [
1952  surfaces.globalRegion
1953  (
1954  hitSurface[i],
1955  hitRegion[i]
1956  )
1957  ];
1958  }
1959  }
1960  }
1961 
1962 
1963  // Change those faces for which there is a different closest patch
1964  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1965 
1966  labelList ownPatch(mesh.nFaces(), -1);
1967  labelList nbrPatch(mesh.nFaces(), -1);
1968 
1969  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1970 
1971  forAll(patches, patchi)
1972  {
1973  const polyPatch& pp = patches[patchi];
1974 
1975  forAll(pp, i)
1976  {
1977  ownPatch[pp.start()+i] = patchi;
1978  nbrPatch[pp.start()+i] = patchi;
1979  }
1980  }
1981 
1982  label nChanged = 0;
1983  forAll(closestPatch, i)
1984  {
1985  label facei = pp.addressing()[i];
1986 
1987  if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[facei])
1988  {
1989  ownPatch[facei] = closestPatch[i];
1990  nbrPatch[facei] = closestPatch[i];
1991  nChanged++;
1992  }
1993  }
1994 
1995  Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
1996  << " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
1997  << endl;
1998 
1999  return meshRefiner_.createBaffles(ownPatch, nbrPatch);
2000 }
2001 
2002 
2003 void Foam::snappySnapDriver::detectWarpedFaces
2004 (
2005  const scalar featureCos,
2006  const indirectPrimitivePatch& pp,
2007 
2008  DynamicList<label>& splitFaces,
2009  DynamicList<labelPair>& splits
2010 ) const
2011 {
2012  const fvMesh& mesh = meshRefiner_.mesh();
2013  const faceList& localFaces = pp.localFaces();
2014  const pointField& localPoints = pp.localPoints();
2015  const labelList& bFaces = pp.addressing();
2016 
2017  splitFaces.clear();
2018  splitFaces.setCapacity(bFaces.size());
2019  splits.clear();
2020  splits.setCapacity(bFaces.size());
2021 
2022  // Determine parallel consistent normals on points
2023  const vectorField pointNormals(PatchTools::pointNormals(mesh, pp));
2024 
2025  face f0(4);
2026  face f1(4);
2027 
2028  forAll(localFaces, facei)
2029  {
2030  const face& f = localFaces[facei];
2031 
2032  if (f.size() >= 4)
2033  {
2034  // See if splitting face across diagonal would make two faces with
2035  // biggish normal angle
2036 
2037  labelPair minDiag(-1, -1);
2038  scalar minCos(great);
2039 
2040  for (label startFp = 0; startFp < f.size()-2; startFp++)
2041  {
2042  label minFp = f.rcIndex(startFp);
2043 
2044  for
2045  (
2046  label endFp = f.fcIndex(f.fcIndex(startFp));
2047  endFp < f.size() && endFp != minFp;
2048  endFp++
2049  )
2050  {
2051  // Form two faces
2052  f0.setSize(endFp-startFp+1);
2053  label i0 = 0;
2054  for (label fp = startFp; fp <= endFp; fp++)
2055  {
2056  f0[i0++] = f[fp];
2057  }
2058  f1.setSize(f.size()+2-f0.size());
2059  label i1 = 0;
2060  for (label fp = endFp; fp != startFp; fp = f.fcIndex(fp))
2061  {
2062  f1[i1++] = f[fp];
2063  }
2064  f1[i1++] = f[startFp];
2065 
2066  // Info<< "Splitting face:" << f << " into f0:" << f0
2067  // << " f1:" << f1 << endl;
2068 
2069  const vector n0 = f0.area(localPoints);
2070  const scalar n0Mag = mag(n0);
2071  const vector n1 = f1.area(localPoints);
2072  const scalar n1Mag = mag(n1);
2073 
2074  if (n0Mag > rootVSmall && n1Mag > rootVSmall)
2075  {
2076  scalar cosAngle = (n0/n0Mag) & (n1/n1Mag);
2077  if (cosAngle < minCos)
2078  {
2079  minCos = cosAngle;
2080  minDiag = labelPair(startFp, endFp);
2081  }
2082  }
2083  }
2084  }
2085 
2086 
2087  if (minCos < featureCos)
2088  {
2089  splitFaces.append(bFaces[facei]);
2090  splits.append(minDiag);
2091  }
2092  }
2093  }
2094 }
2095 
2096 
2099  const dictionary& snapDict,
2100  const dictionary& motionDict,
2101  const scalar featureCos,
2102  const scalar planarAngle,
2103  const snapParameters& snapParams
2104 )
2105 {
2106  fvMesh& mesh = meshRefiner_.mesh();
2107 
2108  Info<< nl
2109  << "Morphing phase" << nl
2110  << "--------------" << nl
2111  << endl;
2112 
2113  // Get the labels of added patches.
2114  labelList adaptPatchIDs(meshRefiner_.meshedPatches());
2115 
2116 
2117  // faceZone handling
2118  // ~~~~~~~~~~~~~~~~~
2119  //
2120  // We convert all faceZones into baffles during snapping so we can use
2121  // a standard mesh motion (except for the mesh checking which for baffles
2122  // created from internal faces should check across the baffles). The state
2123  // is stored in two variables:
2124  // baffles : pairs of boundary faces
2125  // duplicateFace : from mesh face to its baffle colleague (or -1 for
2126  // normal faces)
2127  // There are three types of faceZones according to the faceType property:
2128  //
2129  // internal
2130  // --------
2131  // - baffles: contains all faces on faceZone so
2132  // - mesh checks check across baffles
2133  // - they get back merged into internal faces
2134  // - duplicateFace: from face to duplicate face. Contains
2135  // all faces on faceZone to prevents merging patch faces.
2136  //
2137  // baffle
2138  // ------
2139  // - baffles: contains no faces on faceZone since need not be merged/checked
2140  // across
2141  // - duplicateFace: contains all faces on faceZone to prevent
2142  // merging patch faces.
2143  //
2144  // boundary
2145  // --------
2146  // - baffles: contains no faces on faceZone since need not be merged/checked
2147  // across
2148  // - duplicateFace: contains no faces on faceZone since both sides can
2149  // merge faces independently.
2150 
2151 
2152  // Create baffles (pairs of faces that share the same points)
2153  // Baffles stored as owner and neighbour face that have been created.
2154  List<labelPair> baffles;
2155  meshRefiner_.createZoneBaffles
2156  (
2157  globalToMasterPatch_,
2158  globalToSlavePatch_,
2159  baffles
2160  );
2161 
2162  // Maintain map from face to baffle face (-1 for non-baffle faces). Used
2163  // later on to prevent patchface merging if faceType=baffle
2164  labelList duplicateFace(mesh.nFaces(), -1);
2165 
2166  forAll(baffles, i)
2167  {
2168  const labelPair& baffle = baffles[i];
2169  duplicateFace[baffle.first()] = baffle.second();
2170  duplicateFace[baffle.second()] = baffle.first();
2171  }
2172 
2173  // Selectively 'forget' about the baffles, i.e. not check across them
2174  // or merge across them.
2175  {
2176  const meshFaceZones& fZones = mesh.faceZones();
2177  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
2178  const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
2179 
2180  // Determine which
2181  // - faces to remove from list of baffles (so not merge)
2182  // - points to duplicate
2183 
2184  // Per face if is on faceType 'baffle' or 'boundary'
2185  labelList filterFace(mesh.nFaces(), -1);
2186  label nFilterFaces = 0;
2187  // Per point whether it need to be duplicated
2188  PackedBoolList duplicatePoint(mesh.nPoints());
2189  label nDuplicatePoints = 0;
2190  forAll(surfZones, surfI)
2191  {
2192  const word& faceZoneName = surfZones[surfI].faceZoneName();
2193 
2194  if (faceZoneName.size())
2195  {
2196  const surfaceZonesInfo::faceZoneType& faceType =
2197  surfZones[surfI].faceType();
2198 
2199  if
2200  (
2201  faceType == surfaceZonesInfo::BAFFLE
2202  || faceType == surfaceZonesInfo::BOUNDARY
2203  )
2204  {
2205  // Filter out all faces for this zone.
2206  label zoneI = fZones.findZoneID(faceZoneName);
2207  const faceZone& fZone = fZones[zoneI];
2208  forAll(fZone, i)
2209  {
2210  label facei = fZone[i];
2211  filterFace[facei] = zoneI;
2212  nFilterFaces++;
2213  }
2214 
2215  if (faceType == surfaceZonesInfo::BOUNDARY)
2216  {
2217  forAll(fZone, i)
2218  {
2219  label facei = fZone[i];
2220 
2221  // Allow combining patch faces across this face
2222  duplicateFace[facei] = -1;
2223 
2224  const face& f = mesh.faces()[facei];
2225  forAll(f, fp)
2226  {
2227  if (!duplicatePoint[f[fp]])
2228  {
2229  duplicatePoint[f[fp]] = 1;
2230  nDuplicatePoints++;
2231  }
2232  }
2233  }
2234  }
2235 
2236  Info<< "Surface : " << surfaces.names()[surfI] << nl
2237  << " faces to become baffle : "
2238  << returnReduce(nFilterFaces, sumOp<label>()) << nl
2239  << " points to duplicate : "
2240  << returnReduce(nDuplicatePoints, sumOp<label>())
2241  << endl;
2242  }
2243  }
2244  }
2245 
2246  // Duplicate points if any processor says it needs duplication
2248  (
2249  mesh,
2250  duplicatePoint,
2251  orEqOp<unsigned int>(), // combine op
2252  0u // null value
2253  );
2254  // Mark as duplicate (avoids combining patch faces) if one or both
2255  syncTools::syncFaceList(mesh, duplicateFace, maxEqOp<label>());
2256  // Mark as resulting from baffle/boundary face zone only if both agree
2257  syncTools::syncFaceList(mesh, filterFace, minEqOp<label>());
2258 
2259  // Duplicate points
2260  if (returnReduce(nDuplicatePoints, sumOp<label>()) > 0)
2261  {
2262  // Collect all points (recount since syncPointList might have
2263  // increased set)
2264  nDuplicatePoints = 0;
2265  forAll(duplicatePoint, pointi)
2266  {
2267  if (duplicatePoint[pointi])
2268  {
2269  nDuplicatePoints++;
2270  }
2271  }
2272  labelList candidatePoints(nDuplicatePoints);
2273  nDuplicatePoints = 0;
2274  forAll(duplicatePoint, pointi)
2275  {
2276  if (duplicatePoint[pointi])
2277  {
2278  candidatePoints[nDuplicatePoints++] = pointi;
2279  }
2280  }
2281 
2282 
2283  localPointRegion regionSide(mesh, candidatePoints);
2284  autoPtr<mapPolyMesh> mapPtr = meshRefiner_.dupNonManifoldPoints
2285  (
2286  regionSide
2287  );
2289  (
2290  mapPtr().faceMap(),
2291  label(-1),
2292  filterFace
2293  );
2295  (
2296  mapPtr().faceMap(),
2297  label(-1),
2298  duplicateFace
2299  );
2300 
2301  // Update baffles and baffle-to-baffle addressing
2302 
2303  const labelList& reverseFaceMap = mapPtr().reverseFaceMap();
2304 
2305  forAll(baffles, i)
2306  {
2307  labelPair& baffle = baffles[i];
2308  baffle.first() = reverseFaceMap[baffle.first()];
2309  baffle.second() = reverseFaceMap[baffle.second()];
2310  }
2311 
2312  if (debug&meshRefinement::MESH)
2313  {
2314  const_cast<Time&>(mesh.time())++;
2315  Pout<< "Writing duplicatedPoints mesh to time "
2316  << meshRefiner_.timeName()
2317  << endl;
2318  meshRefiner_.write
2319  (
2322  (
2325  ),
2326  mesh.time().path()/"duplicatedPoints"
2327  );
2328  }
2329  }
2330 
2331 
2332  // Forget about baffles in a BAFFLE/BOUNDARY type zone
2333  DynamicList<labelPair> newBaffles(baffles.size());
2334  forAll(baffles, i)
2335  {
2336  const labelPair& baffle = baffles[i];
2337  if
2338  (
2339  filterFace[baffle.first()] == -1
2340  && filterFace[baffles[i].second()] == -1
2341  )
2342  {
2343  newBaffles.append(baffle);
2344  }
2345  }
2346 
2347  if (newBaffles.size() < baffles.size())
2348  {
2349  // Info<< "Splitting baffles into" << nl
2350  // << " internal : " << newBaffles.size() << nl
2351  // << " baffle : " << baffles.size()-newBaffles.size()
2352  // << nl << endl;
2353  baffles.transfer(newBaffles);
2354  }
2355  Info<< endl;
2356  }
2357 
2358 
2359  bool doFeatures = false;
2360  label nFeatIter = 1;
2361  if (snapParams.nFeatureSnap() > 0)
2362  {
2363  doFeatures = true;
2364  nFeatIter = snapParams.nFeatureSnap();
2365 
2366  Info<< "Snapping to features in " << nFeatIter
2367  << " iterations ..." << endl;
2368  }
2369 
2370 
2371  bool meshOk = false;
2372 
2373  {
2375  (
2377  (
2378  mesh,
2379  adaptPatchIDs
2380  )
2381  );
2382 
2383 
2384  // Distance to attract to nearest feature on surface
2385  const scalarField snapDist(calcSnapDistance(mesh, snapParams, ppPtr()));
2386 
2387 
2388  // Construct iterative mesh mover.
2389  Info<< "Constructing mesh displacer ..." << endl;
2390  Info<< "Using mesh parameters " << motionDict << nl << endl;
2391 
2392  autoPtr<motionSmoother> meshMoverPtr
2393  (
2394  new motionSmoother
2395  (
2396  mesh,
2397  ppPtr(),
2398  adaptPatchIDs,
2400  (
2401  pointMesh::New(mesh),
2402  adaptPatchIDs
2403  ),
2404  motionDict
2405  )
2406  );
2407 
2408 
2409  // Check initial mesh
2410  Info<< "Checking initial mesh ..." << endl;
2411  labelHashSet wrongFaces(mesh.nFaces()/100);
2412  motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
2413  const label nInitErrors = returnReduce
2414  (
2415  wrongFaces.size(),
2416  sumOp<label>()
2417  );
2418 
2419  Info<< "Detected " << nInitErrors << " illegal faces"
2420  << " (concave, zero area or negative cell pyramid volume)"
2421  << endl;
2422 
2423 
2424  Info<< "Checked initial mesh in = "
2425  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
2426 
2427  // Pre-smooth patch vertices (so before determining nearest)
2428  preSmoothPatch
2429  (
2430  meshRefiner_,
2431  snapParams,
2432  nInitErrors,
2433  baffles,
2434  meshMoverPtr()
2435  );
2436 
2437 
2438 
2439  //- Only if in feature attraction mode:
2440  // Nearest feature
2441  vectorField patchAttraction;
2442  // Constraints at feature
2443  List<pointConstraint> patchConstraints;
2444 
2445 
2446  for (label iter = 0; iter < nFeatIter; iter++)
2447  {
2448  // if (doFeatures && (iter == 0 || iter == nFeatIter/2))
2449  //{
2450  // Info<< "Splitting diagonal attractions" << endl;
2451  //
2452  // indirectPrimitivePatch& pp = ppPtr();
2453  // motionSmoother& meshMover = meshMoverPtr();
2454  //
2455  // // Calculate displacement at every patch point. Insert into
2456  // // meshMover.
2457  // // Calculate displacement at every patch point
2458  // pointField nearestPoint;
2459  // vectorField nearestNormal;
2460  //
2461  // if (snapParams.detectNearSurfacesSnap())
2462  // {
2463  // nearestPoint.setSize(pp.nPoints(), vector::max);
2464  // nearestNormal.setSize(pp.nPoints(), Zero);
2465  // }
2466  //
2467  // vectorField disp = calcNearestSurface
2468  // (
2469  // meshRefiner_,
2470  // snapDist,
2471  // pp,
2472  // nearestPoint,
2473  // nearestNormal
2474  // );
2475  //
2476  //
2477  // // Override displacement at thin gaps
2478  // if (snapParams.detectNearSurfacesSnap())
2479  // {
2480  // detectNearSurfaces
2481  // (
2482  // Foam::cos(degToRad(planarAngle)),// planar gaps
2483  // pp,
2484  // nearestPoint, // surfacepoint from nearest test
2485  // nearestNormal, // surfacenormal from nearest test
2486  //
2487  // disp
2488  // );
2489  // }
2490  //
2491  // // Override displacement with feature edge attempt
2492  // const label iter = 0;
2493  // calcNearestSurfaceFeature
2494  // (
2495  // snapParams,
2496  // false, // avoidSnapProblems
2497  // iter,
2498  // featureCos,
2499  // scalar(iter+1)/nFeatIter,
2500  // snapDist,
2501  // disp,
2502  // meshMover,
2503  // patchAttraction,
2504  // patchConstraints
2505  // );
2506  //
2507  //
2508  // const labelList& bFaces = ppPtr().addressing();
2509  // DynamicList<label> splitFaces(bFaces.size());
2510  // DynamicList<labelPair> splits(bFaces.size());
2511  //
2512  // forAll(bFaces, facei)
2513  // {
2514  // const labelPair split
2515  // (
2516  // findDiagonalAttraction
2517  // (
2518  // ppPtr(),
2519  // patchAttraction,
2520  // patchConstraints,
2521  // facei
2522  // )
2523  // );
2524  //
2525  // if (split != labelPair(-1, -1))
2526  // {
2527  // splitFaces.append(bFaces[facei]);
2528  // splits.append(split);
2529  // }
2530  // }
2531  //
2532  // Info<< "Splitting "
2533  // << returnReduce(splitFaces.size(), sumOp<label>())
2534  // << " faces along diagonal attractions" << endl;
2535  //
2536  // autoPtr<mapPolyMesh> mapPtr = meshRefiner_.splitFaces
2537  // (
2538  // splitFaces,
2539  // splits
2540  // );
2541  //
2542  // const labelList& faceMap = mapPtr().faceMap();
2543  // meshRefinement::updateList(faceMap, -1, duplicateFace);
2544  // const labelList& reverseFaceMap = mapPtr().reverseFaceMap();
2545  // forAll(baffles, i)
2546  // {
2547  // labelPair& baffle = baffles[i];
2548  // baffle.first() = reverseFaceMap[baffle.first()];
2549  // baffle.second() = reverseFaceMap[baffle.second()];
2550  // }
2551  //
2552  // meshMoverPtr.clear();
2553  // ppPtr.clear();
2554  //
2555  // ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
2556  // meshMoverPtr.reset
2557  // (
2558  // new motionSmoother
2559  // (
2560  // mesh,
2561  // ppPtr(),
2562  // adaptPatchIDs,
2563  // meshRefinement::makeDisplacementField
2564  // (
2565  // pointMesh::New(mesh),
2566  // adaptPatchIDs
2567  // ),
2568  // motionDict
2569  // )
2570  // );
2571  //
2572  // if (debug&meshRefinement::MESH)
2573  // {
2574  // const_cast<Time&>(mesh.time())++;
2575  // Info<< "Writing split diagonal mesh to time "
2576  // << meshRefiner_.timeName() << endl;
2577  // meshRefiner_.write
2578  // (
2579  // meshRefinement::debugType(debug),
2580  // meshRefinement::writeType
2581  // (
2582  // meshRefinement::writeLevel()
2583  // | meshRefinement::WRITEMESH
2584  // ),
2585  // mesh.time().path()/meshRefiner_.timeName()
2586  // );
2587  // }
2588  //}
2589  // else
2590  // if
2591  //(
2592  // doFeatures
2593  // && (iter == 1 || iter == nFeatIter/2+1 || iter == nFeatIter-1)
2594  //)
2595  //{
2596  // Info<< "Splitting warped faces" << endl;
2597  //
2598  // const labelList& bFaces = ppPtr().addressing();
2599  // DynamicList<label> splitFaces(bFaces.size());
2600  // DynamicList<labelPair> splits(bFaces.size());
2601  //
2602  // detectWarpedFaces
2603  // (
2604  // featureCos,
2605  // ppPtr(),
2606  //
2607  // splitFaces,
2608  // splits
2609  // );
2610  //
2611  // Info<< "Splitting "
2612  // << returnReduce(splitFaces.size(), sumOp<label>())
2613  // << " faces along diagonal to avoid warpage" << endl;
2614  //
2615  // autoPtr<mapPolyMesh> mapPtr = meshRefiner_.splitFaces
2616  // (
2617  // splitFaces,
2618  // splits
2619  // );
2620  //
2621  // const labelList& faceMap = mapPtr().faceMap();
2622  // meshRefinement::updateList(faceMap, -1, duplicateFace);
2623  // const labelList& reverseFaceMap = mapPtr().reverseFaceMap();
2624  // forAll(baffles, i)
2625  // {
2626  // labelPair& baffle = baffles[i];
2627  // baffle.first() = reverseFaceMap[baffle.first()];
2628  // baffle.second() = reverseFaceMap[baffle.second()];
2629  // }
2630  //
2631  // meshMoverPtr.clear();
2632  // ppPtr.clear();
2633  //
2634  // ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
2635  // meshMoverPtr.reset
2636  // (
2637  // new motionSmoother
2638  // (
2639  // mesh,
2640  // ppPtr(),
2641  // adaptPatchIDs,
2642  // meshRefinement::makeDisplacementField
2643  // (
2644  // pointMesh::New(mesh),
2645  // adaptPatchIDs
2646  // ),
2647  // motionDict
2648  // )
2649  // );
2650  //
2651  // if (debug&meshRefinement::MESH)
2652  // {
2653  // const_cast<Time&>(mesh.time())++;
2654  // Info<< "Writing split warped mesh to time "
2655  // << meshRefiner_.timeName() << endl;
2656  // meshRefiner_.write
2657  // (
2658  // meshRefinement::debugType(debug),
2659  // meshRefinement::writeType
2660  // (
2661  // meshRefinement::writeLevel()
2662  // | meshRefinement::WRITEMESH
2663  // ),
2664  // mesh.time().path()/meshRefiner_.timeName()
2665  // );
2666  // }
2667  //}
2668 
2669 
2670 
2671  Info<< nl
2672  << "Morph iteration " << iter << nl
2673  << "-----------------" << endl;
2674 
2675  indirectPrimitivePatch& pp = ppPtr();
2676  motionSmoother& meshMover = meshMoverPtr();
2677 
2678 
2679  // Calculate displacement at every patch point. Insert into
2680  // meshMover.
2681  // Calculate displacement at every patch point
2682  pointField nearestPoint;
2683  vectorField nearestNormal;
2684 
2685  if (snapParams.detectNearSurfacesSnap())
2686  {
2687  nearestPoint.setSize(pp.nPoints(), vector::max);
2688  nearestNormal.setSize(pp.nPoints(), Zero);
2689  }
2690 
2691  vectorField disp = calcNearestSurface
2692  (
2693  meshRefiner_,
2694  snapDist,
2695  pp,
2696  nearestPoint,
2697  nearestNormal
2698  );
2699 
2700 
2701  // Override displacement at thin gaps
2702  if (snapParams.detectNearSurfacesSnap())
2703  {
2704  detectNearSurfaces
2705  (
2706  Foam::cos(degToRad(planarAngle)),// planar cos for gaps
2707  pp,
2708  nearestPoint, // surfacepoint from nearest test
2709  nearestNormal, // surfacenormal from nearest test
2710 
2711  disp
2712  );
2713  }
2714 
2715  // Override displacement with feature edge attempt
2716  if (doFeatures)
2717  {
2718  disp = calcNearestSurfaceFeature
2719  (
2720  snapParams,
2721  true, // avoidSnapProblems
2722  iter,
2723  featureCos,
2724  scalar(iter+1)/nFeatIter,
2725  snapDist,
2726  disp,
2727  meshMover,
2728  patchAttraction,
2729  patchConstraints
2730  );
2731  }
2732 
2733  // Check for displacement being outwards.
2734  outwardsDisplacement(pp, disp);
2735 
2736  // Set initial distribution of displacement field (on patches)
2737  // from patchDisp and make displacement consistent with b.c.
2738  // on displacement pointVectorField.
2739  meshMover.setDisplacement(disp);
2740 
2741 
2742  if (debug&meshRefinement::ATTRACTION)
2743  {
2744  dumpMove
2745  (
2746  mesh.time().path()
2747  / "patchDisplacement_" + name(iter) + ".obj",
2748  pp.localPoints(),
2749  pp.localPoints() + disp
2750  );
2751  }
2752 
2753  // Get smoothly varying internal displacement field.
2754  smoothDisplacement(snapParams, meshMover);
2755 
2756  // Apply internal displacement to mesh.
2757  meshOk = scaleMesh
2758  (
2759  snapParams,
2760  nInitErrors,
2761  baffles,
2762  meshMover
2763  );
2764 
2765  if (!meshOk)
2766  {
2768  << "Did not successfully snap mesh."
2769  << " Continuing to snap to resolve easy" << nl
2770  << " surfaces but the"
2771  << " resulting mesh will not satisfy your quality"
2772  << " constraints" << nl << endl;
2773  }
2774 
2775  if (debug&meshRefinement::MESH)
2776  {
2777  const_cast<Time&>(mesh.time())++;
2778  Info<< "Writing scaled mesh to time "
2779  << meshRefiner_.timeName() << endl;
2780  meshRefiner_.write
2781  (
2784  (
2787  ),
2788  mesh.time().path()/meshRefiner_.timeName()
2789  );
2790  Info<< "Writing displacement field ..." << endl;
2791  meshMover.displacement().write();
2792  tmp<pointScalarField> magDisp
2793  (
2794  mag(meshMover.displacement())
2795  );
2796  magDisp().write();
2797  }
2798 
2799  // Use current mesh as base mesh
2800  meshMover.correct();
2801  }
2802  }
2803 
2804 
2805  // Merge any introduced baffles (from faceZones of faceType 'internal')
2806  {
2807  autoPtr<mapPolyMesh> mapPtr = mergeZoneBaffles(baffles);
2808 
2809  if (mapPtr.valid())
2810  {
2811  forAll(duplicateFace, facei)
2812  {
2813  if (duplicateFace[facei] != -1)
2814  {
2815  duplicateFace[facei] = mapPtr().reverseFaceMap()[facei];
2816  }
2817  }
2818  }
2819  }
2820 
2821  // Repatch faces according to nearest. Do not repatch baffle faces.
2822  {
2823  autoPtr<mapPolyMesh> mapPtr = repatchToSurface
2824  (
2825  snapParams,
2826  adaptPatchIDs,
2827  duplicateFace
2828  );
2830  (
2831  mapPtr().faceMap(),
2832  label(-1),
2833  duplicateFace
2834  );
2835  }
2836 
2837  // Repatching might have caused faces to be on same patch and hence
2838  // mergeable so try again to merge coplanar faces. Do not merge baffle
2839  // faces to ensure they both stay the same.
2840  label nChanged = meshRefiner_.mergePatchFacesUndo
2841  (
2842  featureCos, // minCos
2843  featureCos, // concaveCos
2844  meshRefiner_.meshedPatches(),
2845  motionDict,
2846  duplicateFace // faces not to merge
2847  );
2848 
2849  nChanged += meshRefiner_.mergeEdgesUndo(featureCos, motionDict);
2850 
2851  if (nChanged > 0 && debug & meshRefinement::MESH)
2852  {
2853  const_cast<Time&>(mesh.time())++;
2854  Info<< "Writing patchFace merged mesh to time "
2855  << meshRefiner_.timeName() << endl;
2856  meshRefiner_.write
2857  (
2860  (
2863  ),
2864  meshRefiner_.timeName()
2865  );
2866  }
2867 }
2868 
2869 
2870 // ************************************************************************* //
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
Given a displacement moves the mesh by scaling the displacement back until there are no more mesh err...
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
const labelListList & pointEdges() const
Return point-edge addressing.
Switch detectNearSurfacesSnap() const
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
label nPoints() const
Return number of points supporting patch faces.
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: MeshZones.C:341
const refinementSurfaces & surfaces() const
Reference to surface search engines.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
label nSmoothDispl() const
Number of mesh displacement smoothing iterations.
virtual Ostream & write(const char)=0
Write character.
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
const labelList & surfaces() const
bool scaleMesh(const snapParameters &snapParams, const label nInitErrors, const List< labelPair > &baffles, motionSmoother &)
Do the hard work: move the mesh according to displacement,.
snappySnapDriver(meshRefinement &meshRefiner, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch)
Construct from components.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
void reset(T *=nullptr)
If object pointer already set, delete object and set to given.
Definition: autoPtrI.H:114
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
Type gMin(const FieldField< Field, Type > &f)
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
static void setDisplacement(const labelList &patchIDs, const indirectPrimitivePatch &pp, pointField &patchDisp, pointVectorField &displacement)
Set displacement field from displacement on patch points.
autoPtr< mapPolyMesh > mergeZoneBaffles(const List< labelPair > &)
Merge baffles.
label nFaces() const
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
bool scaleMesh(labelList &checkFaces, const bool smoothMesh=true, const label nAllow=0)
Move mesh with given scale. Return true if mesh ok or has.
void doSnap(const dictionary &snapDict, const dictionary &motionDict, const scalar featureCos, const scalar planarAngle, const snapParameters &snapParams)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const pointMesh & pMesh() const
Reference to pointMesh.
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
List< face > faceList
Definition: faceListFwd.H:43
label nSmoothPatch() const
Number of patch smoothing iterations before finding.
void detectNearSurfaces(const scalar planarCos, const indirectPrimitivePatch &, const pointField &nearestPoint, const vectorField &nearestNormal, vectorField &disp) const
Per patch point override displacement if in gap situation.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static writeType writeLevel()
Get/set write level.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
static vector area(const PointField &ps)
Return vector area given face points.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const Field< PointType > & localPoints() const
Return pointField of points in patch.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:239
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurrences of given element. Linear search.
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
virtual bool write(const bool write=true) const
Write mesh using IO settings from time.
Definition: fvMesh.C:1093
void findNearest(const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &) const
Find nearest point on surfaces.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
label nFeatureSnap() const
static PackedBoolList getMasterFaces(const polyMesh &)
Get per face whether it is uncoupled or a master of a.
Definition: syncTools.C:153
scalar f1
Definition: createFields.H:15
static labelList getZoneSurfacePoints(const fvMesh &mesh, const indirectPrimitivePatch &, const word &zoneName)
Get points both on patch and facezone.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Takes mesh with &#39;baffles&#39; (= boundary faces sharing points). Determines for selected points on bounda...
const fvMesh & mesh() const
Reference to mesh.
Container for data on surfaces used for surface-driven refinement. Contains all the data about the le...
void correct()
Take over existing mesh position.
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
static void updateList(const labelList &newToOld, const T &nullValue, List< T > &elems)
Helper: reorder list according to map.
A list of faces which address into the list of points.
Determines the &#39;side&#39; for every face and connected to a singly-connected (through edges) region of fa...
Definition: regionSide.H:61
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
scalar setErrorReduction(const scalar)
Set the errorReduction (by how much to scale the displacement.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
dimensionedScalar cos(const dimensionedScalar &ds)
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
void findNearestRegion(const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, labelList &hitSurface, labelList &hitRegion) const
Find nearest point on surfaces. Return surface and region on.
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Simple container to keep together snap specific information.
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:74
autoPtr< mapPolyMesh > repatchToSurface(const snapParameters &snapParams, const labelList &adaptPatchIDs, const labelList &preserveFaces)
Repatch faces according to surface nearest the face centre.
void transfer(const FixedList< T, Size > &)
Copy (not transfer) the argument contents.
Definition: FixedListI.H:191
void setCapacity(const label)
Alter the size of the underlying storage.
Definition: DynamicListI.H:130
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
void smoothDisplacement(const snapParameters &snapParams, motionSmoother &) const
Smooth the displacement field to the internal.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static void preSmoothPatch(const meshRefinement &meshRefiner, const snapParameters &snapParams, const label nInitErrors, const List< labelPair > &baffles, motionSmoother &)
Smooth the mesh (patch and internal) to increase visibility.
static const zero Zero
Definition: zero.H:97
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: edgeI.H:181
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
const PtrList< surfaceZonesInfo > & surfZones() const
const vectorField & cellCentres() const
bool visNormal(const vector &n, const vectorField &faceNormals, const labelList &faceLabels)
Check if n is in same direction as normals of all faceLabels.
Definition: meshTools.C:33
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void smooth(const GeometricField< Type, pointPatchField, pointMesh > &fld, const scalarField &edgeWeight, GeometricField< Type, pointPatchField, pointMesh > &newFld) const
Fully explicit smoothing of fields (not positions)
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
static scalarField calcSnapDistance(const fvMesh &mesh, const snapParameters &snapParams, const indirectPrimitivePatch &)
Calculate edge length per patch point.
Foam::polyBoundaryMesh.
const labelListList & pointFaces() const
Return point-face addressing.
OFstream which keeps track of vertices.
Definition: OBJstream.H:53
static const char nl
Definition: Ostream.H:260
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
const indirectPrimitivePatch & patch() const
Reference to patch.
static tmp< pointField > avgCellCentres(const fvMesh &mesh, const indirectPrimitivePatch &)
Helper: calculate average cell centre per point.
Merge points. See below.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static autoPtr< indirectPrimitivePatch > makePatch(const polyMesh &, const labelList &)
Create patch from set of patches.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
messageStream Warning
const Type & second() const
Return second.
Definition: Pair.H:110
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
static tmp< pointField > pointNormals(const polyMesh &, const PrimitivePatch< FaceList, PointField > &)
Return parallel consistent point normals for patches using mesh points.
const vectorField & faceCentres() const
label nSnap() const
Maximum number of snapping relaxation iterations. Should stop.
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronise values on all mesh faces.
Definition: syncTools.H:381
void setSize(const label)
Reset size of List.
Definition: List.C:281
A bit-packed bool list.
const wordList & names() const
Names of surfaces.
scalar snapTol() const
Relative distance for points to be attracted by surface.
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
#define WarningInFunction
Report a warning using Foam::Warning.
label newPointi
Definition: readKivaGrid.H:501
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:476
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:309
label mergePoints(const UList< Type > &points, const scalar mergeTol, const bool verbose, labelList &pointMap, const Type &origin=Type::zero)
Sorts and merges points. All points closer than/equal mergeTol get merged.
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh points.
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
messageStream Info
SubField< vector > subField
Declare type of subField.
Definition: Field.H:100
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
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)
faceZoneType
What to do with faceZone faces.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
bool write() const
Write mesh and all data.
virtual bool write(const bool write=true) const
Write using setting from DB.
static tmp< pointVectorField > makeDisplacementField(const pointMesh &pMesh, const labelList &adaptPatchIDs)
Helper function to make a pointVectorField with correct.
A class for managing temporary objects.
Definition: PtrList.H:53
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
static PackedBoolList getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
void findNearestIntersection(const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surface1, List< pointIndexHit > &hit1, labelList &region1, labelList &surface2, List< pointIndexHit > &hit2, labelList &region2) const
Find intersection nearest to the endpoints. surface1,2 are.
fileName path() const
Return path.
Definition: TimePaths.H:139
const Type & first() const
Return first.
Definition: Pair.H:98
Namespace for OpenFOAM.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49
pointVectorField & displacement()
Reference to displacement field.
static vectorField calcNearestSurface(const meshRefinement &meshRefiner, const scalarField &snapDist, const indirectPrimitivePatch &, pointField &nearestPoint, vectorField &nearestNormal)
Per patch point calculate point on nearest surface. Set as.