00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "autoSnapDriver.H"
00030 #include <OpenFOAM/Time.H>
00031 #include <dynamicMesh/motionSmoother.H>
00032 #include <dynamicMesh/polyTopoChange.H>
00033 #include <OpenFOAM/OFstream.H>
00034 #include <OpenFOAM/syncTools.H>
00035 #include <finiteVolume/fvMesh.H>
00036 #include <OpenFOAM/Time.H>
00037 #include <OpenFOAM/OFstream.H>
00038 #include <OpenFOAM/mapPolyMesh.H>
00039 #include <dynamicMesh/motionSmoother.H>
00040 #include <meshTools/pointEdgePoint.H>
00041 #include <meshTools/PointEdgeWave.H>
00042 #include <OpenFOAM/mergePoints.H>
00043 #include <autoMesh/snapParameters.H>
00044 #include <autoMesh/refinementSurfaces.H>
00045
00046
00047
00048
00049 namespace Foam
00050 {
00051
00052 defineTypeNameAndDebug(autoSnapDriver, 0);
00053
00054 }
00055
00056
00057
00058
00059
00060 Foam::Map<Foam::label> Foam::autoSnapDriver::getZoneBafflePatches
00061 (
00062 const bool allowBoundary
00063 ) const
00064 {
00065 const fvMesh& mesh = meshRefiner_.mesh();
00066 const refinementSurfaces& surfaces = meshRefiner_.surfaces();
00067
00068 Map<label> bafflePatch(mesh.nFaces()/1000);
00069
00070 const wordList& faceZoneNames = surfaces.faceZoneNames();
00071 const faceZoneMesh& fZones = mesh.faceZones();
00072
00073 forAll(faceZoneNames, surfI)
00074 {
00075 if (faceZoneNames[surfI].size())
00076 {
00077
00078 label zoneI = fZones.findZoneID(faceZoneNames[surfI]);
00079
00080 const faceZone& fZone = fZones[zoneI];
00081
00083
00084
00085 label patchI = globalToPatch_[surfaces.globalRegion(surfI, 0)];
00086
00087 Info<< "For surface "
00088 << surfaces.names()[surfI]
00089 << " found faceZone " << fZone.name()
00090 << " and patch " << mesh.boundaryMesh()[patchI].name()
00091 << endl;
00092
00093
00094 forAll(fZone, i)
00095 {
00096 label faceI = fZone[i];
00097
00098 if (allowBoundary || mesh.isInternalFace(faceI))
00099 {
00100 if (!bafflePatch.insert(faceI, patchI))
00101 {
00102 label oldPatchI = bafflePatch[faceI];
00103
00104 if (oldPatchI != patchI)
00105 {
00106 FatalErrorIn("getZoneBafflePatches(const bool)")
00107 << "Face " << faceI
00108 << " fc:" << mesh.faceCentres()[faceI]
00109 << " in zone " << fZone.name()
00110 << " is in patch "
00111 << mesh.boundaryMesh()[oldPatchI].name()
00112 << " and in patch "
00113 << mesh.boundaryMesh()[patchI].name()
00114 << abort(FatalError);
00115 }
00116 }
00117 }
00118 }
00119 }
00120 }
00121 return bafflePatch;
00122 }
00123
00124
00125
00126
00127 Foam::label Foam::autoSnapDriver::getCollocatedPoints
00128 (
00129 const scalar tol,
00130 const pointField& points,
00131 PackedBoolList& isCollocatedPoint
00132 )
00133 {
00134 labelList pointMap;
00135 pointField newPoints;
00136 bool hasMerged = mergePoints
00137 (
00138 points,
00139 tol,
00140 false,
00141 pointMap,
00142 newPoints
00143 );
00144
00145 if (!returnReduce(hasMerged, orOp<bool>()))
00146 {
00147 return 0;
00148 }
00149
00150
00151 label nCollocated = 0;
00152
00153
00154
00155 labelList firstOldPoint(newPoints.size(), -1);
00156 forAll(pointMap, oldPointI)
00157 {
00158 label newPointI = pointMap[oldPointI];
00159
00160 if (firstOldPoint[newPointI] == -1)
00161 {
00162
00163 firstOldPoint[newPointI] = oldPointI;
00164 }
00165 else if (firstOldPoint[newPointI] == -2)
00166 {
00167
00168 isCollocatedPoint.set(oldPointI, 1u);
00169 nCollocated++;
00170 }
00171 else
00172 {
00173
00174 isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
00175 nCollocated++;
00176
00177 isCollocatedPoint.set(oldPointI, 1u);
00178 nCollocated++;
00179
00180
00181 firstOldPoint[newPointI] = -2;
00182 }
00183 }
00184 return returnReduce(nCollocated, sumOp<label>());
00185 }
00186
00187
00188
00189 Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
00190 (
00191 const motionSmoother& meshMover,
00192 const List<labelPair>& baffles
00193 ) const
00194 {
00195 const indirectPrimitivePatch& pp = meshMover.patch();
00196
00197
00198 PackedBoolList nonManifoldPoint(pp.nPoints());
00199 label nNonManifoldPoints = getCollocatedPoints
00200 (
00201 SMALL,
00202 pp.localPoints(),
00203 nonManifoldPoint
00204 );
00205 Info<< "Found " << nNonManifoldPoints << " non-mainfold point(s)."
00206 << endl;
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222 const labelListList& pointFaces = pp.pointFaces();
00223 const labelList& meshPoints = pp.meshPoints();
00224 const pointField& points = pp.points();
00225 const polyMesh& mesh = meshMover.mesh();
00226
00227
00228 PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
00229
00230 {
00231 forAll(baffles, i)
00232 {
00233 label f0 = baffles[i].first();
00234 label f1 = baffles[i].second();
00235
00236 if (isMasterFace.get(f0) == 1)
00237 {
00238
00239 isMasterFace.set(f1, 0);
00240 }
00241 else if (isMasterFace.get(f1) == 1)
00242 {
00243 isMasterFace.set(f0, 0);
00244 }
00245 else
00246 {
00247 FatalErrorIn("autoSnapDriver::smoothPatchDisplacement(..)")
00248 << "Both sides of baffle consisting of faces " << f0
00249 << " and " << f1 << " are already slave faces."
00250 << abort(FatalError);
00251 }
00252 }
00253 }
00254
00255
00256
00257
00258
00259 vectorField avgBoundary(pointFaces.size(), vector::zero);
00260 labelList nBoundary(pointFaces.size(), 0);
00261
00262 forAll(pointFaces, patchPointI)
00263 {
00264 const labelList& pFaces = pointFaces[patchPointI];
00265
00266 forAll(pFaces, pfI)
00267 {
00268 label faceI = pFaces[pfI];
00269
00270 if (isMasterFace.get(pp.addressing()[faceI]) == 1)
00271 {
00272 avgBoundary[patchPointI] += pp[faceI].centre(points);
00273 nBoundary[patchPointI]++;
00274 }
00275 }
00276 }
00277
00278 syncTools::syncPointList
00279 (
00280 mesh,
00281 pp.meshPoints(),
00282 avgBoundary,
00283 plusEqOp<point>(),
00284 vector::zero,
00285 false
00286 );
00287 syncTools::syncPointList
00288 (
00289 mesh,
00290 pp.meshPoints(),
00291 nBoundary,
00292 plusEqOp<label>(),
00293 0,
00294 false
00295 );
00296
00297 forAll(avgBoundary, i)
00298 {
00299 avgBoundary[i] /= nBoundary[i];
00300 }
00301
00302
00303
00304
00305
00306 vectorField avgInternal;
00307 labelList nInternal;
00308 {
00309 vectorField globalSum(mesh.nPoints(), vector::zero);
00310 labelList globalNum(mesh.nPoints(), 0);
00311
00312
00313 const faceList& faces = mesh.faces();
00314
00315 for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
00316 {
00317 const face& f = faces[faceI];
00318 const point& fc = mesh.faceCentres()[faceI];
00319
00320 forAll(f, fp)
00321 {
00322 globalSum[f[fp]] += fc;
00323 globalNum[f[fp]]++;
00324 }
00325 }
00326
00327
00328 const polyBoundaryMesh& patches = mesh.boundaryMesh();
00329
00330 forAll(patches, patchI)
00331 {
00332 if (Pstream::parRun() && isA<processorPolyPatch>(patches[patchI]))
00333 {
00334 const processorPolyPatch& pp =
00335 refCast<const processorPolyPatch>(patches[patchI]);
00336
00337 if (pp.myProcNo() < pp.neighbProcNo())
00338 {
00339 const vectorField::subField faceCentres = pp.faceCentres();
00340
00341 forAll(pp, i)
00342 {
00343 const face& f = pp[i];
00344 const point& fc = faceCentres[i];
00345
00346 forAll(f, fp)
00347 {
00348 globalSum[f[fp]] += fc;
00349 globalNum[f[fp]]++;
00350 }
00351 }
00352 }
00353 }
00354 else if (isA<cyclicPolyPatch>(patches[patchI]))
00355 {
00356 const cyclicPolyPatch& pp =
00357 refCast<const cyclicPolyPatch>(patches[patchI]);
00358
00359 const vectorField::subField faceCentres = pp.faceCentres();
00360
00361 for (label i = 0; i < pp.size()/2; i++)
00362 {
00363 const face& f = pp[i];
00364 const point& fc = faceCentres[i];
00365
00366 forAll(f, fp)
00367 {
00368 globalSum[f[fp]] += fc;
00369 globalNum[f[fp]]++;
00370 }
00371 }
00372 }
00373 }
00374
00375 syncTools::syncPointList
00376 (
00377 mesh,
00378 globalSum,
00379 plusEqOp<vector>(),
00380 vector::zero,
00381 false
00382 );
00383 syncTools::syncPointList
00384 (
00385 mesh,
00386 globalNum,
00387 plusEqOp<label>(),
00388 0,
00389 false
00390 );
00391
00392 avgInternal.setSize(meshPoints.size());
00393 nInternal.setSize(meshPoints.size());
00394
00395 forAll(avgInternal, patchPointI)
00396 {
00397 label meshPointI = meshPoints[patchPointI];
00398
00399 nInternal[patchPointI] = globalNum[meshPointI];
00400
00401 if (nInternal[patchPointI] == 0)
00402 {
00403 avgInternal[patchPointI] = globalSum[meshPointI];
00404 }
00405 else
00406 {
00407 avgInternal[patchPointI] =
00408 globalSum[meshPointI]
00409 / nInternal[patchPointI];
00410 }
00411 }
00412 }
00413
00414
00415
00416
00417
00418 labelList anyCell(mesh.nPoints(), -1);
00419 forAll(mesh.faceNeighbour(), faceI)
00420 {
00421 label own = mesh.faceOwner()[faceI];
00422 const face& f = mesh.faces()[faceI];
00423
00424 forAll(f, fp)
00425 {
00426 anyCell[f[fp]] = own;
00427 }
00428 }
00429 for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
00430 {
00431 label own = mesh.faceOwner()[faceI];
00432
00433 const face& f = mesh.faces()[faceI];
00434
00435 forAll(f, fp)
00436 {
00437 anyCell[f[fp]] = own;
00438 }
00439 }
00440
00441
00442
00443 pointField patchDisp(meshPoints.size(), vector::zero);
00444
00445 forAll(pointFaces, i)
00446 {
00447 label meshPointI = meshPoints[i];
00448 const point& currentPos = pp.points()[meshPointI];
00449
00450
00451
00452
00453
00454
00455
00456
00457 point newPos;
00458
00459 if (nonManifoldPoint.get(i) == 0u)
00460 {
00461
00462
00463 scalar internalBlend = 0.1;
00464 scalar blend = 0.1;
00465
00466 point avgPos =
00467 (
00468 internalBlend*nInternal[i]*avgInternal[i]
00469 +(1-internalBlend)*nBoundary[i]*avgBoundary[i]
00470 )
00471 / (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
00472
00473 newPos = (1-blend)*avgPos + blend*currentPos;
00474 }
00475 else if (nInternal[i] == 0)
00476 {
00477
00478
00479
00480
00481 const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
00482
00483 scalar cellCBlend = 0.8;
00484 scalar blend = 0.1;
00485
00486 point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
00487
00488 newPos = (1-blend)*avgPos + blend*currentPos;
00489 }
00490 else
00491 {
00492
00493 scalar internalBlend = 0.9;
00494 scalar blend = 0.1;
00495
00496 point avgPos =
00497 internalBlend*avgInternal[i]
00498 + (1-internalBlend)*avgBoundary[i];
00499
00500 newPos = (1-blend)*avgPos + blend*currentPos;
00501 }
00502
00503 patchDisp[i] = newPos - currentPos;
00504 }
00505
00506 return patchDisp;
00507 }
00508
00509
00510 Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::edgePatchDist
00511 (
00512 const pointMesh& pMesh,
00513 const indirectPrimitivePatch& pp
00514 )
00515 {
00516 const polyMesh& mesh = pMesh();
00517
00518
00519 List<pointEdgePoint> wallInfo(pp.nPoints());
00520
00521 forAll(pp.localPoints(), ppI)
00522 {
00523 wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
00524 }
00525
00526
00527 List<pointEdgePoint> allPointInfo(mesh.nPoints());
00528
00529
00530 List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
00531
00532 PointEdgeWave<pointEdgePoint> wallCalc
00533 (
00534 mesh,
00535 pp.meshPoints(),
00536 wallInfo,
00537
00538 allPointInfo,
00539 allEdgeInfo,
00540 mesh.globalData().nTotalPoints()
00541 );
00542
00543
00544 tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
00545 scalarField& edgeDist = tedgeDist();
00546
00547 forAll(allEdgeInfo, edgeI)
00548 {
00549 edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
00550 }
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588 return tedgeDist;
00589 }
00590
00591
00592 void Foam::autoSnapDriver::dumpMove
00593 (
00594 const fileName& fName,
00595 const pointField& meshPts,
00596 const pointField& surfPts
00597 )
00598 {
00599
00600 Pout<< nl << "Dumping move direction to " << fName << nl
00601 << "View this Lightwave-OBJ file with e.g. javaview" << nl
00602 << endl;
00603
00604 OFstream nearestStream(fName);
00605
00606 label vertI = 0;
00607
00608 forAll(meshPts, ptI)
00609 {
00610 meshTools::writeOBJ(nearestStream, meshPts[ptI]);
00611 vertI++;
00612
00613 meshTools::writeOBJ(nearestStream, surfPts[ptI]);
00614 vertI++;
00615
00616 nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
00617 }
00618 }
00619
00620
00621
00622
00623 bool Foam::autoSnapDriver::outwardsDisplacement
00624 (
00625 const indirectPrimitivePatch& pp,
00626 const vectorField& patchDisp
00627 )
00628 {
00629 const vectorField& faceNormals = pp.faceNormals();
00630 const labelListList& pointFaces = pp.pointFaces();
00631
00632 forAll(pointFaces, pointI)
00633 {
00634 const labelList& pFaces = pointFaces[pointI];
00635
00636 vector disp(patchDisp[pointI]);
00637
00638 scalar magDisp = mag(disp);
00639
00640 if (magDisp > SMALL)
00641 {
00642 disp /= magDisp;
00643
00644 bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
00645
00646 if (!outwards)
00647 {
00648 Warning<< "Displacement " << patchDisp[pointI]
00649 << " at mesh point " << pp.meshPoints()[pointI]
00650 << " coord " << pp.points()[pp.meshPoints()[pointI]]
00651 << " points through the surrounding patch faces" << endl;
00652 return false;
00653 }
00654 }
00655 else
00656 {
00657
00658 }
00659 }
00660 return true;
00661 }
00662
00663
00664
00665
00666 Foam::autoSnapDriver::autoSnapDriver
00667 (
00668 meshRefinement& meshRefiner,
00669 const labelList& globalToPatch
00670 )
00671 :
00672 meshRefiner_(meshRefiner),
00673 globalToPatch_(globalToPatch)
00674 {}
00675
00676
00677
00678
00679 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::createZoneBaffles
00680 (
00681 List<labelPair>& baffles
00682 )
00683 {
00684 labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
00685
00686 autoPtr<mapPolyMesh> map;
00687
00688
00689 if (zonedSurfaces.size())
00690 {
00691 fvMesh& mesh = meshRefiner_.mesh();
00692
00693
00694 Info<< "Converting zoned faces into baffles ..." << endl;
00695
00696
00697
00698 Map<label> faceToPatch(getZoneBafflePatches(false));
00699
00700 label nZoneFaces = returnReduce(faceToPatch.size(), sumOp<label>());
00701 if (nZoneFaces > 0)
00702 {
00703
00704 labelList ownPatch(mesh.nFaces(), -1);
00705 forAllConstIter(Map<label>, faceToPatch, iter)
00706 {
00707 ownPatch[iter.key()] = iter();
00708 }
00709
00710
00711 map = meshRefiner_.createBaffles(ownPatch, ownPatch);
00712
00713
00714
00715
00716
00717 baffles.setSize(faceToPatch.size());
00718 label baffleI = 0;
00719
00720 const labelList& faceMap = map().faceMap();
00721 const labelList& reverseFaceMap = map().reverseFaceMap();
00722
00723 forAll(faceMap, faceI)
00724 {
00725 label oldFaceI = faceMap[faceI];
00726
00727
00728 Map<label>::const_iterator iter = faceToPatch.find(oldFaceI);
00729
00730 if (iter != faceToPatch.end())
00731 {
00732 label masterFaceI = reverseFaceMap[oldFaceI];
00733 if (faceI != masterFaceI)
00734 {
00735 baffles[baffleI++] = labelPair(masterFaceI, faceI);
00736 }
00737 }
00738 }
00739
00740 if (baffleI != faceToPatch.size())
00741 {
00742 FatalErrorIn("autoSnapDriver::createZoneBaffles(..)")
00743 << "Had " << faceToPatch.size() << " patches to create "
00744 << " but encountered " << baffleI
00745 << " slave faces originating from patcheable faces."
00746 << abort(FatalError);
00747 }
00748
00749 if (debug)
00750 {
00751 const_cast<Time&>(mesh.time())++;
00752 Pout<< "Writing baffled mesh to time "
00753 << meshRefiner_.timeName() << endl;
00754 mesh.write();
00755 }
00756 }
00757 Info<< "Created " << nZoneFaces << " baffles in = "
00758 << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
00759 }
00760 return map;
00761 }
00762
00763
00764 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
00765 (
00766 const List<labelPair>& baffles
00767 )
00768 {
00769 labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
00770
00771 autoPtr<mapPolyMesh> map;
00772
00773
00774 label nBaffles = returnReduce(baffles.size(), sumOp<label>());
00775 if (zonedSurfaces.size() && nBaffles > 0)
00776 {
00777
00778 Info<< "Converting " << nBaffles << " baffles back into zoned faces ..."
00779 << endl;
00780
00781 map = meshRefiner_.mergeBaffles(baffles);
00782
00783 Info<< "Converted baffles in = "
00784 << meshRefiner_.mesh().time().cpuTimeIncrement()
00785 << " s\n" << nl << endl;
00786 }
00787
00788 return map;
00789 }
00790
00791
00792 Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
00793 (
00794 const snapParameters& snapParams,
00795 const indirectPrimitivePatch& pp
00796 ) const
00797 {
00798 const edgeList& edges = pp.edges();
00799 const labelListList& pointEdges = pp.pointEdges();
00800 const pointField& localPoints = pp.localPoints();
00801 const fvMesh& mesh = meshRefiner_.mesh();
00802
00803 scalarField maxEdgeLen(localPoints.size(), -GREAT);
00804
00805 forAll(pointEdges, pointI)
00806 {
00807 const labelList& pEdges = pointEdges[pointI];
00808
00809 forAll(pEdges, pEdgeI)
00810 {
00811 const edge& e = edges[pEdges[pEdgeI]];
00812
00813 scalar len = e.mag(localPoints);
00814
00815 maxEdgeLen[pointI] = max(maxEdgeLen[pointI], len);
00816 }
00817 }
00818
00819 syncTools::syncPointList
00820 (
00821 mesh,
00822 pp.meshPoints(),
00823 maxEdgeLen,
00824 maxEqOp<scalar>(),
00825 -GREAT,
00826 false
00827 );
00828
00829 return snapParams.snapTol()*maxEdgeLen;
00830 }
00831
00832
00833 void Foam::autoSnapDriver::preSmoothPatch
00834 (
00835 const snapParameters& snapParams,
00836 const label nInitErrors,
00837 const List<labelPair>& baffles,
00838 motionSmoother& meshMover
00839 ) const
00840 {
00841 const fvMesh& mesh = meshRefiner_.mesh();
00842
00843 labelList checkFaces;
00844
00845 Info<< "Smoothing patch points ..." << endl;
00846 for
00847 (
00848 label smoothIter = 0;
00849 smoothIter < snapParams.nSmoothPatch();
00850 smoothIter++
00851 )
00852 {
00853 Info<< "Smoothing iteration " << smoothIter << endl;
00854 checkFaces.setSize(mesh.nFaces());
00855 forAll(checkFaces, faceI)
00856 {
00857 checkFaces[faceI] = faceI;
00858 }
00859
00860 pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
00861
00862
00863 meshMover.setDisplacement(patchDisp);
00864 meshMover.correct();
00865
00866 scalar oldErrorReduction = -1;
00867
00868 for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
00869 {
00870 Info<< nl << "Scaling iteration " << snapIter << endl;
00871
00872 if (snapIter == snapParams.nSnap())
00873 {
00874 Info<< "Displacement scaling for error reduction set to 0."
00875 << endl;
00876 oldErrorReduction = meshMover.setErrorReduction(0.0);
00877 }
00878
00879
00880
00881 if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
00882 {
00883 Info<< "Successfully moved mesh" << endl;
00884 break;
00885 }
00886 }
00887
00888 if (oldErrorReduction >= 0)
00889 {
00890 meshMover.setErrorReduction(oldErrorReduction);
00891 }
00892 Info<< endl;
00893 }
00894
00895
00896
00897 meshMover.correct();
00898
00899 if (debug)
00900 {
00901 const_cast<Time&>(mesh.time())++;
00902 Pout<< "Writing patch smoothed mesh to time " << meshRefiner_.timeName()
00903 << endl;
00904
00905 mesh.write();
00906 }
00907
00908 Info<< "Patch points smoothed in = "
00909 << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
00910 }
00911
00912
00913
00914 Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
00915 (
00916 const indirectPrimitivePatch& pp,
00917 const word& zoneName
00918 ) const
00919 {
00920 const fvMesh& mesh = meshRefiner_.mesh();
00921
00922 label zoneI = mesh.faceZones().findZoneID(zoneName);
00923
00924 if (zoneI == -1)
00925 {
00926 FatalErrorIn
00927 (
00928 "autoSnapDriver::getZoneSurfacePoints"
00929 "(const indirectPrimitivePatch&, const word&)"
00930 ) << "Cannot find zone " << zoneName
00931 << exit(FatalError);
00932 }
00933
00934 const faceZone& fZone = mesh.faceZones()[zoneI];
00935
00936
00937
00938
00939
00940 boolList pointOnZone(pp.nPoints(), false);
00941
00942 forAll(fZone, i)
00943 {
00944 const face& f = mesh.faces()[fZone[i]];
00945
00946 forAll(f, fp)
00947 {
00948 label meshPointI = f[fp];
00949
00950 Map<label>::const_iterator iter =
00951 pp.meshPointMap().find(meshPointI);
00952
00953 if (iter != pp.meshPointMap().end())
00954 {
00955 label pointI = iter();
00956 pointOnZone[pointI] = true;
00957 }
00958 }
00959 }
00960
00961 return findIndices(pointOnZone, true);
00962 }
00963
00964
00965 Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
00966 (
00967 const scalarField& snapDist,
00968 motionSmoother& meshMover
00969 ) const
00970 {
00971 Info<< "Calculating patchDisplacement as distance to nearest surface"
00972 << " point ..." << endl;
00973
00974 const indirectPrimitivePatch& pp = meshMover.patch();
00975 const pointField& localPoints = pp.localPoints();
00976 const refinementSurfaces& surfaces = meshRefiner_.surfaces();
00977 const fvMesh& mesh = meshRefiner_.mesh();
00978
00979
00980 vectorField patchDisp(localPoints.size(), vector::zero);
00981
00982 if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
00983 {
00984
00985 labelList snapSurf(localPoints.size(), -1);
00986
00987
00988 labelList zonedSurfaces =
00989 meshRefiner_.surfaces().getNamedSurfaces();
00990 labelList unzonedSurfaces =
00991 meshRefiner_.surfaces().getUnnamedSurfaces();
00992
00993
00994
00995
00996
00997 {
00998 List<pointIndexHit> hitInfo;
00999 labelList hitSurface;
01000 surfaces.findNearest
01001 (
01002 unzonedSurfaces,
01003 localPoints,
01004 sqr(4*snapDist),
01005 hitSurface,
01006 hitInfo
01007 );
01008
01009 forAll(hitInfo, pointI)
01010 {
01011 if (hitInfo[pointI].hit())
01012 {
01013 patchDisp[pointI] =
01014 hitInfo[pointI].hitPoint()
01015 - localPoints[pointI];
01016
01017 snapSurf[pointI] = hitSurface[pointI];
01018 }
01019 }
01020 }
01021
01022
01023
01024
01025
01026
01027
01028 const wordList& faceZoneNames = surfaces.faceZoneNames();
01029
01030
01031 scalarField minSnapDist(snapDist);
01032
01033 forAll(zonedSurfaces, i)
01034 {
01035 label zoneSurfI = zonedSurfaces[i];
01036
01037 const labelList surfacesToTest(1, zoneSurfI);
01038
01039
01040 labelList zonePointIndices
01041 (
01042 getZoneSurfacePoints
01043 (
01044 pp,
01045 faceZoneNames[zoneSurfI]
01046 )
01047 );
01048
01049
01050 List<pointIndexHit> hitInfo;
01051 labelList hitSurface;
01052 surfaces.findNearest
01053 (
01054 labelList(1, zoneSurfI),
01055 pointField(localPoints, zonePointIndices),
01056 sqr(4*scalarField(minSnapDist, zonePointIndices)),
01057 hitSurface,
01058 hitInfo
01059 );
01060
01061 forAll(hitInfo, i)
01062 {
01063 label pointI = zonePointIndices[i];
01064
01065 if (hitInfo[i].hit())
01066 {
01067 patchDisp[pointI] =
01068 hitInfo[i].hitPoint()
01069 - localPoints[pointI];
01070
01071 minSnapDist[pointI] = min
01072 (
01073 minSnapDist[pointI],
01074 mag(patchDisp[pointI])
01075 );
01076
01077 snapSurf[pointI] = zoneSurfI;
01078 }
01079 }
01080 }
01081
01082
01083 forAll(snapSurf, pointI)
01084 {
01085 if (snapSurf[pointI] == -1)
01086 {
01087 WarningIn("autoSnapDriver::calcNearestSurface(..)")
01088 << "For point:" << pointI
01089 << " coordinate:" << localPoints[pointI]
01090 << " did not find any surface within:"
01091 << minSnapDist[pointI]
01092 << " meter." << endl;
01093 }
01094 }
01095
01096 {
01097 scalarField magDisp(mag(patchDisp));
01098
01099 Info<< "Wanted displacement : average:"
01100 << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
01101 << " min:" << gMin(magDisp)
01102 << " max:" << gMax(magDisp) << endl;
01103 }
01104 }
01105
01106 Info<< "Calculated surface displacement in = "
01107 << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
01108
01109
01110
01111 forAll(patchDisp, patchPointI)
01112 {
01113 scalar magDisp = mag(patchDisp[patchPointI]);
01114
01115 if (magDisp > snapDist[patchPointI])
01116 {
01117 patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
01118
01119 Pout<< "Limiting displacement for " << patchPointI
01120 << " from " << magDisp << " to " << snapDist[patchPointI]
01121 << endl;
01122 }
01123 }
01124
01125
01126
01127 syncTools::syncPointList
01128 (
01129 mesh,
01130 pp.meshPoints(),
01131 patchDisp,
01132 minMagEqOp(),
01133 vector(GREAT, GREAT, GREAT),
01134 false
01135 );
01136
01137
01138
01139 outwardsDisplacement(pp, patchDisp);
01140
01141
01142
01143
01144 meshMover.setDisplacement(patchDisp);
01145
01146 if (debug)
01147 {
01148 dumpMove
01149 (
01150 mesh.time().path()/"patchDisplacement.obj",
01151 pp.localPoints(),
01152 pp.localPoints() + patchDisp
01153 );
01154 }
01155
01156 return patchDisp;
01157 }
01158
01159
01160 void Foam::autoSnapDriver::smoothDisplacement
01161 (
01162 const snapParameters& snapParams,
01163 motionSmoother& meshMover
01164 ) const
01165 {
01166 const fvMesh& mesh = meshRefiner_.mesh();
01167 const pointMesh& pMesh = meshMover.pMesh();
01168 const indirectPrimitivePatch& pp = meshMover.patch();
01169
01170 Info<< "Smoothing displacement ..." << endl;
01171
01172
01173 scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL));
01174
01175
01176
01177
01178 pointVectorField& disp = meshMover.displacement();
01179
01180 for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
01181 {
01182 if ((iter % 10) == 0)
01183 {
01184 Info<< "Iteration " << iter << endl;
01185 }
01186 pointVectorField oldDisp(disp);
01187
01188 meshMover.smooth(oldDisp, edgeGamma, false, disp);
01189 }
01190 Info<< "Displacement smoothed in = "
01191 << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
01192
01193 if (debug)
01194 {
01195 const_cast<Time&>(mesh.time())++;
01196 Pout<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
01197 << endl;
01198
01199
01200
01201
01202
01203 mesh.write();
01204
01205 Pout<< "Writing displacement field ..." << endl;
01206 disp.write();
01207 tmp<pointScalarField> magDisp(mag(disp));
01208 magDisp().write();
01209
01210 Pout<< "Writing actual patch displacement ..." << endl;
01211 vectorField actualPatchDisp(disp, pp.meshPoints());
01212 dumpMove
01213 (
01214 mesh.time().path()/"actualPatchDisplacement.obj",
01215 pp.localPoints(),
01216 pp.localPoints() + actualPatchDisp
01217 );
01218 }
01219 }
01220
01221
01222 void Foam::autoSnapDriver::scaleMesh
01223 (
01224 const snapParameters& snapParams,
01225 const label nInitErrors,
01226 const List<labelPair>& baffles,
01227 motionSmoother& meshMover
01228 )
01229 {
01230 const fvMesh& mesh = meshRefiner_.mesh();
01231
01232
01233
01234 labelList checkFaces(identity(mesh.nFaces()));
01235
01236 scalar oldErrorReduction = -1;
01237
01238 Info<< "Moving mesh ..." << endl;
01239 for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
01240 {
01241 Info<< nl << "Iteration " << iter << endl;
01242
01243 if (iter == snapParams.nSnap())
01244 {
01245 Info<< "Displacement scaling for error reduction set to 0." << endl;
01246 oldErrorReduction = meshMover.setErrorReduction(0.0);
01247 }
01248
01249 if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
01250 {
01251 Info<< "Successfully moved mesh" << endl;
01252
01253 break;
01254 }
01255 if (debug)
01256 {
01257 const_cast<Time&>(mesh.time())++;
01258 Pout<< "Writing scaled mesh to time " << meshRefiner_.timeName()
01259 << endl;
01260 mesh.write();
01261
01262 Pout<< "Writing displacement field ..." << endl;
01263 meshMover.displacement().write();
01264 tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
01265 magDisp().write();
01266 }
01267 }
01268
01269 if (oldErrorReduction >= 0)
01270 {
01271 meshMover.setErrorReduction(oldErrorReduction);
01272 }
01273 Info<< "Moved mesh in = "
01274 << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
01275 }
01276
01277
01278
01279
01280
01281
01282
01283 Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface
01284 (
01285 const snapParameters& snapParams,
01286 const labelList& adaptPatchIDs
01287 )
01288 {
01289 const fvMesh& mesh = meshRefiner_.mesh();
01290 const refinementSurfaces& surfaces = meshRefiner_.surfaces();
01291
01292 Info<< "Repatching faces according to nearest surface ..." << endl;
01293
01294
01295 autoPtr<indirectPrimitivePatch> ppPtr
01296 (
01297 meshRefinement::makePatch
01298 (
01299 mesh,
01300 adaptPatchIDs
01301 )
01302 );
01303 indirectPrimitivePatch& pp = ppPtr();
01304
01305
01306 labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
01307 labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
01308
01309
01310
01311 PackedBoolList isZonedFace(mesh.nFaces(), 0);
01312 {
01313
01314 const wordList& faceZoneNames = surfaces.faceZoneNames();
01315 const faceZoneMesh& fZones = mesh.faceZones();
01316
01317 forAll(zonedSurfaces, i)
01318 {
01319 label zoneSurfI = zonedSurfaces[i];
01320
01321 label zoneI = fZones.findZoneID(faceZoneNames[zoneSurfI]);
01322
01323 const faceZone& fZone = fZones[zoneI];
01324
01325 forAll(fZone, i)
01326 {
01327 isZonedFace.set(fZone[i], 1);
01328 }
01329 }
01330 }
01331
01332
01333
01334
01335
01336
01337 labelList closestPatch(pp.size(), -1);
01338 {
01339
01340 scalarField faceSnapDist(pp.size(), -GREAT);
01341 {
01342
01343 const scalarField snapDist(calcSnapDistance(snapParams, pp));
01344
01345 const faceList& localFaces = pp.localFaces();
01346
01347 forAll(localFaces, faceI)
01348 {
01349 const face& f = localFaces[faceI];
01350
01351 forAll(f, fp)
01352 {
01353 faceSnapDist[faceI] = max
01354 (
01355 faceSnapDist[faceI],
01356 snapDist[f[fp]]
01357 );
01358 }
01359 }
01360 }
01361
01362 pointField localFaceCentres(mesh.faceCentres(), pp.addressing());
01363
01364
01365 labelList hitSurface;
01366 labelList hitRegion;
01367 surfaces.findNearestRegion
01368 (
01369 unzonedSurfaces,
01370 localFaceCentres,
01371 sqr(4*faceSnapDist),
01372 hitSurface,
01373 hitRegion
01374 );
01375
01376
01377 forAll(pp, i)
01378 {
01379 label faceI = pp.addressing()[i];
01380
01381 if (hitSurface[i] != -1 && (isZonedFace.get(faceI) == 0))
01382 {
01383 closestPatch[i] = globalToPatch_
01384 [
01385 surfaces.globalRegion
01386 (
01387 hitSurface[i],
01388 hitRegion[i]
01389 )
01390 ];
01391 }
01392 }
01393 }
01394
01395
01396
01397
01398
01399 labelList ownPatch(mesh.nFaces(), -1);
01400 labelList neiPatch(mesh.nFaces(), -1);
01401
01402 const polyBoundaryMesh& patches = mesh.boundaryMesh();
01403
01404 forAll(patches, patchI)
01405 {
01406 const polyPatch& pp = patches[patchI];
01407
01408 forAll(pp, i)
01409 {
01410 ownPatch[pp.start()+i] = patchI;
01411 neiPatch[pp.start()+i] = patchI;
01412 }
01413 }
01414
01415 label nChanged = 0;
01416 forAll(closestPatch, i)
01417 {
01418 label faceI = pp.addressing()[i];
01419
01420 if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
01421 {
01422 ownPatch[faceI] = closestPatch[i];
01423 neiPatch[faceI] = closestPatch[i];
01424 nChanged++;
01425 }
01426 }
01427
01428 Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
01429 << " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
01430 << endl;
01431
01432 return meshRefiner_.createBaffles(ownPatch, neiPatch);
01433 }
01434
01435
01436 void Foam::autoSnapDriver::doSnap
01437 (
01438 const dictionary& snapDict,
01439 const dictionary& motionDict,
01440 const snapParameters& snapParams
01441 )
01442 {
01443 fvMesh& mesh = meshRefiner_.mesh();
01444
01445 Info<< nl
01446 << "Morphing phase" << nl
01447 << "--------------" << nl
01448 << endl;
01449
01450
01451 labelList adaptPatchIDs(meshRefiner_.meshedPatches());
01452
01453
01454
01455 List<labelPair> baffles;
01456 createZoneBaffles(baffles);
01457
01458 {
01459 autoPtr<indirectPrimitivePatch> ppPtr
01460 (
01461 meshRefinement::makePatch
01462 (
01463 mesh,
01464 adaptPatchIDs
01465 )
01466 );
01467 indirectPrimitivePatch& pp = ppPtr();
01468
01469
01470 const scalarField snapDist(calcSnapDistance(snapParams, pp));
01471
01472
01473
01474 Info<< "Constructing mesh displacer ..." << endl;
01475 Info<< "Using mesh parameters " << motionDict << nl << endl;
01476
01477 const pointMesh& pMesh = pointMesh::New(mesh);
01478
01479 motionSmoother meshMover
01480 (
01481 mesh,
01482 pp,
01483 adaptPatchIDs,
01484 meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
01485 motionDict
01486 );
01487
01488
01489
01490 Info<< "Checking initial mesh ..." << endl;
01491 labelHashSet wrongFaces(mesh.nFaces()/100);
01492 motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
01493 const label nInitErrors = returnReduce
01494 (
01495 wrongFaces.size(),
01496 sumOp<label>()
01497 );
01498
01499 Info<< "Detected " << nInitErrors << " illegal faces"
01500 << " (concave, zero area or negative cell pyramid volume)"
01501 << endl;
01502
01503
01504 Info<< "Checked initial mesh in = "
01505 << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
01506
01507
01508 preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
01509
01510
01511
01512 calcNearestSurface(snapDist, meshMover);
01513
01515
01516
01517
01518
01519
01520
01521
01522 scaleMesh(snapParams, nInitErrors, baffles, meshMover);
01523 }
01524
01525
01526 mergeZoneBaffles(baffles);
01527
01528
01529 repatchToSurface(snapParams, adaptPatchIDs);
01530 }
01531
01532
01533