FreeFOAM The Cross-Platform CFD Toolkit
autoSnapDriver.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) 1991-2010 OpenCFD Ltd.
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 "autoSnapDriver.H"
30 #include <OpenFOAM/Time.H>
33 #include <OpenFOAM/OFstream.H>
34 #include <OpenFOAM/syncTools.H>
35 #include <finiteVolume/fvMesh.H>
36 #include <OpenFOAM/Time.H>
37 #include <OpenFOAM/OFstream.H>
38 #include <OpenFOAM/mapPolyMesh.H>
42 #include <OpenFOAM/mergePoints.H>
45 
46 
47 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
48 
49 namespace Foam
50 {
51 
52 defineTypeNameAndDebug(autoSnapDriver, 0);
53 
54 } // End namespace Foam
55 
56 
57 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
58 
59 // Get faces to repatch. Returns map from face to patch.
60 Foam::Map<Foam::label> Foam::autoSnapDriver::getZoneBafflePatches
61 (
62  const bool allowBoundary
63 ) const
64 {
65  const fvMesh& mesh = meshRefiner_.mesh();
66  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
67 
68  Map<label> bafflePatch(mesh.nFaces()/1000);
69 
70  const wordList& faceZoneNames = surfaces.faceZoneNames();
71  const faceZoneMesh& fZones = mesh.faceZones();
72 
73  forAll(faceZoneNames, surfI)
74  {
75  if (faceZoneNames[surfI].size())
76  {
77  // Get zone
78  label zoneI = fZones.findZoneID(faceZoneNames[surfI]);
79 
80  const faceZone& fZone = fZones[zoneI];
81 
83  //label patchI = surfaceToCyclicPatch_[surfI];
84  // Get patch of (first region) of surface
85  label patchI = globalToPatch_[surfaces.globalRegion(surfI, 0)];
86 
87  Info<< "For surface "
88  << surfaces.names()[surfI]
89  << " found faceZone " << fZone.name()
90  << " and patch " << mesh.boundaryMesh()[patchI].name()
91  << endl;
92 
93 
94  forAll(fZone, i)
95  {
96  label faceI = fZone[i];
97 
98  if (allowBoundary || mesh.isInternalFace(faceI))
99  {
100  if (!bafflePatch.insert(faceI, patchI))
101  {
102  label oldPatchI = bafflePatch[faceI];
103 
104  if (oldPatchI != patchI)
105  {
106  FatalErrorIn("getZoneBafflePatches(const bool)")
107  << "Face " << faceI
108  << " fc:" << mesh.faceCentres()[faceI]
109  << " in zone " << fZone.name()
110  << " is in patch "
111  << mesh.boundaryMesh()[oldPatchI].name()
112  << " and in patch "
113  << mesh.boundaryMesh()[patchI].name()
114  << abort(FatalError);
115  }
116  }
117  }
118  }
119  }
120  }
121  return bafflePatch;
122 }
123 
124 
125 // Calculate geometrically collocated points, Requires PackedList to be
126 // sized and initalised!
127 Foam::label Foam::autoSnapDriver::getCollocatedPoints
128 (
129  const scalar tol,
130  const pointField& points,
131  PackedBoolList& isCollocatedPoint
132 )
133 {
134  labelList pointMap;
135  pointField newPoints;
136  bool hasMerged = mergePoints
137  (
138  points, // points
139  tol, // mergeTol
140  false, // verbose
141  pointMap,
142  newPoints
143  );
144 
145  if (!returnReduce(hasMerged, orOp<bool>()))
146  {
147  return 0;
148  }
149 
150  // Determine which newPoints are referenced more than once
151  label nCollocated = 0;
152 
153  // Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
154  // twice)
155  labelList firstOldPoint(newPoints.size(), -1);
156  forAll(pointMap, oldPointI)
157  {
158  label newPointI = pointMap[oldPointI];
159 
160  if (firstOldPoint[newPointI] == -1)
161  {
162  // First use of oldPointI. Store.
163  firstOldPoint[newPointI] = oldPointI;
164  }
165  else if (firstOldPoint[newPointI] == -2)
166  {
167  // Third or more reference of oldPointI -> non-manifold
168  isCollocatedPoint.set(oldPointI, 1u);
169  nCollocated++;
170  }
171  else
172  {
173  // Second reference of oldPointI -> non-manifold
174  isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
175  nCollocated++;
176 
177  isCollocatedPoint.set(oldPointI, 1u);
178  nCollocated++;
179 
180  // Mark with special value to save checking next time round
181  firstOldPoint[newPointI] = -2;
182  }
183  }
184  return returnReduce(nCollocated, sumOp<label>());
185 }
186 
187 
188 // Calculate displacement as average of patch points.
189 Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
190 (
191  const motionSmoother& meshMover,
192  const List<labelPair>& baffles
193 ) const
194 {
195  const indirectPrimitivePatch& pp = meshMover.patch();
196 
197  // Calculate geometrically non-manifold points on the patch to be moved.
198  PackedBoolList nonManifoldPoint(pp.nPoints());
199  label nNonManifoldPoints = getCollocatedPoints
200  (
201  SMALL,
202  pp.localPoints(),
203  nonManifoldPoint
204  );
205  Info<< "Found " << nNonManifoldPoints << " non-mainfold point(s)."
206  << endl;
207 
208 
209  // Average points
210  // ~~~~~~~~~~~~~~
211 
212  // We determine three points:
213  // - average of (centres of) connected patch faces
214  // - average of (centres of) connected internal mesh faces
215  // - as fallback: centre of any connected cell
216  // so we can do something moderately sensible for non/manifold points.
217 
218  // Note: the averages are calculated properly parallel. This is
219  // necessary to get the points shared by processors correct.
220 
221 
222  const labelListList& pointFaces = pp.pointFaces();
223  const labelList& meshPoints = pp.meshPoints();
224  const pointField& points = pp.points();
225  const polyMesh& mesh = meshMover.mesh();
226 
227  // Get labels of faces to count (master of coupled faces and baffle pairs)
228  PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
229 
230  {
231  forAll(baffles, i)
232  {
233  label f0 = baffles[i].first();
234  label f1 = baffles[i].second();
235 
236  if (isMasterFace.get(f0) == 1)
237  {
238  // Make f1 a slave
239  isMasterFace.set(f1, 0);
240  }
241  else if (isMasterFace.get(f1) == 1)
242  {
243  isMasterFace.set(f0, 0);
244  }
245  else
246  {
247  FatalErrorIn("autoSnapDriver::smoothPatchDisplacement(..)")
248  << "Both sides of baffle consisting of faces " << f0
249  << " and " << f1 << " are already slave faces."
250  << abort(FatalError);
251  }
252  }
253  }
254 
255 
256  // Get average position of boundary face centres
257  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
258 
259  vectorField avgBoundary(pointFaces.size(), vector::zero);
260  labelList nBoundary(pointFaces.size(), 0);
261 
262  forAll(pointFaces, patchPointI)
263  {
264  const labelList& pFaces = pointFaces[patchPointI];
265 
266  forAll(pFaces, pfI)
267  {
268  label faceI = pFaces[pfI];
269 
270  if (isMasterFace.get(pp.addressing()[faceI]) == 1)
271  {
272  avgBoundary[patchPointI] += pp[faceI].centre(points);
273  nBoundary[patchPointI]++;
274  }
275  }
276  }
277 
279  (
280  mesh,
281  pp.meshPoints(),
282  avgBoundary,
283  plusEqOp<point>(), // combine op
284  vector::zero, // null value
285  false // no separation
286  );
288  (
289  mesh,
290  pp.meshPoints(),
291  nBoundary,
292  plusEqOp<label>(), // combine op
293  0, // null value
294  false // no separation
295  );
296 
297  forAll(avgBoundary, i)
298  {
299  avgBoundary[i] /= nBoundary[i];
300  }
301 
302 
303  // Get average position of internal face centres
304  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
305 
306  vectorField avgInternal;
307  labelList nInternal;
308  {
309  vectorField globalSum(mesh.nPoints(), vector::zero);
310  labelList globalNum(mesh.nPoints(), 0);
311 
312  // Note: no use of pointFaces
313  const faceList& faces = mesh.faces();
314 
315  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
316  {
317  const face& f = faces[faceI];
318  const point& fc = mesh.faceCentres()[faceI];
319 
320  forAll(f, fp)
321  {
322  globalSum[f[fp]] += fc;
323  globalNum[f[fp]]++;
324  }
325  }
326 
327  // Count coupled faces as internal ones (but only once)
328  const polyBoundaryMesh& patches = mesh.boundaryMesh();
329 
330  forAll(patches, patchI)
331  {
332  if (Pstream::parRun() && isA<processorPolyPatch>(patches[patchI]))
333  {
334  const processorPolyPatch& pp =
335  refCast<const processorPolyPatch>(patches[patchI]);
336 
337  if (pp.myProcNo() < pp.neighbProcNo())
338  {
339  const vectorField::subField faceCentres = pp.faceCentres();
340 
341  forAll(pp, i)
342  {
343  const face& f = pp[i];
344  const point& fc = faceCentres[i];
345 
346  forAll(f, fp)
347  {
348  globalSum[f[fp]] += fc;
349  globalNum[f[fp]]++;
350  }
351  }
352  }
353  }
354  else if (isA<cyclicPolyPatch>(patches[patchI]))
355  {
356  const cyclicPolyPatch& pp =
357  refCast<const cyclicPolyPatch>(patches[patchI]);
358 
359  const vectorField::subField faceCentres = pp.faceCentres();
360 
361  for (label i = 0; i < pp.size()/2; i++)
362  {
363  const face& f = pp[i];
364  const point& fc = faceCentres[i];
365 
366  forAll(f, fp)
367  {
368  globalSum[f[fp]] += fc;
369  globalNum[f[fp]]++;
370  }
371  }
372  }
373  }
374 
376  (
377  mesh,
378  globalSum,
379  plusEqOp<vector>(), // combine op
380  vector::zero, // null value
381  false // no separation
382  );
384  (
385  mesh,
386  globalNum,
387  plusEqOp<label>(), // combine op
388  0, // null value
389  false // no separation
390  );
391 
392  avgInternal.setSize(meshPoints.size());
393  nInternal.setSize(meshPoints.size());
394 
395  forAll(avgInternal, patchPointI)
396  {
397  label meshPointI = meshPoints[patchPointI];
398 
399  nInternal[patchPointI] = globalNum[meshPointI];
400 
401  if (nInternal[patchPointI] == 0)
402  {
403  avgInternal[patchPointI] = globalSum[meshPointI];
404  }
405  else
406  {
407  avgInternal[patchPointI] =
408  globalSum[meshPointI]
409  / nInternal[patchPointI];
410  }
411  }
412  }
413 
414 
415  // Precalculate any cell using mesh point (replacement of pointCells()[])
416  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
417 
418  labelList anyCell(mesh.nPoints(), -1);
419  forAll(mesh.faceNeighbour(), faceI)
420  {
421  label own = mesh.faceOwner()[faceI];
422  const face& f = mesh.faces()[faceI];
423 
424  forAll(f, fp)
425  {
426  anyCell[f[fp]] = own;
427  }
428  }
429  for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
430  {
431  label own = mesh.faceOwner()[faceI];
432 
433  const face& f = mesh.faces()[faceI];
434 
435  forAll(f, fp)
436  {
437  anyCell[f[fp]] = own;
438  }
439  }
440 
441 
442  // Displacement to calculate.
443  pointField patchDisp(meshPoints.size(), vector::zero);
444 
445  forAll(pointFaces, i)
446  {
447  label meshPointI = meshPoints[i];
448  const point& currentPos = pp.points()[meshPointI];
449 
450  // Now we have the two average points: avgBoundary and avgInternal
451  // and how many boundary/internal faces connect to the point
452  // (nBoundary, nInternal)
453  // Do some blending between the two.
454  // Note: the following section has some reasoning behind it but the
455  // blending factors can be experimented with.
456 
457  point newPos;
458 
459  if (nonManifoldPoint.get(i) == 0u)
460  {
461  // Points that are manifold. Weight the internal and boundary
462  // by their number of faces and blend with
463  scalar internalBlend = 0.1;
464  scalar blend = 0.1;
465 
466  point avgPos =
467  (
468  internalBlend*nInternal[i]*avgInternal[i]
469  +(1-internalBlend)*nBoundary[i]*avgBoundary[i]
470  )
471  / (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
472 
473  newPos = (1-blend)*avgPos + blend*currentPos;
474  }
475  else if (nInternal[i] == 0)
476  {
477  // Non-manifold without internal faces. Use any connected cell
478  // as internal point instead. Use precalculated any cell to avoid
479  // e.g. pointCells()[meshPointI][0]
480 
481  const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
482 
483  scalar cellCBlend = 0.8;
484  scalar blend = 0.1;
485 
486  point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
487 
488  newPos = (1-blend)*avgPos + blend*currentPos;
489  }
490  else
491  {
492  // Non-manifold point with internal faces connected to them
493  scalar internalBlend = 0.9;
494  scalar blend = 0.1;
495 
496  point avgPos =
497  internalBlend*avgInternal[i]
498  + (1-internalBlend)*avgBoundary[i];
499 
500  newPos = (1-blend)*avgPos + blend*currentPos;
501  }
502 
503  patchDisp[i] = newPos - currentPos;
504  }
505 
506  return patchDisp;
507 }
508 
509 
510 Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::edgePatchDist
511 (
512  const pointMesh& pMesh,
513  const indirectPrimitivePatch& pp
514 )
515 {
516  const polyMesh& mesh = pMesh();
517 
518  // Set initial changed points to all the patch points
519  List<pointEdgePoint> wallInfo(pp.nPoints());
520 
521  forAll(pp.localPoints(), ppI)
522  {
523  wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
524  }
525 
526  // Current info on points
527  List<pointEdgePoint> allPointInfo(mesh.nPoints());
528 
529  // Current info on edges
530  List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
531 
532  PointEdgeWave<pointEdgePoint> wallCalc
533  (
534  mesh,
535  pp.meshPoints(),
536  wallInfo,
537 
538  allPointInfo,
539  allEdgeInfo,
540  mesh.globalData().nTotalPoints() // max iterations
541  );
542 
543  // Copy edge values into scalarField
544  tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
545  scalarField& edgeDist = tedgeDist();
546 
547  forAll(allEdgeInfo, edgeI)
548  {
549  edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
550  }
551 
552 
553  //{
554  // // For debugging: dump to file
555  // pointScalarField pointDist
556  // (
557  // IOobject
558  // (
559  // "pointDist",
560  // meshRefiner_.timeName(),
561  // mesh.DB(),
562  // IOobject::NO_READ,
563  // IOobject::AUTO_WRITE
564  // ),
565  // pMesh,
566  // dimensionedScalar("pointDist", dimless, 0.0)
567  // );
568  //
569  // forAll(allEdgeInfo, edgeI)
570  // {
571  // scalar d = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
572  //
573  // const edge& e = mesh.edges()[edgeI];
574  //
575  // pointDist[e[0]] += d;
576  // pointDist[e[1]] += d;
577  // }
578  // forAll(pointDist, pointI)
579  // {
580  // pointDist[pointI] /= mesh.pointEdges()[pointI].size();
581  // }
582  // Info<< "Writing patch distance to " << pointDist.name()
583  // << " at time " << meshRefiner_.timeName() << endl;
584  //
585  // pointDist.write();
586  //}
587 
588  return tedgeDist;
589 }
590 
591 
592 void Foam::autoSnapDriver::dumpMove
593 (
594  const fileName& fName,
595  const pointField& meshPts,
596  const pointField& surfPts
597 )
598 {
599  // Dump direction of growth into file
600  Pout<< nl << "Dumping move direction to " << fName << nl
601  << "View this Lightwave-OBJ file with e.g. javaview" << nl
602  << endl;
603 
604  OFstream nearestStream(fName);
605 
606  label vertI = 0;
607 
608  forAll(meshPts, ptI)
609  {
610  meshTools::writeOBJ(nearestStream, meshPts[ptI]);
611  vertI++;
612 
613  meshTools::writeOBJ(nearestStream, surfPts[ptI]);
614  vertI++;
615 
616  nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
617  }
618 }
619 
620 
621 // Check whether all displacement vectors point outwards of patch. Return true
622 // if so.
623 bool Foam::autoSnapDriver::outwardsDisplacement
624 (
625  const indirectPrimitivePatch& pp,
626  const vectorField& patchDisp
627 )
628 {
629  const vectorField& faceNormals = pp.faceNormals();
630  const labelListList& pointFaces = pp.pointFaces();
631 
632  forAll(pointFaces, pointI)
633  {
634  const labelList& pFaces = pointFaces[pointI];
635 
636  vector disp(patchDisp[pointI]);
637 
638  scalar magDisp = mag(disp);
639 
640  if (magDisp > SMALL)
641  {
642  disp /= magDisp;
643 
644  bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
645 
646  if (!outwards)
647  {
648  Warning<< "Displacement " << patchDisp[pointI]
649  << " at mesh point " << pp.meshPoints()[pointI]
650  << " coord " << pp.points()[pp.meshPoints()[pointI]]
651  << " points through the surrounding patch faces" << endl;
652  return false;
653  }
654  }
655  else
656  {
657  //? Displacement small but in wrong direction. Would probably be ok.
658  }
659  }
660  return true;
661 }
662 
663 
664 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
665 
666 Foam::autoSnapDriver::autoSnapDriver
667 (
668  meshRefinement& meshRefiner,
669  const labelList& globalToPatch
670 )
671 :
672  meshRefiner_(meshRefiner),
673  globalToPatch_(globalToPatch)
674 {}
675 
676 
677 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
678 
680 (
681  List<labelPair>& baffles
682 )
683 {
684  labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
685 
687 
688  // No need to sync; all processors will have all same zonedSurfaces.
689  if (zonedSurfaces.size())
690  {
691  fvMesh& mesh = meshRefiner_.mesh();
692 
693  // Split internal faces on interface surfaces
694  Info<< "Converting zoned faces into baffles ..." << endl;
695 
696  // Get faces (internal only) to be baffled. Map from face to patch
697  // label.
698  Map<label> faceToPatch(getZoneBafflePatches(false));
699 
700  label nZoneFaces = returnReduce(faceToPatch.size(), sumOp<label>());
701  if (nZoneFaces > 0)
702  {
703  // Convert into labelLists
704  labelList ownPatch(mesh.nFaces(), -1);
705  forAllConstIter(Map<label>, faceToPatch, iter)
706  {
707  ownPatch[iter.key()] = iter();
708  }
709 
710  // Create baffles. both sides same patch.
711  map = meshRefiner_.createBaffles(ownPatch, ownPatch);
712 
713  // Get pairs of faces created.
714  // Just loop over faceMap and store baffle if we encounter a slave
715  // face.
716 
717  baffles.setSize(faceToPatch.size());
718  label baffleI = 0;
719 
720  const labelList& faceMap = map().faceMap();
721  const labelList& reverseFaceMap = map().reverseFaceMap();
722 
723  forAll(faceMap, faceI)
724  {
725  label oldFaceI = faceMap[faceI];
726 
727  // Does face originate from face-to-patch
728  Map<label>::const_iterator iter = faceToPatch.find(oldFaceI);
729 
730  if (iter != faceToPatch.end())
731  {
732  label masterFaceI = reverseFaceMap[oldFaceI];
733  if (faceI != masterFaceI)
734  {
735  baffles[baffleI++] = labelPair(masterFaceI, faceI);
736  }
737  }
738  }
739 
740  if (baffleI != faceToPatch.size())
741  {
742  FatalErrorIn("autoSnapDriver::createZoneBaffles(..)")
743  << "Had " << faceToPatch.size() << " patches to create "
744  << " but encountered " << baffleI
745  << " slave faces originating from patcheable faces."
746  << abort(FatalError);
747  }
748 
749  if (debug)
750  {
751  const_cast<Time&>(mesh.time())++;
752  Pout<< "Writing baffled mesh to time "
753  << meshRefiner_.timeName() << endl;
754  mesh.write();
755  }
756  }
757  Info<< "Created " << nZoneFaces << " baffles in = "
758  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
759  }
760  return map;
761 }
762 
763 
765 (
766  const List<labelPair>& baffles
767 )
768 {
769  labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
770 
772 
773  // No need to sync; all processors will have all same zonedSurfaces.
774  label nBaffles = returnReduce(baffles.size(), sumOp<label>());
775  if (zonedSurfaces.size() && nBaffles > 0)
776  {
777  // Merge any baffles
778  Info<< "Converting " << nBaffles << " baffles back into zoned faces ..."
779  << endl;
780 
781  map = meshRefiner_.mergeBaffles(baffles);
782 
783  Info<< "Converted baffles in = "
784  << meshRefiner_.mesh().time().cpuTimeIncrement()
785  << " s\n" << nl << endl;
786  }
787 
788  return map;
789 }
790 
791 
793 (
794  const snapParameters& snapParams,
795  const indirectPrimitivePatch& pp
796 ) const
797 {
798  const edgeList& edges = pp.edges();
799  const labelListList& pointEdges = pp.pointEdges();
800  const pointField& localPoints = pp.localPoints();
801  const fvMesh& mesh = meshRefiner_.mesh();
802 
803  scalarField maxEdgeLen(localPoints.size(), -GREAT);
804 
805  forAll(pointEdges, pointI)
806  {
807  const labelList& pEdges = pointEdges[pointI];
808 
809  forAll(pEdges, pEdgeI)
810  {
811  const edge& e = edges[pEdges[pEdgeI]];
812 
813  scalar len = e.mag(localPoints);
814 
815  maxEdgeLen[pointI] = max(maxEdgeLen[pointI], len);
816  }
817  }
818 
820  (
821  mesh,
822  pp.meshPoints(),
823  maxEdgeLen,
824  maxEqOp<scalar>(), // combine op
825  -GREAT, // null value
826  false // no separation
827  );
828 
829  return snapParams.snapTol()*maxEdgeLen;
830 }
831 
832 
834 (
835  const snapParameters& snapParams,
836  const label nInitErrors,
837  const List<labelPair>& baffles,
838  motionSmoother& meshMover
839 ) const
840 {
841  const fvMesh& mesh = meshRefiner_.mesh();
842 
843  labelList checkFaces;
844 
845  Info<< "Smoothing patch points ..." << endl;
846  for
847  (
848  label smoothIter = 0;
849  smoothIter < snapParams.nSmoothPatch();
850  smoothIter++
851  )
852  {
853  Info<< "Smoothing iteration " << smoothIter << endl;
854  checkFaces.setSize(mesh.nFaces());
855  forAll(checkFaces, faceI)
856  {
857  checkFaces[faceI] = faceI;
858  }
859 
860  pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
861 
862  // The current mesh is the starting mesh to smooth from.
863  meshMover.setDisplacement(patchDisp);
864  meshMover.correct();
865 
866  scalar oldErrorReduction = -1;
867 
868  for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
869  {
870  Info<< nl << "Scaling iteration " << snapIter << endl;
871 
872  if (snapIter == snapParams.nSnap())
873  {
874  Info<< "Displacement scaling for error reduction set to 0."
875  << endl;
876  oldErrorReduction = meshMover.setErrorReduction(0.0);
877  }
878 
879  // Try to adapt mesh to obtain displacement by smoothly
880  // decreasing displacement at error locations.
881  if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
882  {
883  Info<< "Successfully moved mesh" << endl;
884  break;
885  }
886  }
887 
888  if (oldErrorReduction >= 0)
889  {
890  meshMover.setErrorReduction(oldErrorReduction);
891  }
892  Info<< endl;
893  }
894 
895 
896  // The current mesh is the starting mesh to smooth from.
897  meshMover.correct();
898 
899  if (debug)
900  {
901  const_cast<Time&>(mesh.time())++;
902  Pout<< "Writing patch smoothed mesh to time " << meshRefiner_.timeName()
903  << endl;
904 
905  mesh.write();
906  }
907 
908  Info<< "Patch points smoothed in = "
909  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
910 }
911 
912 
913 // Get (pp-local) indices of points that are both on zone and on patched surface
915 (
916  const indirectPrimitivePatch& pp,
917  const word& zoneName
918 ) const
919 {
920  const fvMesh& mesh = meshRefiner_.mesh();
921 
922  label zoneI = mesh.faceZones().findZoneID(zoneName);
923 
924  if (zoneI == -1)
925  {
927  (
928  "autoSnapDriver::getZoneSurfacePoints"
929  "(const indirectPrimitivePatch&, const word&)"
930  ) << "Cannot find zone " << zoneName
931  << exit(FatalError);
932  }
933 
934  const faceZone& fZone = mesh.faceZones()[zoneI];
935 
936 
937  // Could use PrimitivePatch & localFaces to extract points but might just
938  // as well do it ourselves.
939 
940  boolList pointOnZone(pp.nPoints(), false);
941 
942  forAll(fZone, i)
943  {
944  const face& f = mesh.faces()[fZone[i]];
945 
946  forAll(f, fp)
947  {
948  label meshPointI = f[fp];
949 
951  pp.meshPointMap().find(meshPointI);
952 
953  if (iter != pp.meshPointMap().end())
954  {
955  label pointI = iter();
956  pointOnZone[pointI] = true;
957  }
958  }
959  }
960 
961  return findIndices(pointOnZone, true);
962 }
963 
964 
966 (
967  const scalarField& snapDist,
968  motionSmoother& meshMover
969 ) const
970 {
971  Info<< "Calculating patchDisplacement as distance to nearest surface"
972  << " point ..." << endl;
973 
974  const indirectPrimitivePatch& pp = meshMover.patch();
975  const pointField& localPoints = pp.localPoints();
976  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
977  const fvMesh& mesh = meshRefiner_.mesh();
978 
979  // Displacement per patch point
980  vectorField patchDisp(localPoints.size(), vector::zero);
981 
982  if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
983  {
984  // Current surface snapped to
985  labelList snapSurf(localPoints.size(), -1);
986 
987  // Divide surfaces into zoned and unzoned
988  labelList zonedSurfaces =
989  meshRefiner_.surfaces().getNamedSurfaces();
990  labelList unzonedSurfaces =
991  meshRefiner_.surfaces().getUnnamedSurfaces();
992 
993 
994  // 1. All points to non-interface surfaces
995  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
996 
997  {
998  List<pointIndexHit> hitInfo;
999  labelList hitSurface;
1000  surfaces.findNearest
1001  (
1002  unzonedSurfaces,
1003  localPoints,
1004  sqr(4*snapDist), // sqr of attract distance
1005  hitSurface,
1006  hitInfo
1007  );
1008 
1009  forAll(hitInfo, pointI)
1010  {
1011  if (hitInfo[pointI].hit())
1012  {
1013  patchDisp[pointI] =
1014  hitInfo[pointI].hitPoint()
1015  - localPoints[pointI];
1016 
1017  snapSurf[pointI] = hitSurface[pointI];
1018  }
1019  }
1020  }
1021 
1022 
1023 
1024  // 2. All points on zones to their respective surface
1025  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1026 
1027  // Surfaces with zone information
1028  const wordList& faceZoneNames = surfaces.faceZoneNames();
1029 
1030  // Current best snap distance
1031  scalarField minSnapDist(snapDist);
1032 
1033  forAll(zonedSurfaces, i)
1034  {
1035  label zoneSurfI = zonedSurfaces[i];
1036 
1037  const labelList surfacesToTest(1, zoneSurfI);
1038 
1039  // Get indices of points both on faceZone and on pp.
1040  labelList zonePointIndices
1041  (
1042  getZoneSurfacePoints
1043  (
1044  pp,
1045  faceZoneNames[zoneSurfI]
1046  )
1047  );
1048 
1049  // Find nearest for points both on faceZone and pp.
1050  List<pointIndexHit> hitInfo;
1051  labelList hitSurface;
1052  surfaces.findNearest
1053  (
1054  labelList(1, zoneSurfI),
1055  pointField(localPoints, zonePointIndices),
1056  sqr(4*scalarField(minSnapDist, zonePointIndices)),
1057  hitSurface,
1058  hitInfo
1059  );
1060 
1061  forAll(hitInfo, i)
1062  {
1063  label pointI = zonePointIndices[i];
1064 
1065  if (hitInfo[i].hit())
1066  {
1067  patchDisp[pointI] =
1068  hitInfo[i].hitPoint()
1069  - localPoints[pointI];
1070 
1071  minSnapDist[pointI] = min
1072  (
1073  minSnapDist[pointI],
1074  mag(patchDisp[pointI])
1075  );
1076 
1077  snapSurf[pointI] = zoneSurfI;
1078  }
1079  }
1080  }
1081 
1082  // Check if all points are being snapped
1083  forAll(snapSurf, pointI)
1084  {
1085  if (snapSurf[pointI] == -1)
1086  {
1087  WarningIn("autoSnapDriver::calcNearestSurface(..)")
1088  << "For point:" << pointI
1089  << " coordinate:" << localPoints[pointI]
1090  << " did not find any surface within:"
1091  << minSnapDist[pointI]
1092  << " meter." << endl;
1093  }
1094  }
1095 
1096  {
1097  scalarField magDisp(mag(patchDisp));
1098 
1099  Info<< "Wanted displacement : average:"
1100  << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
1101  << " min:" << gMin(magDisp)
1102  << " max:" << gMax(magDisp) << endl;
1103  }
1104  }
1105 
1106  Info<< "Calculated surface displacement in = "
1107  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1108 
1109 
1110  // Limit amount of movement.
1111  forAll(patchDisp, patchPointI)
1112  {
1113  scalar magDisp = mag(patchDisp[patchPointI]);
1114 
1115  if (magDisp > snapDist[patchPointI])
1116  {
1117  patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
1118 
1119  Pout<< "Limiting displacement for " << patchPointI
1120  << " from " << magDisp << " to " << snapDist[patchPointI]
1121  << endl;
1122  }
1123  }
1124 
1125  // Points on zones in one domain but only present as point on other
1126  // will not do condition 2 on all. Sync explicitly.
1128  (
1129  mesh,
1130  pp.meshPoints(),
1131  patchDisp,
1132  minMagEqOp(), // combine op
1133  vector(GREAT, GREAT, GREAT), // null value
1134  false // no separation
1135  );
1136 
1137 
1138  // Check for displacement being outwards.
1139  outwardsDisplacement(pp, patchDisp);
1140 
1141  // Set initial distribution of displacement field (on patches) from
1142  // patchDisp and make displacement consistent with b.c. on displacement
1143  // pointVectorField.
1144  meshMover.setDisplacement(patchDisp);
1145 
1146  if (debug)
1147  {
1148  dumpMove
1149  (
1150  mesh.time().path()/"patchDisplacement.obj",
1151  pp.localPoints(),
1152  pp.localPoints() + patchDisp
1153  );
1154  }
1155 
1156  return patchDisp;
1157 }
1158 
1159 
1162  const snapParameters& snapParams,
1163  motionSmoother& meshMover
1164 ) const
1165 {
1166  const fvMesh& mesh = meshRefiner_.mesh();
1167  const pointMesh& pMesh = meshMover.pMesh();
1168  const indirectPrimitivePatch& pp = meshMover.patch();
1169 
1170  Info<< "Smoothing displacement ..." << endl;
1171 
1172  // Set edge diffusivity as inverse of distance to patch
1173  scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL));
1174  //scalarField edgeGamma(mesh.nEdges(), 1.0);
1175  //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
1176 
1177  // Get displacement field
1178  pointVectorField& disp = meshMover.displacement();
1179 
1180  for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
1181  {
1182  if ((iter % 10) == 0)
1183  {
1184  Info<< "Iteration " << iter << endl;
1185  }
1186  pointVectorField oldDisp(disp);
1187 
1188  meshMover.smooth(oldDisp, edgeGamma, false, disp);
1189  }
1190  Info<< "Displacement smoothed in = "
1191  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1192 
1193  if (debug)
1194  {
1195  const_cast<Time&>(mesh.time())++;
1196  Pout<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
1197  << endl;
1198 
1199  // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
1200  // but this will also delete all pointMesh but not pointFields which
1201  // gives an illegal situation.
1202 
1203  mesh.write();
1204 
1205  Pout<< "Writing displacement field ..." << endl;
1206  disp.write();
1207  tmp<pointScalarField> magDisp(mag(disp));
1208  magDisp().write();
1209 
1210  Pout<< "Writing actual patch displacement ..." << endl;
1211  vectorField actualPatchDisp(disp, pp.meshPoints());
1212  dumpMove
1213  (
1214  mesh.time().path()/"actualPatchDisplacement.obj",
1215  pp.localPoints(),
1216  pp.localPoints() + actualPatchDisp
1217  );
1218  }
1219 }
1220 
1221 
1224  const snapParameters& snapParams,
1225  const label nInitErrors,
1226  const List<labelPair>& baffles,
1227  motionSmoother& meshMover
1228 )
1229 {
1230  const fvMesh& mesh = meshRefiner_.mesh();
1231 
1232  // Relax displacement until correct mesh
1233  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1234  labelList checkFaces(identity(mesh.nFaces()));
1235 
1236  scalar oldErrorReduction = -1;
1237 
1238  Info<< "Moving mesh ..." << endl;
1239  for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
1240  {
1241  Info<< nl << "Iteration " << iter << endl;
1242 
1243  if (iter == snapParams.nSnap())
1244  {
1245  Info<< "Displacement scaling for error reduction set to 0." << endl;
1246  oldErrorReduction = meshMover.setErrorReduction(0.0);
1247  }
1248 
1249  if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
1250  {
1251  Info<< "Successfully moved mesh" << endl;
1252 
1253  break;
1254  }
1255  if (debug)
1256  {
1257  const_cast<Time&>(mesh.time())++;
1258  Pout<< "Writing scaled mesh to time " << meshRefiner_.timeName()
1259  << endl;
1260  mesh.write();
1261 
1262  Pout<< "Writing displacement field ..." << endl;
1263  meshMover.displacement().write();
1264  tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
1265  magDisp().write();
1266  }
1267  }
1268 
1269  if (oldErrorReduction >= 0)
1270  {
1271  meshMover.setErrorReduction(oldErrorReduction);
1272  }
1273  Info<< "Moved mesh in = "
1274  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1275 }
1276 
1277 
1278 // After snapping: correct patching according to nearest surface.
1279 // Code is very similar to calcNearestSurface.
1280 // - calculate face-wise snap distance as max of point-wise
1281 // - calculate face-wise nearest surface point
1282 // - repatch face according to patch for surface point.
1285  const snapParameters& snapParams,
1286  const labelList& adaptPatchIDs
1287 )
1288 {
1289  const fvMesh& mesh = meshRefiner_.mesh();
1290  const refinementSurfaces& surfaces = meshRefiner_.surfaces();
1291 
1292  Info<< "Repatching faces according to nearest surface ..." << endl;
1293 
1294  // Get the labels of added patches.
1296  (
1298  (
1299  mesh,
1300  adaptPatchIDs
1301  )
1302  );
1303  indirectPrimitivePatch& pp = ppPtr();
1304 
1305  // Divide surfaces into zoned and unzoned
1306  labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
1307  labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
1308 
1309 
1310  // Faces that do not move
1311  PackedBoolList isZonedFace(mesh.nFaces(), 0);
1312  {
1313  // 1. All faces on zoned surfaces
1314  const wordList& faceZoneNames = surfaces.faceZoneNames();
1315  const faceZoneMesh& fZones = mesh.faceZones();
1316 
1317  forAll(zonedSurfaces, i)
1318  {
1319  label zoneSurfI = zonedSurfaces[i];
1320 
1321  label zoneI = fZones.findZoneID(faceZoneNames[zoneSurfI]);
1322 
1323  const faceZone& fZone = fZones[zoneI];
1324 
1325  forAll(fZone, i)
1326  {
1327  isZonedFace.set(fZone[i], 1);
1328  }
1329  }
1330  }
1331 
1332 
1333  // Determine per pp face which patch it should be in
1334  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1335 
1336  // Patch that face should be in
1337  labelList closestPatch(pp.size(), -1);
1338  {
1339  // face snap distance as max of point snap distance
1340  scalarField faceSnapDist(pp.size(), -GREAT);
1341  {
1342  // Distance to attract to nearest feature on surface
1343  const scalarField snapDist(calcSnapDistance(snapParams, pp));
1344 
1345  const faceList& localFaces = pp.localFaces();
1346 
1347  forAll(localFaces, faceI)
1348  {
1349  const face& f = localFaces[faceI];
1350 
1351  forAll(f, fp)
1352  {
1353  faceSnapDist[faceI] = max
1354  (
1355  faceSnapDist[faceI],
1356  snapDist[f[fp]]
1357  );
1358  }
1359  }
1360  }
1361 
1362  pointField localFaceCentres(mesh.faceCentres(), pp.addressing());
1363 
1364  // Get nearest surface and region
1365  labelList hitSurface;
1366  labelList hitRegion;
1367  surfaces.findNearestRegion
1368  (
1369  unzonedSurfaces,
1370  localFaceCentres,
1371  sqr(4*faceSnapDist), // sqr of attract distance
1372  hitSurface,
1373  hitRegion
1374  );
1375 
1376  // Get patch
1377  forAll(pp, i)
1378  {
1379  label faceI = pp.addressing()[i];
1380 
1381  if (hitSurface[i] != -1 && (isZonedFace.get(faceI) == 0))
1382  {
1383  closestPatch[i] = globalToPatch_
1384  [
1385  surfaces.globalRegion
1386  (
1387  hitSurface[i],
1388  hitRegion[i]
1389  )
1390  ];
1391  }
1392  }
1393  }
1394 
1395 
1396  // Change those faces for which there is a different closest patch
1397  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1398 
1399  labelList ownPatch(mesh.nFaces(), -1);
1400  labelList neiPatch(mesh.nFaces(), -1);
1401 
1402  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1403 
1404  forAll(patches, patchI)
1405  {
1406  const polyPatch& pp = patches[patchI];
1407 
1408  forAll(pp, i)
1409  {
1410  ownPatch[pp.start()+i] = patchI;
1411  neiPatch[pp.start()+i] = patchI;
1412  }
1413  }
1414 
1415  label nChanged = 0;
1416  forAll(closestPatch, i)
1417  {
1418  label faceI = pp.addressing()[i];
1419 
1420  if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
1421  {
1422  ownPatch[faceI] = closestPatch[i];
1423  neiPatch[faceI] = closestPatch[i];
1424  nChanged++;
1425  }
1426  }
1427 
1428  Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
1429  << " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
1430  << endl;
1431 
1432  return meshRefiner_.createBaffles(ownPatch, neiPatch);
1433 }
1434 
1435 
1438  const dictionary& snapDict,
1439  const dictionary& motionDict,
1440  const snapParameters& snapParams
1441 )
1442 {
1443  fvMesh& mesh = meshRefiner_.mesh();
1444 
1445  Info<< nl
1446  << "Morphing phase" << nl
1447  << "--------------" << nl
1448  << endl;
1449 
1450  // Get the labels of added patches.
1451  labelList adaptPatchIDs(meshRefiner_.meshedPatches());
1452 
1453  // Create baffles (pairs of faces that share the same points)
1454  // Baffles stored as owner and neighbour face that have been created.
1455  List<labelPair> baffles;
1456  createZoneBaffles(baffles);
1457 
1458  {
1460  (
1462  (
1463  mesh,
1464  adaptPatchIDs
1465  )
1466  );
1467  indirectPrimitivePatch& pp = ppPtr();
1468 
1469  // Distance to attract to nearest feature on surface
1470  const scalarField snapDist(calcSnapDistance(snapParams, pp));
1471 
1472 
1473  // Construct iterative mesh mover.
1474  Info<< "Constructing mesh displacer ..." << endl;
1475  Info<< "Using mesh parameters " << motionDict << nl << endl;
1476 
1477  const pointMesh& pMesh = pointMesh::New(mesh);
1478 
1479  motionSmoother meshMover
1480  (
1481  mesh,
1482  pp,
1483  adaptPatchIDs,
1484  meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
1485  motionDict
1486  );
1487 
1488 
1489  // Check initial mesh
1490  Info<< "Checking initial mesh ..." << endl;
1491  labelHashSet wrongFaces(mesh.nFaces()/100);
1492  motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
1493  const label nInitErrors = returnReduce
1494  (
1495  wrongFaces.size(),
1496  sumOp<label>()
1497  );
1498 
1499  Info<< "Detected " << nInitErrors << " illegal faces"
1500  << " (concave, zero area or negative cell pyramid volume)"
1501  << endl;
1502 
1503 
1504  Info<< "Checked initial mesh in = "
1505  << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
1506 
1507  // Pre-smooth patch vertices (so before determining nearest)
1508  preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
1509 
1510  // Calculate displacement at every patch point. Insert into
1511  // meshMover.
1512  calcNearestSurface(snapDist, meshMover);
1513 
1515  //- 2009-12-16 : was not found to be beneficial. Keeping internal
1516  // fields fixed slightly increases skewness (on boundary)
1517  // but lowers non-orthogonality quite a bit (e.g. 65->59 degrees).
1518  // Maybe if better smoother?
1519  //smoothDisplacement(snapParams, meshMover);
1520 
1521  // Apply internal displacement to mesh.
1522  scaleMesh(snapParams, nInitErrors, baffles, meshMover);
1523  }
1524 
1525  // Merge any introduced baffles.
1526  mergeZoneBaffles(baffles);
1527 
1528  // Repatch faces according to nearest.
1529  repatchToSurface(snapParams, adaptPatchIDs);
1530 }
1531 
1532 
1533 // ************************ vim: set sw=4 sts=4 et: ************************ //