Newer
Older
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
All to do with snapping to the surface
\*----------------------------------------------------------------------------*/
#include "autoSnapDriver.H"
#include "motionSmoother.H"
#include "polyTopoChange.H"
#include "OFstream.H"
#include "syncTools.H"
#include "fvMesh.H"
#include "Time.H"
#include "OFstream.H"
#include "mapPolyMesh.H"
#include "pointEdgePoint.H"
#include "PointEdgeWave.H"
#include "mergePoints.H"
#include "snapParameters.H"
#include "refinementSurfaces.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
defineTypeNameAndDebug(autoSnapDriver, 0);
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Calculate geometrically collocated points, Requires PackedList to be
// sized and initalised!
Foam::label Foam::autoSnapDriver::getCollocatedPoints
(
const scalar tol,
const pointField& points,
PackedBoolList& isCollocatedPoint
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
)
{
labelList pointMap;
pointField newPoints;
bool hasMerged = mergePoints
(
points, // points
tol, // mergeTol
false, // verbose
pointMap,
newPoints
);
if (!returnReduce(hasMerged, orOp<bool>()))
{
return 0;
}
// Determine which newPoints are referenced more than once
label nCollocated = 0;
// Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
// twice)
labelList firstOldPoint(newPoints.size(), -1);
forAll(pointMap, oldPointI)
{
label newPointI = pointMap[oldPointI];
if (firstOldPoint[newPointI] == -1)
{
// First use of oldPointI. Store.
firstOldPoint[newPointI] = oldPointI;
}
else if (firstOldPoint[newPointI] == -2)
{
// Third or more reference of oldPointI -> non-manifold
isCollocatedPoint.set(oldPointI, 1u);
nCollocated++;
}
else
{
// Second reference of oldPointI -> non-manifold
isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
nCollocated++;
isCollocatedPoint.set(oldPointI, 1u);
nCollocated++;
// Mark with special value to save checking next time round
firstOldPoint[newPointI] = -2;
}
}
return returnReduce(nCollocated, sumOp<label>());
}
// Calculate displacement as average of patch points.
Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
(
const motionSmoother& meshMover,
const List<labelPair>& baffles
) const
{
const indirectPrimitivePatch& pp = meshMover.patch();
// Calculate geometrically non-manifold points on the patch to be moved.
PackedBoolList nonManifoldPoint(pp.nPoints());
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
label nNonManifoldPoints = getCollocatedPoints
(
SMALL,
pp.localPoints(),
nonManifoldPoint
);
Info<< "Found " << nNonManifoldPoints << " non-mainfold point(s)."
<< endl;
// Average points
// ~~~~~~~~~~~~~~
// We determine three points:
// - average of (centres of) connected patch faces
// - average of (centres of) connected internal mesh faces
// - as fallback: centre of any connected cell
// so we can do something moderately sensible for non/manifold points.
// Note: the averages are calculated properly parallel. This is
// necessary to get the points shared by processors correct.
const labelListList& pointFaces = pp.pointFaces();
const labelList& meshPoints = pp.meshPoints();
const pointField& points = pp.points();
const polyMesh& mesh = meshMover.mesh();
// Get labels of faces to count (master of coupled faces and baffle pairs)
PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
{
forAll(baffles, i)
{
label f0 = baffles[i].first();
label f1 = baffles[i].second();
if (isMasterFace.get(f0))
else if (isMasterFace.get(f1))
}
else
{
FatalErrorIn("autoSnapDriver::smoothPatchDisplacement(..)")
<< "Both sides of baffle consisting of faces " << f0
<< " and " << f1 << " are already slave faces."
<< abort(FatalError);
}
}
}
// Get average position of boundary face centres
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vectorField avgBoundary(pointFaces.size(), vector::zero);
labelList nBoundary(pointFaces.size(), 0);
forAll(pointFaces, patchPointI)
{
const labelList& pFaces = pointFaces[patchPointI];
forAll(pFaces, pfI)
{
if (isMasterFace.get(pp.addressing()[faceI]))
{
avgBoundary[patchPointI] += pp[faceI].centre(points);
nBoundary[patchPointI]++;
}
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
avgBoundary,
plusEqOp<point>(), // combine op
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
nBoundary,
plusEqOp<label>(), // combine op
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
);
forAll(avgBoundary, i)
{
avgBoundary[i] /= nBoundary[i];
}
// Get average position of internal face centres
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vectorField avgInternal;
labelList nInternal;
{
vectorField globalSum(mesh.nPoints(), vector::zero);
labelList globalNum(mesh.nPoints(), 0);
// Note: no use of pointFaces
const faceList& faces = mesh.faces();
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
{
const face& f = faces[faceI];
const point& fc = mesh.faceCentres()[faceI];
forAll(f, fp)
{
globalSum[f[fp]] += fc;
globalNum[f[fp]]++;
}
}
// Count coupled faces as internal ones (but only once)
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patches, patchI)
{
if
(
patches[patchI].coupled()
&& refCast<const coupledPolyPatch>(patches[patchI]).owner()
)
const coupledPolyPatch& pp =
refCast<const coupledPolyPatch>(patches[patchI]);
const vectorField::subField faceCentres = pp.faceCentres();
{
const face& f = pp[i];
const point& fc = faceCentres[i];
forAll(f, fp)
{
globalSum[f[fp]] += fc;
globalNum[f[fp]]++;
}
}
}
}
syncTools::syncPointList
(
mesh,
globalSum,
plusEqOp<vector>(), // combine op
);
syncTools::syncPointList
(
mesh,
globalNum,
plusEqOp<label>(), // combine op
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
);
avgInternal.setSize(meshPoints.size());
nInternal.setSize(meshPoints.size());
forAll(avgInternal, patchPointI)
{
label meshPointI = meshPoints[patchPointI];
nInternal[patchPointI] = globalNum[meshPointI];
if (nInternal[patchPointI] == 0)
{
avgInternal[patchPointI] = globalSum[meshPointI];
}
else
{
avgInternal[patchPointI] =
globalSum[meshPointI]
/ nInternal[patchPointI];
}
}
}
// Precalculate any cell using mesh point (replacement of pointCells()[])
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList anyCell(mesh.nPoints(), -1);
forAll(mesh.faceNeighbour(), faceI)
{
label own = mesh.faceOwner()[faceI];
const face& f = mesh.faces()[faceI];
forAll(f, fp)
{
anyCell[f[fp]] = own;
}
}
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
{
label own = mesh.faceOwner()[faceI];
const face& f = mesh.faces()[faceI];
forAll(f, fp)
{
anyCell[f[fp]] = own;
}
}
// Displacement to calculate.
pointField patchDisp(meshPoints.size(), vector::zero);
forAll(pointFaces, i)
{
label meshPointI = meshPoints[i];
const point& currentPos = pp.points()[meshPointI];
// Now we have the two average points: avgBoundary and avgInternal
// and how many boundary/internal faces connect to the point
// (nBoundary, nInternal)
// Do some blending between the two.
// Note: the following section has some reasoning behind it but the
// blending factors can be experimented with.
point newPos;
if (!nonManifoldPoint.get(i))
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
{
// Points that are manifold. Weight the internal and boundary
// by their number of faces and blend with
scalar internalBlend = 0.1;
scalar blend = 0.1;
point avgPos =
(
internalBlend*nInternal[i]*avgInternal[i]
+(1-internalBlend)*nBoundary[i]*avgBoundary[i]
)
/ (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
newPos = (1-blend)*avgPos + blend*currentPos;
}
else if (nInternal[i] == 0)
{
// Non-manifold without internal faces. Use any connected cell
// as internal point instead. Use precalculated any cell to avoid
// e.g. pointCells()[meshPointI][0]
const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
scalar cellCBlend = 0.8;
scalar blend = 0.1;
point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
newPos = (1-blend)*avgPos + blend*currentPos;
}
else
{
// Non-manifold point with internal faces connected to them
scalar internalBlend = 0.9;
scalar blend = 0.1;
point avgPos =
internalBlend*avgInternal[i]
+ (1-internalBlend)*avgBoundary[i];
newPos = (1-blend)*avgPos + blend*currentPos;
}
patchDisp[i] = newPos - currentPos;
}
return patchDisp;
}
Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::edgePatchDist
(
const pointMesh& pMesh,
const indirectPrimitivePatch& pp
)
{
const polyMesh& mesh = pMesh();
// Set initial changed points to all the patch points
List<pointEdgePoint> wallInfo(pp.nPoints());
forAll(pp.localPoints(), ppI)
{
wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
}
// Current info on points
List<pointEdgePoint> allPointInfo(mesh.nPoints());
// Current info on edges
List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
PointEdgeWave<pointEdgePoint> wallCalc
(
pp.meshPoints(),
wallInfo,
allPointInfo,
allEdgeInfo,
mesh.globalData().nTotalPoints() // max iterations
);
// Copy edge values into scalarField
tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
scalarField& edgeDist = tedgeDist();
forAll(allEdgeInfo, edgeI)
{
edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
}
//{
// // For debugging: dump to file
// pointScalarField pointDist
// (
// IOobject
// (
// "pointDist",
// mesh.DB(),
// IOobject::NO_READ,
// IOobject::AUTO_WRITE
// ),
// pMesh,
// dimensionedScalar("pointDist", dimless, 0.0)
// );
//
// forAll(allEdgeInfo, edgeI)
// {
// scalar d = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
//
// const edge& e = mesh.edges()[edgeI];
//
// pointDist[e[0]] += d;
// pointDist[e[1]] += d;
// }
// forAll(pointDist, pointI)
// {
// pointDist[pointI] /= mesh.pointEdges()[pointI].size();
// }
// Info<< "Writing patch distance to " << pointDist.name()
// << " at time " << meshRefiner_.timeName() << endl;
//
// pointDist.write();
//}
return tedgeDist;
}
void Foam::autoSnapDriver::dumpMove
(
const fileName& fName,
const pointField& meshPts,
const pointField& surfPts
)
{
// Dump direction of growth into file
Pout<< nl << "Dumping move direction to " << fName << endl;
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
OFstream nearestStream(fName);
label vertI = 0;
forAll(meshPts, ptI)
{
meshTools::writeOBJ(nearestStream, meshPts[ptI]);
vertI++;
meshTools::writeOBJ(nearestStream, surfPts[ptI]);
vertI++;
nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
}
}
// Check whether all displacement vectors point outwards of patch. Return true
// if so.
bool Foam::autoSnapDriver::outwardsDisplacement
(
const indirectPrimitivePatch& pp,
const vectorField& patchDisp
)
{
const vectorField& faceNormals = pp.faceNormals();
const labelListList& pointFaces = pp.pointFaces();
forAll(pointFaces, pointI)
{
const labelList& pFaces = pointFaces[pointI];
vector disp(patchDisp[pointI]);
scalar magDisp = mag(disp);
if (magDisp > SMALL)
{
disp /= magDisp;
bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
if (!outwards)
{
Warning<< "Displacement " << patchDisp[pointI]
<< " at mesh point " << pp.meshPoints()[pointI]
<< " coord " << pp.points()[pp.meshPoints()[pointI]]
<< " points through the surrounding patch faces" << endl;
return false;
}
}
else
{
//? Displacement small but in wrong direction. Would probably be ok.
}
}
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::autoSnapDriver::autoSnapDriver
(
meshRefinement& meshRefiner,
const labelList& globalToPatch
)
:
meshRefiner_(meshRefiner),
globalToPatch_(globalToPatch)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
(
const List<labelPair>& baffles
)
{
labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
autoPtr<mapPolyMesh> map;
// No need to sync; all processors will have all same zonedSurfaces.
label nBaffles = returnReduce(baffles.size(), sumOp<label>());
if (zonedSurfaces.size() && nBaffles > 0)
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
{
// Merge any baffles
Info<< "Converting " << nBaffles << " baffles back into zoned faces ..."
<< endl;
map = meshRefiner_.mergeBaffles(baffles);
Info<< "Converted baffles in = "
<< meshRefiner_.mesh().time().cpuTimeIncrement()
<< " s\n" << nl << endl;
}
return map;
}
Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
(
const snapParameters& snapParams,
const indirectPrimitivePatch& pp
) const
{
const edgeList& edges = pp.edges();
const labelListList& pointEdges = pp.pointEdges();
const pointField& localPoints = pp.localPoints();
const fvMesh& mesh = meshRefiner_.mesh();
scalarField maxEdgeLen(localPoints.size(), -GREAT);
forAll(pointEdges, pointI)
{
const labelList& pEdges = pointEdges[pointI];
forAll(pEdges, pEdgeI)
{
const edge& e = edges[pEdges[pEdgeI]];
scalar len = e.mag(localPoints);
maxEdgeLen[pointI] = max(maxEdgeLen[pointI], len);
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
maxEdgeLen,
maxEqOp<scalar>(), // combine op
return scalarField(snapParams.snapTol()*maxEdgeLen);
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
}
void Foam::autoSnapDriver::preSmoothPatch
(
const snapParameters& snapParams,
const label nInitErrors,
const List<labelPair>& baffles,
motionSmoother& meshMover
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
labelList checkFaces;
Info<< "Smoothing patch points ..." << endl;
for
(
label smoothIter = 0;
smoothIter < snapParams.nSmoothPatch();
smoothIter++
)
{
Info<< "Smoothing iteration " << smoothIter << endl;
checkFaces.setSize(mesh.nFaces());
forAll(checkFaces, faceI)
{
checkFaces[faceI] = faceI;
}
pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
// The current mesh is the starting mesh to smooth from.
meshMover.setDisplacement(patchDisp);
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
meshMover.correct();
scalar oldErrorReduction = -1;
for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
{
Info<< nl << "Scaling iteration " << snapIter << endl;
if (snapIter == snapParams.nSnap())
{
Info<< "Displacement scaling for error reduction set to 0."
<< endl;
oldErrorReduction = meshMover.setErrorReduction(0.0);
}
// Try to adapt mesh to obtain displacement by smoothly
// decreasing displacement at error locations.
if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
{
Info<< "Successfully moved mesh" << endl;
break;
}
}
if (oldErrorReduction >= 0)
{
meshMover.setErrorReduction(oldErrorReduction);
}
Info<< endl;
}
// The current mesh is the starting mesh to smooth from.
meshMover.correct();
if (debug)
{
const_cast<Time&>(mesh.time())++;
Pout<< "Writing patch smoothed mesh to time "
<< meshRefiner_.timeName() << '.' << endl;
meshRefiner_.write
(
debug,
mesh.time().path()/meshRefiner_.timeName()
);
Pout<< "Dumped mesh in = "
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
}
Info<< "Patch points smoothed in = "
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
}
// Get (pp-local) indices of points that are both on zone and on patched surface
Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
(
const indirectPrimitivePatch& pp,
const word& zoneName
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
label zoneI = mesh.faceZones().findZoneID(zoneName);
if (zoneI == -1)
{
FatalErrorIn
(
"autoSnapDriver::getZoneSurfacePoints"
"(const indirectPrimitivePatch&, const word&)"
) << "Cannot find zone " << zoneName
<< exit(FatalError);
}
const faceZone& fZone = mesh.faceZones()[zoneI];
// Could use PrimitivePatch & localFaces to extract points but might just
// as well do it ourselves.
boolList pointOnZone(pp.nPoints(), false);
forAll(fZone, i)
{
const face& f = mesh.faces()[fZone[i]];
forAll(f, fp)
{
label meshPointI = f[fp];
Map<label>::const_iterator iter =
pp.meshPointMap().find(meshPointI);
if (iter != pp.meshPointMap().end())
{
label pointI = iter();
pointOnZone[pointI] = true;
}
}
}
return findIndices(pointOnZone, true);
}
Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
(
const scalarField& snapDist,
motionSmoother& meshMover
) const
{
Info<< "Calculating patchDisplacement as distance to nearest surface"
<< " point ..." << endl;
const indirectPrimitivePatch& pp = meshMover.patch();
const pointField& localPoints = pp.localPoints();
const refinementSurfaces& surfaces = meshRefiner_.surfaces();
const fvMesh& mesh = meshRefiner_.mesh();
// Displacement per patch point
vectorField patchDisp(localPoints.size(), vector::zero);
if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
// Current surface snapped to
labelList snapSurf(localPoints.size(), -1);
labelList zonedSurfaces =
meshRefiner_.surfaces().getNamedSurfaces();
labelList unzonedSurfaces =
meshRefiner_.surfaces().getUnnamedSurfaces();
// 1. All points to non-interface surfaces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
List<pointIndexHit> hitInfo;
labelList hitSurface;
surfaces.findNearest
(
unzonedSurfaces,
localPoints,
sqr(snapDist), // sqr of attract distance
hitSurface,
hitInfo
);
forAll(hitInfo, pointI)
if (hitInfo[pointI].hit())
{
patchDisp[pointI] =
hitInfo[pointI].hitPoint()
- localPoints[pointI];
snapSurf[pointI] = hitSurface[pointI];
// 2. All points on zones to their respective surface
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Surfaces with zone information
const wordList& faceZoneNames = surfaces.faceZoneNames();
forAll(zonedSurfaces, i)
{
label zoneSurfI = zonedSurfaces[i];
// Get indices of points both on faceZone and on pp.
labelList zonePointIndices
getZoneSurfacePoints
(
pp,
faceZoneNames[zoneSurfI]
)
);
// Find nearest for points both on faceZone and pp.
List<pointIndexHit> hitInfo;
labelList hitSurface;
surfaces.findNearest
(
labelList(1, zoneSurfI),
sqr(scalarField(minSnapDist, zonePointIndices)),
label pointI = zonePointIndices[i];
minSnapDist[pointI] = min
(
minSnapDist[pointI],
mag(patchDisp[pointI])
);
// Check if all points are being snapped
forAll(snapSurf, pointI)
{
if (snapSurf[pointI] == -1)
{
WarningIn("autoSnapDriver::calcNearestSurface(..)")
<< "For point:" << pointI
<< " coordinate:" << localPoints[pointI]
<< " did not find any surface within:"
<< minSnapDist[pointI]
<< " meter." << endl;
}
}
Info<< "Wanted displacement : average:"
<< gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
<< " min:" << gMin(magDisp)
<< " max:" << gMax(magDisp) << endl;
}
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
}
Info<< "Calculated surface displacement in = "
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
// Limit amount of movement.
forAll(patchDisp, patchPointI)
{
scalar magDisp = mag(patchDisp[patchPointI]);
if (magDisp > snapDist[patchPointI])
{
patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
Pout<< "Limiting displacement for " << patchPointI
<< " from " << magDisp << " to " << snapDist[patchPointI]
<< endl;
}
}
// Points on zones in one domain but only present as point on other
// will not do condition 2 on all. Sync explicitly.
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
patchDisp,
minMagSqrEqOp<point>(), // combine op
vector(GREAT, GREAT, GREAT) // null value (note: cannot use VGREAT)
);
return patchDisp;
}
void Foam::autoSnapDriver::smoothDisplacement
(
const snapParameters& snapParams,
motionSmoother& meshMover
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
const indirectPrimitivePatch& pp = meshMover.patch();
Info<< "Smoothing displacement ..." << endl;
// Set edge diffusivity as inverse of distance to patch
scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + SMALL));
//scalarField edgeGamma(mesh.nEdges(), 1.0);
//scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
// Get displacement field
pointVectorField& disp = meshMover.displacement();
for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
{
if ((iter % 10) == 0)
{
Info<< "Iteration " << iter << endl;
}
pointVectorField oldDisp(disp);
}
Info<< "Displacement smoothed in = "
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
if (debug)
{