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