62 const bool allowBoundary
65 const fvMesh&
mesh = meshRefiner_.
mesh();
66 const refinementSurfaces& surfaces = meshRefiner_.
surfaces();
68 Map<label> bafflePatch(mesh.nFaces()/1000);
70 const wordList& faceZoneNames = surfaces.faceZoneNames();
73 forAll(faceZoneNames, surfI)
75 if (faceZoneNames[surfI].size())
78 label zoneI = fZones.findZoneID(faceZoneNames[surfI]);
80 const faceZone& fZone = fZones[zoneI];
85 label patchI = globalToPatch_[surfaces.globalRegion(surfI, 0)];
88 << surfaces.names()[surfI]
89 <<
" found faceZone " << fZone.name()
90 <<
" and patch " << mesh.boundaryMesh()[patchI].name()
96 label faceI = fZone[i];
98 if (allowBoundary || mesh.isInternalFace(faceI))
100 if (!bafflePatch.insert(faceI, patchI))
102 label oldPatchI = bafflePatch[faceI];
104 if (oldPatchI != patchI)
108 <<
" fc:" << mesh.faceCentres()[faceI]
109 <<
" in zone " << fZone.name()
111 << mesh.boundaryMesh()[oldPatchI].name()
113 << mesh.boundaryMesh()[patchI].name()
127 Foam::label Foam::autoSnapDriver::getCollocatedPoints
151 label nCollocated = 0;
155 labelList firstOldPoint(newPoints.size(), -1);
156 forAll(pointMap, oldPointI)
158 label newPointI = pointMap[oldPointI];
160 if (firstOldPoint[newPointI] == -1)
163 firstOldPoint[newPointI] = oldPointI;
165 else if (firstOldPoint[newPointI] == -2)
168 isCollocatedPoint.set(oldPointI, 1u);
174 isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
177 isCollocatedPoint.set(oldPointI, 1u);
181 firstOldPoint[newPointI] = -2;
191 const motionSmoother& meshMover,
192 const List<labelPair>& baffles
199 label nNonManifoldPoints = getCollocatedPoints
205 Info<<
"Found " << nNonManifoldPoints <<
" non-mainfold point(s)."
223 const labelList& meshPoints = pp.meshPoints();
225 const polyMesh& mesh = meshMover.mesh();
233 label f0 = baffles[i].first();
234 label f1 = baffles[i].second();
236 if (isMasterFace.get(f0) == 1)
239 isMasterFace.set(f1, 0);
241 else if (isMasterFace.get(f1) == 1)
243 isMasterFace.set(f0, 0);
247 FatalErrorIn(
"autoSnapDriver::smoothPatchDisplacement(..)")
248 <<
"Both sides of baffle consisting of faces " << f0
249 <<
" and " << f1 <<
" are already slave faces."
260 labelList nBoundary(pointFaces.size(), 0);
262 forAll(pointFaces, patchPointI)
268 label faceI = pFaces[pfI];
270 if (isMasterFace.get(pp.addressing()[faceI]) == 1)
272 avgBoundary[patchPointI] += pp[faceI].centre(points);
273 nBoundary[patchPointI]++;
299 avgBoundary[i] /= nBoundary[i];
313 const faceList& faces = mesh.faces();
315 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
317 const face&
f = faces[faceI];
318 const point& fc = mesh.faceCentres()[faceI];
322 globalSum[f[fp]] += fc;
328 const polyBoundaryMesh&
patches = mesh.boundaryMesh();
334 const processorPolyPatch& pp =
335 refCast<const processorPolyPatch>(patches[patchI]);
337 if (pp.myProcNo() < pp.neighbProcNo())
343 const face& f = pp[i];
344 const point& fc = faceCentres[i];
348 globalSum[f[fp]] += fc;
354 else if (isA<cyclicPolyPatch>(patches[patchI]))
356 const cyclicPolyPatch& pp =
357 refCast<const cyclicPolyPatch>(patches[patchI]);
361 for (label i = 0; i < pp.size()/2; i++)
363 const face& f = pp[i];
364 const point& fc = faceCentres[i];
368 globalSum[f[fp]] += fc;
392 avgInternal.setSize(meshPoints.size());
393 nInternal.setSize(meshPoints.size());
395 forAll(avgInternal, patchPointI)
397 label meshPointI = meshPoints[patchPointI];
399 nInternal[patchPointI] = globalNum[meshPointI];
401 if (nInternal[patchPointI] == 0)
403 avgInternal[patchPointI] = globalSum[meshPointI];
407 avgInternal[patchPointI] =
408 globalSum[meshPointI]
409 / nInternal[patchPointI];
419 forAll(mesh.faceNeighbour(), faceI)
421 label own = mesh.faceOwner()[faceI];
422 const face& f = mesh.faces()[faceI];
426 anyCell[f[fp]] = own;
429 for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
431 label own = mesh.faceOwner()[faceI];
433 const face& f = mesh.faces()[faceI];
437 anyCell[f[fp]] = own;
447 label meshPointI = meshPoints[i];
448 const point& currentPos = pp.points()[meshPointI];
459 if (nonManifoldPoint.get(i) == 0u)
463 scalar internalBlend = 0.1;
468 internalBlend*nInternal[i]*avgInternal[i]
469 +(1-internalBlend)*nBoundary[i]*avgBoundary[i]
471 / (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
473 newPos = (1-blend)*avgPos + blend*currentPos;
475 else if (nInternal[i] == 0)
481 const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
483 scalar cellCBlend = 0.8;
486 point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
488 newPos = (1-blend)*avgPos + blend*currentPos;
493 scalar internalBlend = 0.9;
497 internalBlend*avgInternal[i]
498 + (1-internalBlend)*avgBoundary[i];
500 newPos = (1-blend)*avgPos + blend*currentPos;
503 patchDisp[i] = newPos - currentPos;
512 const pointMesh& pMesh,
516 const polyMesh& mesh = pMesh();
519 List<pointEdgePoint> wallInfo(pp.nPoints());
521 forAll(pp.localPoints(), ppI)
523 wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
527 List<pointEdgePoint> allPointInfo(mesh.nPoints());
530 List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
532 PointEdgeWave<pointEdgePoint> wallCalc
540 mesh.globalData().nTotalPoints()
544 tmp<scalarField> tedgeDist(
new scalarField(mesh.nEdges()));
547 forAll(allEdgeInfo, edgeI)
549 edgeDist[edgeI] =
Foam::sqrt(allEdgeInfo[edgeI].distSqr());
592 void Foam::autoSnapDriver::dumpMove
594 const fileName& fName,
600 Pout<<
nl <<
"Dumping move direction to " << fName <<
nl
601 <<
"View this Lightwave-OBJ file with e.g. javaview" <<
nl
604 OFstream nearestStream(fName);
616 nearestStream<<
"l " << vertI-1 <<
' ' << vertI <<
nl;
623 bool Foam::autoSnapDriver::outwardsDisplacement
632 forAll(pointFaces, pointI)
634 const labelList& pFaces = pointFaces[pointI];
636 vector disp(patchDisp[pointI]);
638 scalar magDisp =
mag(disp);
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;
666 Foam::autoSnapDriver::autoSnapDriver
672 meshRefiner_(meshRefiner),
673 globalToPatch_(globalToPatch)
684 labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
689 if (zonedSurfaces.
size())
691 fvMesh& mesh = meshRefiner_.mesh();
694 Info<<
"Converting zoned faces into baffles ..." <<
endl;
698 Map<label> faceToPatch(getZoneBafflePatches(
false));
707 ownPatch[iter.key()] = iter();
711 map = meshRefiner_.createBaffles(ownPatch, ownPatch);
720 const labelList& faceMap = map().faceMap();
721 const labelList& reverseFaceMap = map().reverseFaceMap();
725 label oldFaceI = faceMap[faceI];
730 if (iter != faceToPatch.
end())
732 label masterFaceI = reverseFaceMap[oldFaceI];
733 if (faceI != masterFaceI)
735 baffles[baffleI++] =
labelPair(masterFaceI, faceI);
740 if (baffleI != faceToPatch.
size())
743 <<
"Had " << faceToPatch.
size() <<
" patches to create "
744 <<
" but encountered " << baffleI
745 <<
" slave faces originating from patcheable faces."
752 Pout<<
"Writing baffled mesh to time "
753 << meshRefiner_.timeName() <<
endl;
757 Info<<
"Created " << nZoneFaces <<
" baffles in = "
769 labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
775 if (zonedSurfaces.
size() && nBaffles > 0)
778 Info<<
"Converting " << nBaffles <<
" baffles back into zoned faces ..."
781 map = meshRefiner_.mergeBaffles(baffles);
783 Info<<
"Converted baffles in = "
784 << meshRefiner_.mesh().time().cpuTimeIncrement()
801 const fvMesh& mesh = meshRefiner_.mesh();
805 forAll(pointEdges, pointI)
807 const labelList& pEdges = pointEdges[pointI];
811 const edge&
e = edges[pEdges[pEdgeI]];
813 scalar len = e.
mag(localPoints);
815 maxEdgeLen[pointI] =
max(maxEdgeLen[pointI], len);
829 return snapParams.
snapTol()*maxEdgeLen;
836 const label nInitErrors,
841 const fvMesh& mesh = meshRefiner_.mesh();
845 Info<<
"Smoothing patch points ..." <<
endl;
848 label smoothIter = 0;
853 Info<<
"Smoothing iteration " << smoothIter <<
endl;
857 checkFaces[faceI] = faceI;
860 pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
866 scalar oldErrorReduction = -1;
868 for (label snapIter = 0; snapIter < 2*snapParams.
nSnap(); snapIter++)
870 Info<<
nl <<
"Scaling iteration " << snapIter <<
endl;
872 if (snapIter == snapParams.
nSnap())
874 Info<<
"Displacement scaling for error reduction set to 0."
881 if (meshMover.
scaleMesh(checkFaces, baffles,
true, nInitErrors))
883 Info<<
"Successfully moved mesh" <<
endl;
888 if (oldErrorReduction >= 0)
902 Pout<<
"Writing patch smoothed mesh to time " << meshRefiner_.timeName()
908 Info<<
"Patch points smoothed in = "
920 const fvMesh& mesh = meshRefiner_.mesh();
928 "autoSnapDriver::getZoneSurfacePoints"
929 "(const indirectPrimitivePatch&, const word&)"
930 ) <<
"Cannot find zone " << zoneName
948 label meshPointI = f[fp];
955 label pointI = iter();
956 pointOnZone[pointI] =
true;
971 Info<<
"Calculating patchDisplacement as distance to nearest surface"
972 <<
" point ..." <<
endl;
977 const fvMesh& mesh = meshRefiner_.mesh();
989 meshRefiner_.surfaces().getNamedSurfaces();
991 meshRefiner_.surfaces().getUnnamedSurfaces();
1011 if (hitInfo[pointI].hit())
1014 hitInfo[pointI].hitPoint()
1015 - localPoints[pointI];
1017 snapSurf[pointI] = hitSurface[pointI];
1035 label zoneSurfI = zonedSurfaces[i];
1037 const labelList surfacesToTest(1, zoneSurfI);
1042 getZoneSurfacePoints
1045 faceZoneNames[zoneSurfI]
1063 label pointI = zonePointIndices[i];
1065 if (hitInfo[i].hit())
1068 hitInfo[i].hitPoint()
1069 - localPoints[pointI];
1071 minSnapDist[pointI] =
min
1073 minSnapDist[pointI],
1074 mag(patchDisp[pointI])
1077 snapSurf[pointI] = zoneSurfI;
1085 if (snapSurf[pointI] == -1)
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;
1099 Info<<
"Wanted displacement : average:"
1101 <<
" min:" <<
gMin(magDisp)
1102 <<
" max:" <<
gMax(magDisp) <<
endl;
1106 Info<<
"Calculated surface displacement in = "
1111 forAll(patchDisp, patchPointI)
1113 scalar magDisp =
mag(patchDisp[patchPointI]);
1115 if (magDisp > snapDist[patchPointI])
1117 patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
1119 Pout<<
"Limiting displacement for " << patchPointI
1120 <<
" from " << magDisp <<
" to " << snapDist[patchPointI]
1133 vector(GREAT, GREAT, GREAT),
1139 outwardsDisplacement(pp, patchDisp);
1150 mesh.
time().
path()/
"patchDisplacement.obj",
1166 const fvMesh& mesh = meshRefiner_.mesh();
1170 Info<<
"Smoothing displacement ..." <<
endl;
1173 scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL));
1180 for (label iter = 0; iter < snapParams.
nSmoothDispl(); iter++)
1182 if ((iter % 10) == 0)
1184 Info<<
"Iteration " << iter <<
endl;
1188 meshMover.
smooth(oldDisp, edgeGamma,
false, disp);
1190 Info<<
"Displacement smoothed in = "
1196 Pout<<
"Writing smoothed mesh to time " << meshRefiner_.timeName()
1205 Pout<<
"Writing displacement field ..." <<
endl;
1210 Pout<<
"Writing actual patch displacement ..." <<
endl;
1214 mesh.
time().
path()/
"actualPatchDisplacement.obj",
1225 const label nInitErrors,
1230 const fvMesh& mesh = meshRefiner_.mesh();
1236 scalar oldErrorReduction = -1;
1239 for (label iter = 0; iter < 2*snapParams.
nSnap(); iter++)
1243 if (iter == snapParams.
nSnap())
1245 Info<<
"Displacement scaling for error reduction set to 0." <<
endl;
1249 if (meshMover.
scaleMesh(checkFaces, baffles,
true, nInitErrors))
1251 Info<<
"Successfully moved mesh" <<
endl;
1258 Pout<<
"Writing scaled mesh to time " << meshRefiner_.timeName()
1262 Pout<<
"Writing displacement field ..." <<
endl;
1269 if (oldErrorReduction >= 0)
1273 Info<<
"Moved mesh in = "
1289 const fvMesh& mesh = meshRefiner_.mesh();
1292 Info<<
"Repatching faces according to nearest surface ..." <<
endl;
1306 labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
1307 labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
1319 label zoneSurfI = zonedSurfaces[i];
1321 label zoneI = fZones.
findZoneID(faceZoneNames[zoneSurfI]);
1323 const faceZone& fZone = fZones[zoneI];
1327 isZonedFace.set(fZone[i], 1);
1343 const scalarField snapDist(calcSnapDistance(snapParams, pp));
1347 forAll(localFaces, faceI)
1349 const face& f = localFaces[faceI];
1353 faceSnapDist[faceI] =
max
1355 faceSnapDist[faceI],
1371 sqr(4*faceSnapDist),
1379 label faceI = pp.addressing()[i];
1381 if (hitSurface[i] != -1 && (isZonedFace.get(faceI) == 0))
1383 closestPatch[i] = globalToPatch_
1410 ownPatch[pp.
start()+i] = patchI;
1411 neiPatch[pp.
start()+i] = patchI;
1418 label faceI = pp.addressing()[i];
1420 if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
1422 ownPatch[faceI] = closestPatch[i];
1423 neiPatch[faceI] = closestPatch[i];
1432 return meshRefiner_.createBaffles(ownPatch, neiPatch);
1443 fvMesh& mesh = meshRefiner_.mesh();
1446 <<
"Morphing phase" <<
nl
1447 <<
"--------------" <<
nl
1451 labelList adaptPatchIDs(meshRefiner_.meshedPatches());
1456 createZoneBaffles(baffles);
1470 const scalarField snapDist(calcSnapDistance(snapParams, pp));
1474 Info<<
"Constructing mesh displacer ..." <<
endl;
1475 Info<<
"Using mesh parameters " << motionDict <<
nl <<
endl;
1490 Info<<
"Checking initial mesh ..." <<
endl;
1499 Info<<
"Detected " << nInitErrors <<
" illegal faces"
1500 <<
" (concave, zero area or negative cell pyramid volume)"
1504 Info<<
"Checked initial mesh in = "
1508 preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
1512 calcNearestSurface(snapDist, meshMover);
1522 scaleMesh(snapParams, nInitErrors, baffles, meshMover);
1526 mergeZoneBaffles(baffles);
1529 repatchToSurface(snapParams, adaptPatchIDs);