FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

autoSnapDriver.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 Description
00025     All to do with snapping to the surface
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 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00048 
00049 namespace Foam
00050 {
00051 
00052 defineTypeNameAndDebug(autoSnapDriver, 0);
00053 
00054 } // End namespace Foam
00055 
00056 
00057 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00058 
00059 // Get faces to repatch. Returns map from face to patch.
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             // Get zone
00078             label zoneI = fZones.findZoneID(faceZoneNames[surfI]);
00079 
00080             const faceZone& fZone = fZones[zoneI];
00081 
00083             //label patchI = surfaceToCyclicPatch_[surfI];
00084             // Get patch of (first region) of surface
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 // Calculate geometrically collocated points, Requires PackedList to be
00126 // sized and initalised!
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,                         // points
00139         tol,                            // mergeTol
00140         false,                          // verbose
00141         pointMap,
00142         newPoints
00143     );
00144 
00145     if (!returnReduce(hasMerged, orOp<bool>()))
00146     {
00147         return 0;
00148     }
00149 
00150     // Determine which newPoints are referenced more than once
00151     label nCollocated = 0;
00152 
00153     // Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
00154     // twice)
00155     labelList firstOldPoint(newPoints.size(), -1);
00156     forAll(pointMap, oldPointI)
00157     {
00158         label newPointI = pointMap[oldPointI];
00159 
00160         if (firstOldPoint[newPointI] == -1)
00161         {
00162             // First use of oldPointI. Store.
00163             firstOldPoint[newPointI] = oldPointI;
00164         }
00165         else if (firstOldPoint[newPointI] == -2)
00166         {
00167             // Third or more reference of oldPointI -> non-manifold
00168             isCollocatedPoint.set(oldPointI, 1u);
00169             nCollocated++;
00170         }
00171         else
00172         {
00173             // Second reference of oldPointI -> non-manifold
00174             isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
00175             nCollocated++;
00176 
00177             isCollocatedPoint.set(oldPointI, 1u);
00178             nCollocated++;
00179 
00180             // Mark with special value to save checking next time round
00181             firstOldPoint[newPointI] = -2;
00182         }
00183     }
00184     return returnReduce(nCollocated, sumOp<label>());
00185 }
00186 
00187 
00188 // Calculate displacement as average of patch points.
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     // Calculate geometrically non-manifold points on the patch to be moved.
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     // Average points
00210     // ~~~~~~~~~~~~~~
00211 
00212     // We determine three points:
00213     // - average of (centres of) connected patch faces
00214     // - average of (centres of) connected internal mesh faces
00215     // - as fallback: centre of any connected cell
00216     // so we can do something moderately sensible for non/manifold points.
00217 
00218     // Note: the averages are calculated properly parallel. This is
00219     // necessary to get the points shared by processors correct.
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     // Get labels of faces to count (master of coupled faces and baffle pairs)
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                 // Make f1 a slave
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     // Get average position of boundary face centres
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>(),  // combine op
00284         vector::zero,       // null value
00285         false               // no separation
00286     );
00287     syncTools::syncPointList
00288     (
00289         mesh,
00290         pp.meshPoints(),
00291         nBoundary,
00292         plusEqOp<label>(),  // combine op
00293         0,                  // null value
00294         false               // no separation
00295     );
00296 
00297     forAll(avgBoundary, i)
00298     {
00299         avgBoundary[i] /= nBoundary[i];
00300     }
00301 
00302 
00303     // Get average position of internal face centres
00304     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00305 
00306     vectorField avgInternal;
00307     labelList nInternal;
00308     {
00309         vectorField globalSum(mesh.nPoints(), vector::zero);
00310         labelList globalNum(mesh.nPoints(), 0);
00311 
00312         // Note: no use of pointFaces
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         // Count coupled faces as internal ones (but only once)
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>(), // combine op
00380             vector::zero,       // null value
00381             false               // no separation
00382         );
00383         syncTools::syncPointList
00384         (
00385             mesh,
00386             globalNum,
00387             plusEqOp<label>(),  // combine op
00388             0,                  // null value
00389             false               // no separation
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     // Precalculate any cell using mesh point (replacement of pointCells()[])
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     // Displacement to calculate.
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         // Now we have the two average points: avgBoundary and avgInternal
00451         // and how many boundary/internal faces connect to the point
00452         // (nBoundary, nInternal)
00453         // Do some blending between the two.
00454         // Note: the following section has some reasoning behind it but the
00455         // blending factors can be experimented with.
00456 
00457         point newPos;
00458 
00459         if (nonManifoldPoint.get(i) == 0u)
00460         {
00461             // Points that are manifold. Weight the internal and boundary
00462             // by their number of faces and blend with
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             // Non-manifold without internal faces. Use any connected cell
00478             // as internal point instead. Use precalculated any cell to avoid
00479             // e.g. pointCells()[meshPointI][0]
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             // Non-manifold point with internal faces connected to them
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     // Set initial changed points to all the patch points
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     // Current info on points
00527     List<pointEdgePoint> allPointInfo(mesh.nPoints());
00528 
00529     // Current info on edges
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()  // max iterations
00541     );
00542 
00543     // Copy edge values into scalarField
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     //    // For debugging: dump to file
00555     //    pointScalarField pointDist
00556     //    (
00557     //        IOobject
00558     //        (
00559     //            "pointDist",
00560     //            meshRefiner_.timeName(),
00561     //            mesh.DB(),
00562     //            IOobject::NO_READ,
00563     //            IOobject::AUTO_WRITE
00564     //        ),
00565     //        pMesh,
00566     //        dimensionedScalar("pointDist", dimless, 0.0)
00567     //    );
00568     //
00569     //    forAll(allEdgeInfo, edgeI)
00570     //    {
00571     //        scalar d = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
00572     //
00573     //        const edge& e = mesh.edges()[edgeI];
00574     //
00575     //        pointDist[e[0]] += d;
00576     //        pointDist[e[1]] += d;
00577     //    }
00578     //    forAll(pointDist, pointI)
00579     //    {
00580     //        pointDist[pointI] /= mesh.pointEdges()[pointI].size();
00581     //    }
00582     //    Info<< "Writing patch distance to " << pointDist.name()
00583     //        << " at time " << meshRefiner_.timeName() << endl;
00584     //
00585     //    pointDist.write();
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     // Dump direction of growth into file
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 // Check whether all displacement vectors point outwards of patch. Return true
00622 // if so.
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             //? Displacement small but in wrong direction. Would probably be ok.
00658         }
00659     }
00660     return true;
00661 }
00662 
00663 
00664 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
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 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
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     // No need to sync; all processors will have all same zonedSurfaces.
00689     if (zonedSurfaces.size())
00690     {
00691         fvMesh& mesh = meshRefiner_.mesh();
00692 
00693         // Split internal faces on interface surfaces
00694         Info<< "Converting zoned faces into baffles ..." << endl;
00695 
00696         // Get faces (internal only) to be baffled. Map from face to patch
00697         // label.
00698         Map<label> faceToPatch(getZoneBafflePatches(false));
00699 
00700         label nZoneFaces = returnReduce(faceToPatch.size(), sumOp<label>());
00701         if (nZoneFaces > 0)
00702         {
00703             // Convert into labelLists
00704             labelList ownPatch(mesh.nFaces(), -1);
00705             forAllConstIter(Map<label>, faceToPatch, iter)
00706             {
00707                 ownPatch[iter.key()] = iter();
00708             }
00709 
00710             // Create baffles. both sides same patch.
00711             map = meshRefiner_.createBaffles(ownPatch, ownPatch);
00712 
00713             // Get pairs of faces created.
00714             // Just loop over faceMap and store baffle if we encounter a slave
00715             // face.
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                 // Does face originate from face-to-patch
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     // No need to sync; all processors will have all same zonedSurfaces.
00774     label nBaffles = returnReduce(baffles.size(), sumOp<label>());
00775     if (zonedSurfaces.size() && nBaffles > 0)
00776     {
00777         // Merge any baffles
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>(),  // combine op
00825         -GREAT,             // null value
00826         false               // no separation
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         // The current mesh is the starting mesh to smooth from.
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             // Try to adapt mesh to obtain displacement by smoothly
00880             // decreasing displacement at error locations.
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     // The current mesh is the starting mesh to smooth from.
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 // Get (pp-local) indices of points that are both on zone and on patched surface
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     // Could use PrimitivePatch & localFaces to extract points but might just
00938     // as well do it ourselves.
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     // Displacement per patch point
00980     vectorField patchDisp(localPoints.size(), vector::zero);
00981 
00982     if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
00983     {
00984         // Current surface snapped to
00985         labelList snapSurf(localPoints.size(), -1);
00986 
00987         // Divide surfaces into zoned and unzoned
00988         labelList zonedSurfaces =
00989             meshRefiner_.surfaces().getNamedSurfaces();
00990         labelList unzonedSurfaces =
00991             meshRefiner_.surfaces().getUnnamedSurfaces();
00992 
00993 
00994         // 1. All points to non-interface surfaces
00995         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00996 
00997         {
00998             List<pointIndexHit> hitInfo;
00999             labelList hitSurface;
01000             surfaces.findNearest
01001             (
01002                 unzonedSurfaces,
01003                 localPoints,
01004                 sqr(4*snapDist),        // sqr of attract distance
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         // 2. All points on zones to their respective surface
01025         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01026 
01027         // Surfaces with zone information
01028         const wordList& faceZoneNames = surfaces.faceZoneNames();
01029 
01030         // Current best snap distance
01031         scalarField minSnapDist(snapDist);
01032 
01033         forAll(zonedSurfaces, i)
01034         {
01035             label zoneSurfI = zonedSurfaces[i];
01036 
01037             const labelList surfacesToTest(1, zoneSurfI);
01038 
01039             // Get indices of points both on faceZone and on pp.
01040             labelList zonePointIndices
01041             (
01042                 getZoneSurfacePoints
01043                 (
01044                     pp,
01045                     faceZoneNames[zoneSurfI]
01046                 )
01047             );
01048 
01049             // Find nearest for points both on faceZone and pp.
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         // Check if all points are being snapped
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     // Limit amount of movement.
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     // Points on zones in one domain but only present as point on other
01126     // will not do condition 2 on all. Sync explicitly.
01127     syncTools::syncPointList
01128     (
01129         mesh,
01130         pp.meshPoints(),
01131         patchDisp,
01132         minMagEqOp(),                   // combine op
01133         vector(GREAT, GREAT, GREAT),    // null value
01134         false                           // no separation
01135     );
01136 
01137 
01138     // Check for displacement being outwards.
01139     outwardsDisplacement(pp, patchDisp);
01140 
01141     // Set initial distribution of displacement field (on patches) from
01142     // patchDisp and make displacement consistent with b.c. on displacement
01143     // pointVectorField.
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     // Set edge diffusivity as inverse of distance to patch
01173     scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL));
01174     //scalarField edgeGamma(mesh.nEdges(), 1.0);
01175     //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
01176 
01177     // Get displacement field
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         // Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
01200         // but this will also delete all pointMesh but not pointFields which
01201         // gives an illegal situation.
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     // Relax displacement until correct mesh
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 // After snapping: correct patching according to nearest surface.
01279 // Code is very similar to calcNearestSurface.
01280 // - calculate face-wise snap distance as max of point-wise
01281 // - calculate face-wise nearest surface point
01282 // - repatch face according to patch for surface point.
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     // Get the labels of added patches.
01295     autoPtr<indirectPrimitivePatch> ppPtr
01296     (
01297         meshRefinement::makePatch
01298         (
01299             mesh,
01300             adaptPatchIDs
01301         )
01302     );
01303     indirectPrimitivePatch& pp = ppPtr();
01304 
01305     // Divide surfaces into zoned and unzoned
01306     labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
01307     labelList unzonedSurfaces = meshRefiner_.surfaces().getUnnamedSurfaces();
01308 
01309 
01310     // Faces that do not move
01311     PackedBoolList isZonedFace(mesh.nFaces(), 0);
01312     {
01313         // 1. All faces on zoned surfaces
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     // Determine per pp face which patch it should be in
01334     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
01335 
01336     // Patch that face should be in
01337     labelList closestPatch(pp.size(), -1);
01338     {
01339         // face snap distance as max of point snap distance
01340         scalarField faceSnapDist(pp.size(), -GREAT);
01341         {
01342             // Distance to attract to nearest feature on surface
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         // Get nearest surface and region
01365         labelList hitSurface;
01366         labelList hitRegion;
01367         surfaces.findNearestRegion
01368         (
01369             unzonedSurfaces,
01370             localFaceCentres,
01371             sqr(4*faceSnapDist),    // sqr of attract distance
01372             hitSurface,
01373             hitRegion
01374         );
01375 
01376         // Get patch
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     // Change those faces for which there is a different closest patch
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     // Get the labels of added patches.
01451     labelList adaptPatchIDs(meshRefiner_.meshedPatches());
01452 
01453     // Create baffles (pairs of faces that share the same points)
01454     // Baffles stored as owner and neighbour face that have been created.
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         // Distance to attract to nearest feature on surface
01470         const scalarField snapDist(calcSnapDistance(snapParams, pp));
01471 
01472 
01473         // Construct iterative mesh mover.
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         // Check initial mesh
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         // Pre-smooth patch vertices (so before determining nearest)
01508         preSmoothPatch(snapParams, nInitErrors, baffles, meshMover);
01509 
01510         // Calculate displacement at every patch point. Insert into
01511         // meshMover.
01512         calcNearestSurface(snapDist, meshMover);
01513 
01515         //- 2009-12-16 : was not found to be beneficial. Keeping internal
01516         // fields fixed slightly increases skewness (on boundary)
01517         // but lowers non-orthogonality quite a bit (e.g. 65->59 degrees).
01518         // Maybe if better smoother?
01519         //smoothDisplacement(snapParams, meshMover);
01520 
01521         // Apply internal displacement to mesh.
01522         scaleMesh(snapParams, nInitErrors, baffles, meshMover);
01523     }
01524 
01525     // Merge any introduced baffles.
01526     mergeZoneBaffles(baffles);
01527 
01528     // Repatch faces according to nearest.
01529     repatchToSurface(snapParams, adaptPatchIDs);
01530 }
01531 
01532 
01533 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines