mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
3071 lines
86 KiB
C
3071 lines
86 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
|
|
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Description
|
|
All to do with snapping to the surface
|
|
|
|
\*----------------------------------------------------------------------------*/
|
|
|
|
#include "snappySnapDriver.H"
|
|
#include "motionSmoother.H"
|
|
#include "polyTopoChange.H"
|
|
#include "syncTools.H"
|
|
#include "fvMesh.H"
|
|
#include "Time.H"
|
|
#include "OFstream.H"
|
|
#include "OBJstream.H"
|
|
#include "mapPolyMesh.H"
|
|
#include "pointEdgePoint.H"
|
|
#include "PointEdgeWave.H"
|
|
#include "mergePoints.H"
|
|
#include "snapParameters.H"
|
|
#include "refinementSurfaces.H"
|
|
#include "searchableSurfaces.H"
|
|
#include "unitConversion.H"
|
|
#include "localPointRegion.H"
|
|
#include "PatchTools.H"
|
|
#include "refinementFeatures.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
namespace Foam
|
|
{
|
|
|
|
defineTypeNameAndDebug(snappySnapDriver, 0);
|
|
|
|
} // End namespace Foam
|
|
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
// Calculate geometrically collocated points, Requires PackedList to be
|
|
// sized and initalised!
|
|
Foam::label Foam::snappySnapDriver::getCollocatedPoints
|
|
(
|
|
const scalar tol,
|
|
const pointField& points,
|
|
PackedBoolList& isCollocatedPoint
|
|
)
|
|
{
|
|
labelList pointMap;
|
|
label nUnique = mergePoints
|
|
(
|
|
points, // points
|
|
tol, // mergeTol
|
|
false, // verbose
|
|
pointMap
|
|
);
|
|
bool hasMerged = (nUnique < points.size());
|
|
|
|
if (!returnReduce(hasMerged, orOp<bool>()))
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
// Determine which merged points are referenced more than once
|
|
label nCollocated = 0;
|
|
|
|
// Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
|
|
// twice)
|
|
labelList firstOldPoint(nUnique, -1);
|
|
forAll(pointMap, oldPointI)
|
|
{
|
|
label newPointI = pointMap[oldPointI];
|
|
|
|
if (firstOldPoint[newPointI] == -1)
|
|
{
|
|
// First use of oldPointI. Store.
|
|
firstOldPoint[newPointI] = oldPointI;
|
|
}
|
|
else if (firstOldPoint[newPointI] == -2)
|
|
{
|
|
// Third or more reference of oldPointI -> non-manifold
|
|
isCollocatedPoint.set(oldPointI, 1u);
|
|
nCollocated++;
|
|
}
|
|
else
|
|
{
|
|
// Second reference of oldPointI -> non-manifold
|
|
isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
|
|
nCollocated++;
|
|
|
|
isCollocatedPoint.set(oldPointI, 1u);
|
|
nCollocated++;
|
|
|
|
// Mark with special value to save checking next time round
|
|
firstOldPoint[newPointI] = -2;
|
|
}
|
|
}
|
|
return returnReduce(nCollocated, sumOp<label>());
|
|
}
|
|
|
|
|
|
Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothInternalDisplacement
|
|
(
|
|
const meshRefinement& meshRefiner,
|
|
const motionSmoother& meshMover
|
|
)
|
|
{
|
|
const indirectPrimitivePatch& pp = meshMover.patch();
|
|
const polyMesh& mesh = meshMover.mesh();
|
|
|
|
// Get neighbour refinement
|
|
const hexRef8& cutter = meshRefiner.meshCutter();
|
|
const labelList& cellLevel = cutter.cellLevel();
|
|
|
|
|
|
// Get the faces on the boundary
|
|
PackedBoolList isFront(mesh.nFaces());
|
|
forAll(pp.addressing(), i)
|
|
{
|
|
isFront[pp.addressing()[i]] = true;
|
|
}
|
|
|
|
// Walk out from the surface a bit. Poor man's FaceCellWave.
|
|
// Commented out for now - not sure if needed and if so how much
|
|
//for (label iter = 0; iter < 2; iter++)
|
|
//{
|
|
// PackedBoolList newIsFront(mesh.nFaces());
|
|
//
|
|
// forAll(isFront, faceI)
|
|
// {
|
|
// if (isFront[faceI])
|
|
// {
|
|
// label own = mesh.faceOwner()[faceI];
|
|
// const cell& ownFaces = mesh.cells()[own];
|
|
// forAll(ownFaces, i)
|
|
// {
|
|
// newIsFront[ownFaces[i]] = true;
|
|
// }
|
|
//
|
|
// if (mesh.isInternalFace(faceI))
|
|
// {
|
|
// label nei = mesh.faceNeighbour()[faceI];
|
|
// const cell& neiFaces = mesh.cells()[nei];
|
|
// forAll(neiFaces, i)
|
|
// {
|
|
// newIsFront[neiFaces[i]] = true;
|
|
// }
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// syncTools::syncFaceList
|
|
// (
|
|
// mesh,
|
|
// newIsFront,
|
|
// orEqOp<unsigned int>()
|
|
// );
|
|
//
|
|
// isFront = newIsFront;
|
|
//}
|
|
|
|
// Mark all points on faces
|
|
// - not on the boundary
|
|
// - inbetween differing refinement levels
|
|
PackedBoolList isMovingPoint(mesh.nPoints());
|
|
|
|
label nInterface = 0;
|
|
|
|
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
|
|
{
|
|
label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
|
|
label neiLevel = cellLevel[mesh.faceNeighbour()[faceI]];
|
|
|
|
if (!isFront[faceI] && ownLevel != neiLevel)
|
|
{
|
|
const face& f = mesh.faces()[faceI];
|
|
forAll(f, fp)
|
|
{
|
|
isMovingPoint[f[fp]] = true;
|
|
}
|
|
|
|
nInterface++;
|
|
}
|
|
}
|
|
|
|
labelList neiCellLevel;
|
|
syncTools::swapBoundaryCellList(mesh, cellLevel, neiCellLevel);
|
|
|
|
for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
|
|
{
|
|
label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
|
|
label neiLevel = neiCellLevel[faceI-mesh.nInternalFaces()];
|
|
|
|
if (!isFront[faceI] && ownLevel != neiLevel)
|
|
{
|
|
const face& f = mesh.faces()[faceI];
|
|
forAll(f, fp)
|
|
{
|
|
isMovingPoint[f[fp]] = true;
|
|
}
|
|
|
|
nInterface++;
|
|
}
|
|
}
|
|
|
|
if (debug)
|
|
{
|
|
reduce(nInterface, sumOp<label>());
|
|
Info<< "Found " << nInterface << " faces out of "
|
|
<< mesh.globalData().nTotalFaces()
|
|
<< " inbetween refinement regions." << endl;
|
|
}
|
|
|
|
// Make sure that points that are coupled to a moving point are marked
|
|
// as well
|
|
syncTools::syncPointList(mesh, isMovingPoint, maxEqOp<unsigned int>(), 0);
|
|
|
|
// Unmark any point on the boundary. If we're doing zero iterations of
|
|
// face-cell wave we might have coupled points not being unmarked.
|
|
forAll(pp.meshPoints(), pointI)
|
|
{
|
|
isMovingPoint[pp.meshPoints()[pointI]] = false;
|
|
}
|
|
|
|
// Make sure that points that are coupled to meshPoints but not on a patch
|
|
// are unmarked as well
|
|
syncTools::syncPointList(mesh, isMovingPoint, minEqOp<unsigned int>(), 1);
|
|
|
|
|
|
// Calculate average of connected cells
|
|
labelList nCells(mesh.nPoints(), 0);
|
|
pointField sumLocation(mesh.nPoints(), Zero);
|
|
|
|
forAll(isMovingPoint, pointI)
|
|
{
|
|
if (isMovingPoint[pointI])
|
|
{
|
|
const labelList& pCells = mesh.pointCells(pointI);
|
|
|
|
forAll(pCells, i)
|
|
{
|
|
sumLocation[pointI] += mesh.cellCentres()[pCells[i]];
|
|
nCells[pointI]++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Sum
|
|
syncTools::syncPointList(mesh, nCells, plusEqOp<label>(), label(0));
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
sumLocation,
|
|
plusEqOp<point>(),
|
|
vector::zero
|
|
);
|
|
|
|
tmp<pointField> tdisplacement(new pointField(mesh.nPoints(), Zero));
|
|
pointField& displacement = tdisplacement.ref();
|
|
|
|
label nAdapted = 0;
|
|
|
|
forAll(displacement, pointI)
|
|
{
|
|
if (nCells[pointI] > 0)
|
|
{
|
|
displacement[pointI] =
|
|
sumLocation[pointI]/nCells[pointI]-mesh.points()[pointI];
|
|
nAdapted++;
|
|
}
|
|
}
|
|
|
|
reduce(nAdapted, sumOp<label>());
|
|
Info<< "Smoothing " << nAdapted << " points inbetween refinement regions."
|
|
<< endl;
|
|
|
|
return tdisplacement;
|
|
}
|
|
|
|
|
|
// Calculate displacement as average of patch points.
|
|
Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
|
|
(
|
|
const motionSmoother& meshMover,
|
|
const List<labelPair>& baffles
|
|
)
|
|
{
|
|
const indirectPrimitivePatch& pp = meshMover.patch();
|
|
|
|
// Calculate geometrically non-manifold points on the patch to be moved.
|
|
PackedBoolList nonManifoldPoint(pp.nPoints());
|
|
label nNonManifoldPoints = getCollocatedPoints
|
|
(
|
|
SMALL,
|
|
pp.localPoints(),
|
|
nonManifoldPoint
|
|
);
|
|
Info<< "Found " << nNonManifoldPoints << " non-manifold point(s)."
|
|
<< endl;
|
|
|
|
|
|
// Average points
|
|
// ~~~~~~~~~~~~~~
|
|
|
|
// We determine three points:
|
|
// - average of (centres of) connected patch faces
|
|
// - average of (centres of) connected internal mesh faces
|
|
// - as fallback: centre of any connected cell
|
|
// so we can do something moderately sensible for non/manifold points.
|
|
|
|
// Note: the averages are calculated properly parallel. This is
|
|
// necessary to get the points shared by processors correct.
|
|
|
|
|
|
const labelListList& pointFaces = pp.pointFaces();
|
|
const labelList& meshPoints = pp.meshPoints();
|
|
const pointField& points = pp.points();
|
|
const polyMesh& mesh = meshMover.mesh();
|
|
|
|
// Get labels of faces to count (master of coupled faces and baffle pairs)
|
|
PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
|
|
|
|
{
|
|
forAll(baffles, i)
|
|
{
|
|
label f0 = baffles[i].first();
|
|
label f1 = baffles[i].second();
|
|
|
|
if (isMasterFace.get(f0))
|
|
{
|
|
// Make f1 a slave
|
|
isMasterFace.unset(f1);
|
|
}
|
|
else if (isMasterFace.get(f1))
|
|
{
|
|
isMasterFace.unset(f0);
|
|
}
|
|
else
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Both sides of baffle consisting of faces " << f0
|
|
<< " and " << f1 << " are already slave faces."
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Get average position of boundary face centres
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
vectorField avgBoundary(pointFaces.size(), Zero);
|
|
labelList nBoundary(pointFaces.size(), 0);
|
|
|
|
forAll(pointFaces, patchPointI)
|
|
{
|
|
const labelList& pFaces = pointFaces[patchPointI];
|
|
|
|
forAll(pFaces, pfI)
|
|
{
|
|
label faceI = pFaces[pfI];
|
|
|
|
if (isMasterFace.get(pp.addressing()[faceI]))
|
|
{
|
|
avgBoundary[patchPointI] += pp[faceI].centre(points);
|
|
nBoundary[patchPointI]++;
|
|
}
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
avgBoundary,
|
|
plusEqOp<point>(), // combine op
|
|
vector::zero // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
nBoundary,
|
|
plusEqOp<label>(), // combine op
|
|
label(0) // null value
|
|
);
|
|
|
|
forAll(avgBoundary, i)
|
|
{
|
|
avgBoundary[i] /= nBoundary[i];
|
|
}
|
|
|
|
|
|
// Get average position of internal face centres
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
vectorField avgInternal;
|
|
labelList nInternal;
|
|
{
|
|
vectorField globalSum(mesh.nPoints(), Zero);
|
|
labelList globalNum(mesh.nPoints(), 0);
|
|
|
|
// Note: no use of pointFaces
|
|
const faceList& faces = mesh.faces();
|
|
|
|
for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
|
|
{
|
|
const face& f = faces[faceI];
|
|
const point& fc = mesh.faceCentres()[faceI];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
globalSum[f[fp]] += fc;
|
|
globalNum[f[fp]]++;
|
|
}
|
|
}
|
|
|
|
// Count coupled faces as internal ones (but only once)
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
|
|
forAll(patches, patchI)
|
|
{
|
|
if
|
|
(
|
|
patches[patchI].coupled()
|
|
&& refCast<const coupledPolyPatch>(patches[patchI]).owner()
|
|
)
|
|
{
|
|
const coupledPolyPatch& pp =
|
|
refCast<const coupledPolyPatch>(patches[patchI]);
|
|
|
|
const vectorField::subField faceCentres = pp.faceCentres();
|
|
|
|
forAll(pp, i)
|
|
{
|
|
const face& f = pp[i];
|
|
const point& fc = faceCentres[i];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
globalSum[f[fp]] += fc;
|
|
globalNum[f[fp]]++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
globalSum,
|
|
plusEqOp<vector>(), // combine op
|
|
vector::zero // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
globalNum,
|
|
plusEqOp<label>(), // combine op
|
|
label(0) // null value
|
|
);
|
|
|
|
avgInternal.setSize(meshPoints.size());
|
|
nInternal.setSize(meshPoints.size());
|
|
|
|
forAll(avgInternal, patchPointI)
|
|
{
|
|
label meshPointI = meshPoints[patchPointI];
|
|
|
|
nInternal[patchPointI] = globalNum[meshPointI];
|
|
|
|
if (nInternal[patchPointI] == 0)
|
|
{
|
|
avgInternal[patchPointI] = globalSum[meshPointI];
|
|
}
|
|
else
|
|
{
|
|
avgInternal[patchPointI] =
|
|
globalSum[meshPointI]
|
|
/ nInternal[patchPointI];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Precalculate any cell using mesh point (replacement of pointCells()[])
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
labelList anyCell(mesh.nPoints(), -1);
|
|
forAll(mesh.faceOwner(), faceI)
|
|
{
|
|
label own = mesh.faceOwner()[faceI];
|
|
const face& f = mesh.faces()[faceI];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
anyCell[f[fp]] = own;
|
|
}
|
|
}
|
|
|
|
|
|
// Displacement to calculate.
|
|
tmp<pointField> tpatchDisp(new pointField(meshPoints.size(), Zero));
|
|
pointField& patchDisp = tpatchDisp.ref();
|
|
|
|
forAll(pointFaces, i)
|
|
{
|
|
label meshPointI = meshPoints[i];
|
|
const point& currentPos = pp.points()[meshPointI];
|
|
|
|
// Now we have the two average points: avgBoundary and avgInternal
|
|
// and how many boundary/internal faces connect to the point
|
|
// (nBoundary, nInternal)
|
|
// Do some blending between the two.
|
|
// Note: the following section has some reasoning behind it but the
|
|
// blending factors can be experimented with.
|
|
|
|
point newPos;
|
|
|
|
if (!nonManifoldPoint.get(i))
|
|
{
|
|
// Points that are manifold. Weight the internal and boundary
|
|
// by their number of faces and blend with
|
|
scalar internalBlend = 0.1;
|
|
scalar blend = 0.1;
|
|
|
|
point avgPos =
|
|
(
|
|
internalBlend*nInternal[i]*avgInternal[i]
|
|
+(1-internalBlend)*nBoundary[i]*avgBoundary[i]
|
|
)
|
|
/ (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
|
|
|
|
newPos = (1-blend)*avgPos + blend*currentPos;
|
|
}
|
|
else if (nInternal[i] == 0)
|
|
{
|
|
// Non-manifold without internal faces. Use any connected cell
|
|
// as internal point instead. Use precalculated any cell to avoid
|
|
// e.g. pointCells()[meshPointI][0]
|
|
|
|
const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
|
|
|
|
scalar cellCBlend = 0.8;
|
|
scalar blend = 0.1;
|
|
|
|
point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
|
|
|
|
newPos = (1-blend)*avgPos + blend*currentPos;
|
|
}
|
|
else
|
|
{
|
|
// Non-manifold point with internal faces connected to them
|
|
scalar internalBlend = 0.9;
|
|
scalar blend = 0.1;
|
|
|
|
point avgPos =
|
|
internalBlend*avgInternal[i]
|
|
+ (1-internalBlend)*avgBoundary[i];
|
|
|
|
newPos = (1-blend)*avgPos + blend*currentPos;
|
|
}
|
|
|
|
patchDisp[i] = newPos - currentPos;
|
|
}
|
|
|
|
return tpatchDisp;
|
|
}
|
|
//XXXXXXX
|
|
//Foam::tmp<Foam::pointField> Foam::snappySnapDriver::avg
|
|
//(
|
|
// const indirectPrimitivePatch& pp,
|
|
// const pointField& localPoints
|
|
//)
|
|
//{
|
|
// const labelListList& pointEdges = pp.pointEdges();
|
|
// const edgeList& edges = pp.edges();
|
|
//
|
|
// tmp<pointField> tavg(new pointField(pointEdges.size(), Zero));
|
|
// pointField& avg = tavg();
|
|
//
|
|
// forAll(pointEdges, vertI)
|
|
// {
|
|
// vector& avgPos = avg[vertI];
|
|
//
|
|
// const labelList& pEdges = pointEdges[vertI];
|
|
//
|
|
// forAll(pEdges, myEdgeI)
|
|
// {
|
|
// const edge& e = edges[pEdges[myEdgeI]];
|
|
//
|
|
// label otherVertI = e.otherVertex(vertI);
|
|
//
|
|
// avgPos += localPoints[otherVertI];
|
|
// }
|
|
//
|
|
// avgPos /= pEdges.size();
|
|
// }
|
|
// return tavg;
|
|
//}
|
|
//Foam::tmp<Foam::pointField>
|
|
//Foam::snappySnapDriver::smoothLambdaMuPatchDisplacement
|
|
//(
|
|
// const motionSmoother& meshMover,
|
|
// const List<labelPair>& baffles
|
|
//)
|
|
//{
|
|
// const indirectPrimitivePatch& pp = meshMover.patch();
|
|
// pointField newLocalPoints(pp.localPoints());
|
|
//
|
|
// const label iters = 90;
|
|
// const scalar lambda = 0.33;
|
|
// const scalar mu = 0.34;
|
|
//
|
|
// for (label iter = 0; iter < iters; iter++)
|
|
// {
|
|
// // Lambda
|
|
// newLocalPoints =
|
|
// (1 - lambda)*newLocalPoints
|
|
// + lambda*avg(pp, newLocalPoints);
|
|
//
|
|
// // Mu
|
|
// newLocalPoints =
|
|
// (1 + mu)*newLocalPoints
|
|
// - mu*avg(pp, newLocalPoints);
|
|
// }
|
|
// return newLocalPoints-pp.localPoints();
|
|
//}
|
|
//XXXXXXX
|
|
|
|
|
|
Foam::tmp<Foam::scalarField> Foam::snappySnapDriver::edgePatchDist
|
|
(
|
|
const pointMesh& pMesh,
|
|
const indirectPrimitivePatch& pp
|
|
)
|
|
{
|
|
const polyMesh& mesh = pMesh();
|
|
|
|
// Set initial changed points to all the patch points
|
|
List<pointEdgePoint> wallInfo(pp.nPoints());
|
|
|
|
forAll(pp.localPoints(), ppI)
|
|
{
|
|
wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
|
|
}
|
|
|
|
// Current info on points
|
|
List<pointEdgePoint> allPointInfo(mesh.nPoints());
|
|
|
|
// Current info on edges
|
|
List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
|
|
|
|
PointEdgeWave<pointEdgePoint> wallCalc
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
wallInfo,
|
|
|
|
allPointInfo,
|
|
allEdgeInfo,
|
|
mesh.globalData().nTotalPoints() // max iterations
|
|
);
|
|
|
|
// Copy edge values into scalarField
|
|
tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
|
|
scalarField& edgeDist = tedgeDist.ref();
|
|
|
|
forAll(allEdgeInfo, edgeI)
|
|
{
|
|
edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
|
|
}
|
|
|
|
return tedgeDist;
|
|
}
|
|
|
|
|
|
void Foam::snappySnapDriver::dumpMove
|
|
(
|
|
const fileName& fName,
|
|
const pointField& meshPts,
|
|
const pointField& surfPts
|
|
)
|
|
{
|
|
// Dump direction of growth into file
|
|
Info<< "Dumping move direction to " << fName << endl;
|
|
|
|
OFstream nearestStream(fName);
|
|
|
|
label vertI = 0;
|
|
|
|
forAll(meshPts, ptI)
|
|
{
|
|
meshTools::writeOBJ(nearestStream, meshPts[ptI]);
|
|
vertI++;
|
|
|
|
meshTools::writeOBJ(nearestStream, surfPts[ptI]);
|
|
vertI++;
|
|
|
|
nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
|
|
}
|
|
}
|
|
|
|
|
|
// Check whether all displacement vectors point outwards of patch. Return true
|
|
// if so.
|
|
bool Foam::snappySnapDriver::outwardsDisplacement
|
|
(
|
|
const indirectPrimitivePatch& pp,
|
|
const vectorField& patchDisp
|
|
)
|
|
{
|
|
const vectorField& faceNormals = pp.faceNormals();
|
|
const labelListList& pointFaces = pp.pointFaces();
|
|
|
|
forAll(pointFaces, pointI)
|
|
{
|
|
const labelList& pFaces = pointFaces[pointI];
|
|
|
|
vector disp(patchDisp[pointI]);
|
|
|
|
scalar magDisp = mag(disp);
|
|
|
|
if (magDisp > SMALL)
|
|
{
|
|
disp /= magDisp;
|
|
|
|
bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
|
|
|
|
if (!outwards)
|
|
{
|
|
Warning<< "Displacement " << patchDisp[pointI]
|
|
<< " at mesh point " << pp.meshPoints()[pointI]
|
|
<< " coord " << pp.points()[pp.meshPoints()[pointI]]
|
|
<< " points through the surrounding patch faces" << endl;
|
|
return false;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
//? Displacement small but in wrong direction. Would probably be ok.
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::snappySnapDriver::snappySnapDriver
|
|
(
|
|
meshRefinement& meshRefiner,
|
|
const labelList& globalToMasterPatch,
|
|
const labelList& globalToSlavePatch
|
|
)
|
|
:
|
|
meshRefiner_(meshRefiner),
|
|
globalToMasterPatch_(globalToMasterPatch),
|
|
globalToSlavePatch_(globalToSlavePatch)
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
Foam::scalarField Foam::snappySnapDriver::calcSnapDistance
|
|
(
|
|
const fvMesh& mesh,
|
|
const snapParameters& snapParams,
|
|
const indirectPrimitivePatch& pp
|
|
)
|
|
{
|
|
const edgeList& edges = pp.edges();
|
|
const labelListList& pointEdges = pp.pointEdges();
|
|
const pointField& localPoints = pp.localPoints();
|
|
|
|
scalarField maxEdgeLen(localPoints.size(), -GREAT);
|
|
|
|
forAll(pointEdges, pointI)
|
|
{
|
|
const labelList& pEdges = pointEdges[pointI];
|
|
|
|
forAll(pEdges, pEdgeI)
|
|
{
|
|
const edge& e = edges[pEdges[pEdgeI]];
|
|
|
|
scalar len = e.mag(localPoints);
|
|
|
|
maxEdgeLen[pointI] = max(maxEdgeLen[pointI], len);
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
maxEdgeLen,
|
|
maxEqOp<scalar>(), // combine op
|
|
-GREAT // null value
|
|
);
|
|
|
|
return scalarField(snapParams.snapTol()*maxEdgeLen);
|
|
}
|
|
|
|
|
|
void Foam::snappySnapDriver::preSmoothPatch
|
|
(
|
|
const meshRefinement& meshRefiner,
|
|
const snapParameters& snapParams,
|
|
const label nInitErrors,
|
|
const List<labelPair>& baffles,
|
|
motionSmoother& meshMover
|
|
)
|
|
{
|
|
const fvMesh& mesh = meshRefiner.mesh();
|
|
|
|
labelList checkFaces;
|
|
|
|
if (snapParams.nSmoothInternal() > 0)
|
|
{
|
|
Info<< "Smoothing patch and internal points ..." << endl;
|
|
}
|
|
else
|
|
{
|
|
Info<< "Smoothing patch points ..." << endl;
|
|
}
|
|
|
|
for
|
|
(
|
|
label smoothIter = 0;
|
|
smoothIter < snapParams.nSmoothPatch();
|
|
smoothIter++
|
|
)
|
|
{
|
|
Info<< "Smoothing iteration " << smoothIter << endl;
|
|
checkFaces.setSize(mesh.nFaces());
|
|
forAll(checkFaces, faceI)
|
|
{
|
|
checkFaces[faceI] = faceI;
|
|
}
|
|
|
|
// If enabled smooth the internal points
|
|
if (snapParams.nSmoothInternal() > smoothIter)
|
|
{
|
|
// Override values on internal points on refinement interfaces
|
|
meshMover.pointDisplacement().internalField() =
|
|
smoothInternalDisplacement(meshRefiner, meshMover);
|
|
}
|
|
|
|
// Smooth the patch points
|
|
pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
|
|
//pointField patchDisp
|
|
//(
|
|
// smoothLambdaMuPatchDisplacement(meshMover, baffles)
|
|
//);
|
|
|
|
// Take over patch displacement as boundary condition on
|
|
// pointDisplacement
|
|
meshMover.setDisplacement(patchDisp);
|
|
|
|
// Start off from current mesh.points()
|
|
meshMover.correct();
|
|
|
|
scalar oldErrorReduction = -1;
|
|
|
|
for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
|
|
{
|
|
Info<< nl << "Scaling iteration " << snapIter << endl;
|
|
|
|
if (snapIter == snapParams.nSnap())
|
|
{
|
|
Info<< "Displacement scaling for error reduction set to 0."
|
|
<< endl;
|
|
oldErrorReduction = meshMover.setErrorReduction(0.0);
|
|
}
|
|
|
|
// Try to adapt mesh to obtain displacement by smoothly
|
|
// decreasing displacement at error locations.
|
|
if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
|
|
{
|
|
Info<< "Successfully moved mesh" << endl;
|
|
break;
|
|
}
|
|
}
|
|
|
|
if (oldErrorReduction >= 0)
|
|
{
|
|
meshMover.setErrorReduction(oldErrorReduction);
|
|
}
|
|
Info<< endl;
|
|
}
|
|
|
|
|
|
// The current mesh is the starting mesh to smooth from.
|
|
meshMover.correct();
|
|
|
|
if (debug&meshRefinement::MESH)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing patch smoothed mesh to time "
|
|
<< meshRefiner.timeName() << '.' << endl;
|
|
meshRefiner.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
mesh.time().path()/meshRefiner.timeName()
|
|
);
|
|
Info<< "Dumped mesh in = "
|
|
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
|
|
}
|
|
|
|
Info<< "Patch points smoothed in = "
|
|
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
|
|
}
|
|
|
|
|
|
// Get (pp-local) indices of points that are both on zone and on patched surface
|
|
Foam::labelList Foam::snappySnapDriver::getZoneSurfacePoints
|
|
(
|
|
const fvMesh& mesh,
|
|
const indirectPrimitivePatch& pp,
|
|
const word& zoneName
|
|
)
|
|
{
|
|
label zoneI = mesh.faceZones().findZoneID(zoneName);
|
|
|
|
if (zoneI == -1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Cannot find zone " << zoneName
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
const faceZone& fZone = mesh.faceZones()[zoneI];
|
|
|
|
|
|
// Could use PrimitivePatch & localFaces to extract points but might just
|
|
// as well do it ourselves.
|
|
|
|
boolList pointOnZone(pp.nPoints(), false);
|
|
|
|
forAll(fZone, i)
|
|
{
|
|
const face& f = mesh.faces()[fZone[i]];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
label meshPointI = f[fp];
|
|
|
|
Map<label>::const_iterator iter =
|
|
pp.meshPointMap().find(meshPointI);
|
|
|
|
if (iter != pp.meshPointMap().end())
|
|
{
|
|
label pointI = iter();
|
|
pointOnZone[pointI] = true;
|
|
}
|
|
}
|
|
}
|
|
|
|
return findIndices(pointOnZone, true);
|
|
}
|
|
|
|
|
|
Foam::tmp<Foam::pointField> Foam::snappySnapDriver::avgCellCentres
|
|
(
|
|
const fvMesh& mesh,
|
|
const indirectPrimitivePatch& pp
|
|
)
|
|
{
|
|
const labelListList& pointFaces = pp.pointFaces();
|
|
|
|
|
|
tmp<pointField> tavgBoundary
|
|
(
|
|
new pointField(pointFaces.size(), Zero)
|
|
);
|
|
pointField& avgBoundary = tavgBoundary.ref();
|
|
labelList nBoundary(pointFaces.size(), 0);
|
|
|
|
forAll(pointFaces, pointI)
|
|
{
|
|
const labelList& pFaces = pointFaces[pointI];
|
|
|
|
forAll(pFaces, pfI)
|
|
{
|
|
label faceI = pFaces[pfI];
|
|
label meshFaceI = pp.addressing()[faceI];
|
|
|
|
label own = mesh.faceOwner()[meshFaceI];
|
|
avgBoundary[pointI] += mesh.cellCentres()[own];
|
|
nBoundary[pointI]++;
|
|
}
|
|
}
|
|
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
avgBoundary,
|
|
plusEqOp<point>(), // combine op
|
|
vector::zero // null value
|
|
);
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
nBoundary,
|
|
plusEqOp<label>(), // combine op
|
|
label(0) // null value
|
|
);
|
|
|
|
forAll(avgBoundary, i)
|
|
{
|
|
avgBoundary[i] /= nBoundary[i];
|
|
}
|
|
return tavgBoundary;
|
|
}
|
|
|
|
|
|
//Foam::tmp<Foam::scalarField> Foam::snappySnapDriver::calcEdgeLen
|
|
//(
|
|
// const indirectPrimitivePatch& pp
|
|
//) const
|
|
//{
|
|
// // Get local edge length based on refinement level
|
|
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
// // (Ripped from snappyLayerDriver)
|
|
//
|
|
// tmp<scalarField> tedgeLen(new scalarField(pp.nPoints()));
|
|
// scalarField& edgeLen = tedgeLen();
|
|
// {
|
|
// const fvMesh& mesh = meshRefiner_.mesh();
|
|
// const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
|
|
// const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
|
|
//
|
|
// labelList maxPointLevel(pp.nPoints(), labelMin);
|
|
//
|
|
// forAll(pp, i)
|
|
// {
|
|
// label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]];
|
|
// const face& f = pp.localFaces()[i];
|
|
// forAll(f, fp)
|
|
// {
|
|
// maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
|
|
// }
|
|
// }
|
|
//
|
|
// syncTools::syncPointList
|
|
// (
|
|
// mesh,
|
|
// pp.meshPoints(),
|
|
// maxPointLevel,
|
|
// maxEqOp<label>(),
|
|
// labelMin // null value
|
|
// );
|
|
//
|
|
//
|
|
// forAll(maxPointLevel, pointI)
|
|
// {
|
|
// // Find undistorted edge size for this level.
|
|
// edgeLen[pointI] = edge0Len/(1<<maxPointLevel[pointI]);
|
|
// }
|
|
// }
|
|
// return tedgeLen;
|
|
//}
|
|
|
|
|
|
void Foam::snappySnapDriver::detectNearSurfaces
|
|
(
|
|
const scalar planarCos,
|
|
const indirectPrimitivePatch& pp,
|
|
const pointField& nearestPoint,
|
|
const vectorField& nearestNormal,
|
|
|
|
vectorField& disp
|
|
) const
|
|
{
|
|
Info<< "Detecting near surfaces ..." << endl;
|
|
|
|
const pointField& localPoints = pp.localPoints();
|
|
const labelList& meshPoints = pp.meshPoints();
|
|
const refinementSurfaces& surfaces = meshRefiner_.surfaces();
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
//// Get local edge length based on refinement level
|
|
//const scalarField edgeLen(calcEdgeLen(pp));
|
|
//
|
|
//// Generate rays for every surface point
|
|
//// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
//
|
|
//{
|
|
// const scalar cos45 = Foam::cos(degToRad(45));
|
|
// vector n(cos45, cos45, cos45);
|
|
// n /= mag(n);
|
|
//
|
|
// pointField start(14*pp.nPoints());
|
|
// pointField end(start.size());
|
|
//
|
|
// label rayI = 0;
|
|
// forAll(localPoints, pointI)
|
|
// {
|
|
// const point& pt = localPoints[pointI];
|
|
//
|
|
// // Along coordinate axes
|
|
//
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() -= edgeLen[pointI];
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() += edgeLen[pointI];
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.y() -= edgeLen[pointI];
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.y() += edgeLen[pointI];
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.z() -= edgeLen[pointI];
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.z() += edgeLen[pointI];
|
|
// }
|
|
//
|
|
// // At 45 degrees
|
|
//
|
|
// const vector vec(edgeLen[pointI]*n);
|
|
//
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() += vec.x();
|
|
// endPt.y() += vec.y();
|
|
// endPt.z() += vec.z();
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() -= vec.x();
|
|
// endPt.y() += vec.y();
|
|
// endPt.z() += vec.z();
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() += vec.x();
|
|
// endPt.y() -= vec.y();
|
|
// endPt.z() += vec.z();
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() -= vec.x();
|
|
// endPt.y() -= vec.y();
|
|
// endPt.z() += vec.z();
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() += vec.x();
|
|
// endPt.y() += vec.y();
|
|
// endPt.z() -= vec.z();
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() -= vec.x();
|
|
// endPt.y() += vec.y();
|
|
// endPt.z() -= vec.z();
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() += vec.x();
|
|
// endPt.y() -= vec.y();
|
|
// endPt.z() -= vec.z();
|
|
// }
|
|
// {
|
|
// start[rayI] = pt;
|
|
// point& endPt = end[rayI++];
|
|
// endPt = pt;
|
|
// endPt.x() -= vec.x();
|
|
// endPt.y() -= vec.y();
|
|
// endPt.z() -= vec.z();
|
|
// }
|
|
// }
|
|
//
|
|
// labelList surface1;
|
|
// List<pointIndexHit> hit1;
|
|
// labelList region1;
|
|
// vectorField normal1;
|
|
//
|
|
// labelList surface2;
|
|
// List<pointIndexHit> hit2;
|
|
// labelList region2;
|
|
// vectorField normal2;
|
|
// surfaces.findNearestIntersection
|
|
// (
|
|
// unzonedSurfaces, // surfacesToTest,
|
|
// start,
|
|
// end,
|
|
//
|
|
// surface1,
|
|
// hit1,
|
|
// region1,
|
|
// normal1,
|
|
//
|
|
// surface2,
|
|
// hit2,
|
|
// region2,
|
|
// normal2
|
|
// );
|
|
//
|
|
// // All intersections
|
|
// {
|
|
// OBJstream str
|
|
// (
|
|
// mesh.time().path()
|
|
// / "surfaceHits_" + meshRefiner_.timeName() + ".obj"
|
|
// );
|
|
//
|
|
// Info<< "Dumping intersections with rays to " << str.name()
|
|
// << endl;
|
|
//
|
|
// forAll(hit1, i)
|
|
// {
|
|
// if (hit1[i].hit())
|
|
// {
|
|
// str.write(linePointRef(start[i], hit1[i].hitPoint()));
|
|
// }
|
|
// if (hit2[i].hit())
|
|
// {
|
|
// str.write(linePointRef(start[i], hit2[i].hitPoint()));
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// // Co-planar intersections
|
|
// {
|
|
// OBJstream str
|
|
// (
|
|
// mesh.time().path()
|
|
// / "coplanarHits_" + meshRefiner_.timeName() + ".obj"
|
|
// );
|
|
//
|
|
// Info<< "Dumping intersections with co-planar surfaces to "
|
|
// << str.name() << endl;
|
|
//
|
|
// forAll(localPoints, pointI)
|
|
// {
|
|
// bool hasNormal = false;
|
|
// point surfPointA;
|
|
// vector surfNormalA;
|
|
// point surfPointB;
|
|
// vector surfNormalB;
|
|
//
|
|
// bool isCoplanar = false;
|
|
//
|
|
// label rayI = 14*pointI;
|
|
// for (label i = 0; i < 14; i++)
|
|
// {
|
|
// if (hit1[rayI].hit())
|
|
// {
|
|
// const point& pt = hit1[rayI].hitPoint();
|
|
// const vector& n = normal1[rayI];
|
|
//
|
|
// if (!hasNormal)
|
|
// {
|
|
// hasNormal = true;
|
|
// surfPointA = pt;
|
|
// surfNormalA = n;
|
|
// }
|
|
// else
|
|
// {
|
|
// if
|
|
// (
|
|
// meshRefiner_.isGap
|
|
// (
|
|
// planarCos,
|
|
// surfPointA,
|
|
// surfNormalA,
|
|
// pt,
|
|
// n
|
|
// )
|
|
// )
|
|
// {
|
|
// isCoplanar = true;
|
|
// surfPointB = pt;
|
|
// surfNormalB = n;
|
|
// break;
|
|
// }
|
|
// }
|
|
// }
|
|
// if (hit2[rayI].hit())
|
|
// {
|
|
// const point& pt = hit2[rayI].hitPoint();
|
|
// const vector& n = normal2[rayI];
|
|
//
|
|
// if (!hasNormal)
|
|
// {
|
|
// hasNormal = true;
|
|
// surfPointA = pt;
|
|
// surfNormalA = n;
|
|
// }
|
|
// else
|
|
// {
|
|
// if
|
|
// (
|
|
// meshRefiner_.isGap
|
|
// (
|
|
// planarCos,
|
|
// surfPointA,
|
|
// surfNormalA,
|
|
// pt,
|
|
// n
|
|
// )
|
|
// )
|
|
// {
|
|
// isCoplanar = true;
|
|
// surfPointB = pt;
|
|
// surfNormalB = n;
|
|
// break;
|
|
// }
|
|
// }
|
|
// }
|
|
//
|
|
// rayI++;
|
|
// }
|
|
//
|
|
// if (isCoplanar)
|
|
// {
|
|
// str.write(linePointRef(surfPointA, surfPointB));
|
|
// }
|
|
// }
|
|
// }
|
|
//}
|
|
|
|
|
|
const pointField avgCc(avgCellCentres(mesh, pp));
|
|
|
|
// Construct rays through localPoints to beyond cell centre
|
|
pointField start(pp.nPoints());
|
|
pointField end(pp.nPoints());
|
|
forAll(localPoints, pointI)
|
|
{
|
|
const point& pt = localPoints[pointI];
|
|
const vector d = 2*(avgCc[pointI]-pt);
|
|
start[pointI] = pt - d;
|
|
end[pointI] = pt + d;
|
|
}
|
|
|
|
|
|
autoPtr<OBJstream> gapStr;
|
|
if (debug&meshRefinement::ATTRACTION)
|
|
{
|
|
gapStr.reset
|
|
(
|
|
new OBJstream
|
|
(
|
|
mesh.time().path()
|
|
/ "detectNearSurfaces_" + meshRefiner_.timeName() + ".obj"
|
|
)
|
|
);
|
|
}
|
|
|
|
|
|
const PackedBoolList isPatchMasterPoint
|
|
(
|
|
meshRefinement::getMasterPoints
|
|
(
|
|
mesh,
|
|
meshPoints
|
|
)
|
|
);
|
|
|
|
label nOverride = 0;
|
|
|
|
// 1. All points to non-interface surfaces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
{
|
|
const labelList unzonedSurfaces =
|
|
surfaceZonesInfo::getUnnamedSurfaces
|
|
(
|
|
meshRefiner_.surfaces().surfZones()
|
|
);
|
|
|
|
// Do intersection test
|
|
labelList surface1;
|
|
List<pointIndexHit> hit1;
|
|
labelList region1;
|
|
vectorField normal1;
|
|
|
|
labelList surface2;
|
|
List<pointIndexHit> hit2;
|
|
labelList region2;
|
|
vectorField normal2;
|
|
surfaces.findNearestIntersection
|
|
(
|
|
unzonedSurfaces,
|
|
start,
|
|
end,
|
|
|
|
surface1,
|
|
hit1,
|
|
region1,
|
|
normal1,
|
|
|
|
surface2,
|
|
hit2,
|
|
region2,
|
|
normal2
|
|
);
|
|
|
|
|
|
forAll(localPoints, pointI)
|
|
{
|
|
// Current location
|
|
const point& pt = localPoints[pointI];
|
|
|
|
bool override = false;
|
|
|
|
//if (hit1[pointI].hit())
|
|
//{
|
|
// if
|
|
// (
|
|
// meshRefiner_.isGap
|
|
// (
|
|
// planarCos,
|
|
// nearestPoint[pointI],
|
|
// nearestNormal[pointI],
|
|
// hit1[pointI].hitPoint(),
|
|
// normal1[pointI]
|
|
// )
|
|
// )
|
|
// {
|
|
// disp[pointI] = hit1[pointI].hitPoint()-pt;
|
|
// override = true;
|
|
// }
|
|
//}
|
|
//if (hit2[pointI].hit())
|
|
//{
|
|
// if
|
|
// (
|
|
// meshRefiner_.isGap
|
|
// (
|
|
// planarCos,
|
|
// nearestPoint[pointI],
|
|
// nearestNormal[pointI],
|
|
// hit2[pointI].hitPoint(),
|
|
// normal2[pointI]
|
|
// )
|
|
// )
|
|
// {
|
|
// disp[pointI] = hit2[pointI].hitPoint()-pt;
|
|
// override = true;
|
|
// }
|
|
//}
|
|
|
|
if (hit1[pointI].hit() && hit2[pointI].hit())
|
|
{
|
|
if
|
|
(
|
|
meshRefiner_.isGap
|
|
(
|
|
planarCos,
|
|
hit1[pointI].hitPoint(),
|
|
normal1[pointI],
|
|
hit2[pointI].hitPoint(),
|
|
normal2[pointI]
|
|
)
|
|
)
|
|
{
|
|
// TBD: check if the attraction (to nearest) would attract
|
|
// good enough and not override attraction
|
|
|
|
if (gapStr.valid())
|
|
{
|
|
const point& intPt = hit2[pointI].hitPoint();
|
|
gapStr().write(linePointRef(pt, intPt));
|
|
}
|
|
|
|
// Choose hit2 : nearest to end point (so inside the domain)
|
|
disp[pointI] = hit2[pointI].hitPoint()-pt;
|
|
override = true;
|
|
}
|
|
}
|
|
|
|
if (override && isPatchMasterPoint[pointI])
|
|
{
|
|
nOverride++;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// 2. All points on zones to their respective surface
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
{
|
|
// Surfaces with zone information
|
|
const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
|
|
|
|
const labelList zonedSurfaces = surfaceZonesInfo::getNamedSurfaces
|
|
(
|
|
surfZones
|
|
);
|
|
|
|
forAll(zonedSurfaces, i)
|
|
{
|
|
label zoneSurfI = zonedSurfaces[i];
|
|
|
|
const word& faceZoneName = surfZones[zoneSurfI].faceZoneName();
|
|
|
|
const labelList surfacesToTest(1, zoneSurfI);
|
|
|
|
// Get indices of points both on faceZone and on pp.
|
|
labelList zonePointIndices
|
|
(
|
|
getZoneSurfacePoints
|
|
(
|
|
mesh,
|
|
pp,
|
|
faceZoneName
|
|
)
|
|
);
|
|
|
|
// Do intersection test
|
|
labelList surface1;
|
|
List<pointIndexHit> hit1;
|
|
labelList region1;
|
|
vectorField normal1;
|
|
|
|
labelList surface2;
|
|
List<pointIndexHit> hit2;
|
|
labelList region2;
|
|
vectorField normal2;
|
|
surfaces.findNearestIntersection
|
|
(
|
|
surfacesToTest,
|
|
pointField(start, zonePointIndices),
|
|
pointField(end, zonePointIndices),
|
|
|
|
surface1,
|
|
hit1,
|
|
region1,
|
|
normal1,
|
|
|
|
surface2,
|
|
hit2,
|
|
region2,
|
|
normal2
|
|
);
|
|
|
|
|
|
forAll(hit1, i)
|
|
{
|
|
label pointI = zonePointIndices[i];
|
|
|
|
// Current location
|
|
const point& pt = localPoints[pointI];
|
|
|
|
bool override = false;
|
|
|
|
//if (hit1[i].hit())
|
|
//{
|
|
// if
|
|
// (
|
|
// meshRefiner_.isGap
|
|
// (
|
|
// planarCos,
|
|
// nearestPoint[pointI],
|
|
// nearestNormal[pointI],
|
|
// hit1[i].hitPoint(),
|
|
// normal1[i]
|
|
// )
|
|
// )
|
|
// {
|
|
// disp[pointI] = hit1[i].hitPoint()-pt;
|
|
// override = true;
|
|
// }
|
|
//}
|
|
//if (hit2[i].hit())
|
|
//{
|
|
// if
|
|
// (
|
|
// meshRefiner_.isGap
|
|
// (
|
|
// planarCos,
|
|
// nearestPoint[pointI],
|
|
// nearestNormal[pointI],
|
|
// hit2[i].hitPoint(),
|
|
// normal2[i]
|
|
// )
|
|
// )
|
|
// {
|
|
// disp[pointI] = hit2[i].hitPoint()-pt;
|
|
// override = true;
|
|
// }
|
|
//}
|
|
|
|
if (hit1[i].hit() && hit2[i].hit())
|
|
{
|
|
if
|
|
(
|
|
meshRefiner_.isGap
|
|
(
|
|
planarCos,
|
|
hit1[i].hitPoint(),
|
|
normal1[i],
|
|
hit2[i].hitPoint(),
|
|
normal2[i]
|
|
)
|
|
)
|
|
{
|
|
if (gapStr.valid())
|
|
{
|
|
const point& intPt = hit2[i].hitPoint();
|
|
gapStr().write(linePointRef(pt, intPt));
|
|
}
|
|
|
|
disp[pointI] = hit2[i].hitPoint()-pt;
|
|
override = true;
|
|
}
|
|
}
|
|
|
|
if (override && isPatchMasterPoint[pointI])
|
|
{
|
|
nOverride++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< "Overriding nearest with intersection of close gaps at "
|
|
<< returnReduce(nOverride, sumOp<label>())
|
|
<< " out of " << returnReduce(pp.nPoints(), sumOp<label>())
|
|
<< " points." << endl;
|
|
}
|
|
|
|
|
|
void Foam::snappySnapDriver::calcNearestSurface
|
|
(
|
|
const refinementSurfaces& surfaces,
|
|
|
|
const labelList& surfacesToTest,
|
|
const labelListList& regionsToTest,
|
|
|
|
const pointField& localPoints,
|
|
const labelList& zonePointIndices,
|
|
|
|
scalarField& minSnapDist,
|
|
labelList& snapSurf,
|
|
vectorField& patchDisp,
|
|
|
|
// Optional: nearest point, normal
|
|
pointField& nearestPoint,
|
|
vectorField& nearestNormal
|
|
)
|
|
{
|
|
// Find nearest for points both on faceZone and pp.
|
|
List<pointIndexHit> hitInfo;
|
|
labelList hitSurface;
|
|
|
|
if (nearestNormal.size() == localPoints.size())
|
|
{
|
|
labelList hitRegion;
|
|
vectorField hitNormal;
|
|
surfaces.findNearestRegion
|
|
(
|
|
surfacesToTest,
|
|
regionsToTest,
|
|
|
|
pointField(localPoints, zonePointIndices),
|
|
sqr(scalarField(minSnapDist, zonePointIndices)),
|
|
|
|
hitSurface,
|
|
hitInfo,
|
|
hitRegion,
|
|
hitNormal
|
|
);
|
|
|
|
forAll(hitInfo, i)
|
|
{
|
|
if (hitInfo[i].hit())
|
|
{
|
|
label pointI = zonePointIndices[i];
|
|
nearestPoint[pointI] = hitInfo[i].hitPoint();
|
|
nearestNormal[pointI] = hitNormal[i];
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
surfaces.findNearest
|
|
(
|
|
surfacesToTest,
|
|
regionsToTest,
|
|
|
|
pointField(localPoints, zonePointIndices),
|
|
sqr(scalarField(minSnapDist, zonePointIndices)),
|
|
|
|
hitSurface,
|
|
hitInfo
|
|
);
|
|
}
|
|
|
|
forAll(hitInfo, i)
|
|
{
|
|
if (hitInfo[i].hit())
|
|
{
|
|
label pointI = zonePointIndices[i];
|
|
|
|
patchDisp[pointI] = hitInfo[i].hitPoint() - localPoints[pointI];
|
|
minSnapDist[pointI] = mag(patchDisp[pointI]);
|
|
snapSurf[pointI] = hitSurface[i];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::vectorField Foam::snappySnapDriver::calcNearestSurface
|
|
(
|
|
const bool strictRegionSnap,
|
|
const meshRefinement& meshRefiner,
|
|
const labelList& globalToMasterPatch,
|
|
const labelList& globalToSlavePatch,
|
|
const scalarField& snapDist,
|
|
const indirectPrimitivePatch& pp,
|
|
pointField& nearestPoint,
|
|
vectorField& nearestNormal
|
|
)
|
|
{
|
|
Info<< "Calculating patchDisplacement as distance to nearest surface"
|
|
<< " point ..." << endl;
|
|
if (strictRegionSnap)
|
|
{
|
|
Info<< " non-zone points : attract to local region on surface only"
|
|
<< nl
|
|
<< " zone points : attract to local region on surface only"
|
|
<< nl
|
|
<< endl;
|
|
}
|
|
else
|
|
{
|
|
Info<< " non-zone points :"
|
|
<< " attract to nearest of all non-zone surfaces"
|
|
<< nl
|
|
<< " zone points : attract to zone surface only" << nl
|
|
<< endl;
|
|
}
|
|
|
|
|
|
const pointField& localPoints = pp.localPoints();
|
|
const refinementSurfaces& surfaces = meshRefiner.surfaces();
|
|
const fvMesh& mesh = meshRefiner.mesh();
|
|
|
|
// Displacement per patch point
|
|
vectorField patchDisp(localPoints.size(), Zero);
|
|
|
|
if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
|
|
{
|
|
// Current surface snapped to. Used to check whether points have been
|
|
// snapped at all
|
|
labelList snapSurf(localPoints.size(), -1);
|
|
|
|
// Current best snap distance (since point might be on multiple
|
|
// regions)
|
|
scalarField minSnapDist(snapDist);
|
|
|
|
|
|
if (strictRegionSnap)
|
|
{
|
|
// Attract patch points to same region only
|
|
|
|
forAll(surfaces.surfaces(), surfI)
|
|
{
|
|
label geomI = surfaces.surfaces()[surfI];
|
|
label nRegions = surfaces.geometry()[geomI].regions().size();
|
|
|
|
const labelList surfacesToTest(1, surfI);
|
|
|
|
for (label regionI = 0; regionI < nRegions; regionI++)
|
|
{
|
|
label globalI = surfaces.globalRegion(surfI, regionI);
|
|
label masterPatchI = globalToMasterPatch[globalI];
|
|
|
|
// Get indices of points both on patch and on pp
|
|
labelList zonePointIndices
|
|
(
|
|
getFacePoints
|
|
(
|
|
pp,
|
|
mesh.boundaryMesh()[masterPatchI]
|
|
)
|
|
);
|
|
|
|
calcNearestSurface
|
|
(
|
|
surfaces,
|
|
|
|
surfacesToTest,
|
|
labelListList(1, labelList(1, regionI)), //regionsToTest
|
|
|
|
localPoints,
|
|
zonePointIndices,
|
|
|
|
minSnapDist,
|
|
snapSurf,
|
|
patchDisp,
|
|
|
|
// Optional: nearest point, normal
|
|
nearestPoint,
|
|
nearestNormal
|
|
);
|
|
|
|
if (globalToSlavePatch[globalI] != masterPatchI)
|
|
{
|
|
label slavePatchI = globalToSlavePatch[globalI];
|
|
|
|
// Get indices of points both on patch and on pp
|
|
labelList zonePointIndices
|
|
(
|
|
getFacePoints
|
|
(
|
|
pp,
|
|
mesh.boundaryMesh()[slavePatchI]
|
|
)
|
|
);
|
|
|
|
calcNearestSurface
|
|
(
|
|
surfaces,
|
|
|
|
surfacesToTest,
|
|
labelListList(1, labelList(1, regionI)),
|
|
|
|
localPoints,
|
|
zonePointIndices,
|
|
|
|
minSnapDist,
|
|
snapSurf,
|
|
patchDisp,
|
|
|
|
// Optional: nearest point, normal
|
|
nearestPoint,
|
|
nearestNormal
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Divide surfaces into zoned and unzoned
|
|
const labelList unzonedSurfaces =
|
|
surfaceZonesInfo::getUnnamedSurfaces
|
|
(
|
|
meshRefiner.surfaces().surfZones()
|
|
);
|
|
|
|
|
|
// 1. All points to non-interface surfaces
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
List<pointIndexHit> hitInfo;
|
|
labelList hitSurface;
|
|
|
|
if (nearestNormal.size() == localPoints.size())
|
|
{
|
|
labelList hitRegion;
|
|
vectorField hitNormal;
|
|
surfaces.findNearestRegion
|
|
(
|
|
unzonedSurfaces,
|
|
localPoints,
|
|
sqr(snapDist),
|
|
hitSurface,
|
|
hitInfo,
|
|
hitRegion,
|
|
hitNormal
|
|
);
|
|
|
|
forAll(hitInfo, pointI)
|
|
{
|
|
if (hitInfo[pointI].hit())
|
|
{
|
|
nearestPoint[pointI] = hitInfo[pointI].hitPoint();
|
|
nearestNormal[pointI] = hitNormal[pointI];
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
surfaces.findNearest
|
|
(
|
|
unzonedSurfaces,
|
|
localPoints,
|
|
sqr(snapDist), // sqr of attract distance
|
|
hitSurface,
|
|
hitInfo
|
|
);
|
|
}
|
|
|
|
forAll(hitInfo, pointI)
|
|
{
|
|
if (hitInfo[pointI].hit())
|
|
{
|
|
patchDisp[pointI] =
|
|
hitInfo[pointI].hitPoint()
|
|
- localPoints[pointI];
|
|
|
|
snapSurf[pointI] = hitSurface[pointI];
|
|
}
|
|
}
|
|
|
|
|
|
const labelList zonedSurfaces =
|
|
surfaceZonesInfo::getNamedSurfaces
|
|
(
|
|
meshRefiner.surfaces().surfZones()
|
|
);
|
|
|
|
|
|
// 2. All points on zones to their respective surface
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Surfaces with zone information
|
|
const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
|
|
|
|
forAll(zonedSurfaces, i)
|
|
{
|
|
label surfI = zonedSurfaces[i];
|
|
|
|
const word& faceZoneName = surfZones[surfI].faceZoneName();
|
|
|
|
const labelList surfacesToTest(1, surfI);
|
|
|
|
label geomI = surfaces.surfaces()[surfI];
|
|
label nRegions = surfaces.geometry()[geomI].regions().size();
|
|
|
|
|
|
// Get indices of points both on faceZone and on pp.
|
|
labelList zonePointIndices
|
|
(
|
|
getZoneSurfacePoints
|
|
(
|
|
mesh,
|
|
pp,
|
|
faceZoneName
|
|
)
|
|
);
|
|
|
|
|
|
calcNearestSurface
|
|
(
|
|
surfaces,
|
|
|
|
surfacesToTest,
|
|
labelListList(1, identity(nRegions)),
|
|
|
|
localPoints,
|
|
zonePointIndices,
|
|
|
|
minSnapDist,
|
|
snapSurf,
|
|
patchDisp,
|
|
|
|
// Optional: nearest point, normal
|
|
nearestPoint,
|
|
nearestNormal
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// Check if all points are being snapped
|
|
forAll(snapSurf, pointI)
|
|
{
|
|
if (snapSurf[pointI] == -1)
|
|
{
|
|
WarningInFunction
|
|
<< "For point:" << pointI
|
|
<< " coordinate:" << localPoints[pointI]
|
|
<< " did not find any surface within:"
|
|
<< minSnapDist[pointI]
|
|
<< " metre." << endl;
|
|
}
|
|
}
|
|
|
|
{
|
|
const PackedBoolList isPatchMasterPoint
|
|
(
|
|
meshRefinement::getMasterPoints
|
|
(
|
|
mesh,
|
|
pp.meshPoints()
|
|
)
|
|
);
|
|
|
|
scalarField magDisp(mag(patchDisp));
|
|
|
|
Info<< "Wanted displacement : average:"
|
|
<< meshRefinement::gAverage(isPatchMasterPoint, magDisp)
|
|
<< " min:" << gMin(magDisp)
|
|
<< " max:" << gMax(magDisp) << endl;
|
|
}
|
|
}
|
|
|
|
Info<< "Calculated surface displacement in = "
|
|
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
|
|
|
|
|
|
// Limit amount of movement. Can not happen for triSurfaceMesh but
|
|
// can happen for some analytical shapes?
|
|
forAll(patchDisp, patchPointI)
|
|
{
|
|
scalar magDisp = mag(patchDisp[patchPointI]);
|
|
|
|
if (magDisp > snapDist[patchPointI])
|
|
{
|
|
patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
|
|
|
|
Pout<< "Limiting displacement for " << patchPointI
|
|
<< " from " << magDisp << " to " << snapDist[patchPointI]
|
|
<< endl;
|
|
}
|
|
}
|
|
|
|
// Points on zones in one domain but only present as point on other
|
|
// will not do condition 2 on all. Sync explicitly.
|
|
syncTools::syncPointList
|
|
(
|
|
mesh,
|
|
pp.meshPoints(),
|
|
patchDisp,
|
|
minMagSqrEqOp<point>(), // combine op
|
|
vector(GREAT, GREAT, GREAT) // null value (note: cannot use VGREAT)
|
|
);
|
|
|
|
return patchDisp;
|
|
}
|
|
|
|
|
|
void Foam::snappySnapDriver::smoothDisplacement
|
|
(
|
|
const snapParameters& snapParams,
|
|
motionSmoother& meshMover
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
const indirectPrimitivePatch& pp = meshMover.patch();
|
|
|
|
Info<< "Smoothing displacement ..." << endl;
|
|
|
|
// Set edge diffusivity as inverse of distance to patch
|
|
scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + SMALL));
|
|
//scalarField edgeGamma(mesh.nEdges(), 1.0);
|
|
//scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
|
|
|
|
// Get displacement field
|
|
pointVectorField& disp = meshMover.displacement();
|
|
|
|
for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
|
|
{
|
|
if ((iter % 10) == 0)
|
|
{
|
|
Info<< "Iteration " << iter << endl;
|
|
}
|
|
pointVectorField oldDisp(disp);
|
|
meshMover.smooth(oldDisp, edgeGamma, disp);
|
|
}
|
|
Info<< "Displacement smoothed in = "
|
|
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
|
|
|
|
if (debug&meshRefinement::MESH)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
|
|
<< endl;
|
|
|
|
// Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
|
|
// but this will also delete all pointMesh but not pointFields which
|
|
// gives an illegal situation.
|
|
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
mesh.time().path()/meshRefiner_.timeName()
|
|
);
|
|
Info<< "Writing displacement field ..." << endl;
|
|
disp.write();
|
|
tmp<pointScalarField> magDisp(mag(disp));
|
|
magDisp().write();
|
|
|
|
Info<< "Writing actual patch displacement ..." << endl;
|
|
vectorField actualPatchDisp(disp, pp.meshPoints());
|
|
dumpMove
|
|
(
|
|
mesh.time().path()
|
|
/ "actualPatchDisplacement_" + meshRefiner_.timeName() + ".obj",
|
|
pp.localPoints(),
|
|
pp.localPoints() + actualPatchDisp
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
bool Foam::snappySnapDriver::scaleMesh
|
|
(
|
|
const snapParameters& snapParams,
|
|
const label nInitErrors,
|
|
const List<labelPair>& baffles,
|
|
motionSmoother& meshMover
|
|
)
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
// Relax displacement until correct mesh
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
labelList checkFaces(identity(mesh.nFaces()));
|
|
|
|
scalar oldErrorReduction = -1;
|
|
|
|
bool meshOk = false;
|
|
|
|
Info<< "Moving mesh ..." << endl;
|
|
for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
|
|
{
|
|
Info<< nl << "Iteration " << iter << endl;
|
|
|
|
if (iter == snapParams.nSnap())
|
|
{
|
|
Info<< "Displacement scaling for error reduction set to 0." << endl;
|
|
oldErrorReduction = meshMover.setErrorReduction(0.0);
|
|
}
|
|
|
|
meshOk = meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors);
|
|
|
|
if (meshOk)
|
|
{
|
|
Info<< "Successfully moved mesh" << endl;
|
|
break;
|
|
}
|
|
if (debug&meshRefinement::MESH)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing scaled mesh to time " << meshRefiner_.timeName()
|
|
<< endl;
|
|
mesh.write();
|
|
|
|
Info<< "Writing displacement field ..." << endl;
|
|
meshMover.displacement().write();
|
|
tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
|
|
magDisp().write();
|
|
}
|
|
}
|
|
|
|
if (oldErrorReduction >= 0)
|
|
{
|
|
meshMover.setErrorReduction(oldErrorReduction);
|
|
}
|
|
Info<< "Moved mesh in = "
|
|
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
|
|
|
|
return meshOk;
|
|
}
|
|
|
|
|
|
// After snapping: correct patching according to nearest surface.
|
|
// Code is very similar to calcNearestSurface.
|
|
// - calculate face-wise snap distance as max of point-wise
|
|
// - calculate face-wise nearest surface point
|
|
// - repatch face according to patch for surface point.
|
|
Foam::autoPtr<Foam::mapPolyMesh> Foam::snappySnapDriver::repatchToSurface
|
|
(
|
|
const snapParameters& snapParams,
|
|
const labelList& adaptPatchIDs,
|
|
const labelList& preserveFaces
|
|
)
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
const refinementSurfaces& surfaces = meshRefiner_.surfaces();
|
|
|
|
Info<< "Repatching faces according to nearest surface ..." << endl;
|
|
|
|
// Get the labels of added patches.
|
|
autoPtr<indirectPrimitivePatch> ppPtr
|
|
(
|
|
meshRefinement::makePatch
|
|
(
|
|
mesh,
|
|
adaptPatchIDs
|
|
)
|
|
);
|
|
indirectPrimitivePatch& pp = ppPtr();
|
|
|
|
// Divide surfaces into zoned and unzoned
|
|
labelList zonedSurfaces =
|
|
surfaceZonesInfo::getNamedSurfaces(surfaces.surfZones());
|
|
labelList unzonedSurfaces =
|
|
surfaceZonesInfo::getUnnamedSurfaces(surfaces.surfZones());
|
|
|
|
|
|
// Faces that do not move
|
|
PackedBoolList isZonedFace(mesh.nFaces());
|
|
{
|
|
// 1. Preserve faces in preserveFaces list
|
|
forAll(preserveFaces, faceI)
|
|
{
|
|
if (preserveFaces[faceI] != -1)
|
|
{
|
|
isZonedFace.set(faceI, 1);
|
|
}
|
|
}
|
|
|
|
// 2. All faces on zoned surfaces
|
|
const PtrList<surfaceZonesInfo>& surfZones = surfaces.surfZones();
|
|
const faceZoneMesh& fZones = mesh.faceZones();
|
|
|
|
forAll(zonedSurfaces, i)
|
|
{
|
|
const label zoneSurfI = zonedSurfaces[i];
|
|
const faceZone& fZone = fZones[surfZones[zoneSurfI].faceZoneName()];
|
|
|
|
forAll(fZone, i)
|
|
{
|
|
isZonedFace.set(fZone[i], 1);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Determine per pp face which patch it should be in
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// Patch that face should be in
|
|
labelList closestPatch(pp.size(), -1);
|
|
{
|
|
// face snap distance as max of point snap distance
|
|
scalarField faceSnapDist(pp.size(), -GREAT);
|
|
{
|
|
// Distance to attract to nearest feature on surface
|
|
const scalarField snapDist
|
|
(
|
|
calcSnapDistance
|
|
(
|
|
mesh,
|
|
snapParams,
|
|
pp
|
|
)
|
|
);
|
|
|
|
const faceList& localFaces = pp.localFaces();
|
|
|
|
forAll(localFaces, faceI)
|
|
{
|
|
const face& f = localFaces[faceI];
|
|
|
|
forAll(f, fp)
|
|
{
|
|
faceSnapDist[faceI] = max
|
|
(
|
|
faceSnapDist[faceI],
|
|
snapDist[f[fp]]
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
pointField localFaceCentres(mesh.faceCentres(), pp.addressing());
|
|
|
|
// Get nearest surface and region
|
|
labelList hitSurface;
|
|
labelList hitRegion;
|
|
surfaces.findNearestRegion
|
|
(
|
|
unzonedSurfaces,
|
|
localFaceCentres,
|
|
sqr(faceSnapDist), // sqr of attract distance
|
|
hitSurface,
|
|
hitRegion
|
|
);
|
|
|
|
// Get patch
|
|
forAll(pp, i)
|
|
{
|
|
label faceI = pp.addressing()[i];
|
|
|
|
if (hitSurface[i] != -1 && !isZonedFace.get(faceI))
|
|
{
|
|
closestPatch[i] = globalToMasterPatch_
|
|
[
|
|
surfaces.globalRegion
|
|
(
|
|
hitSurface[i],
|
|
hitRegion[i]
|
|
)
|
|
];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Change those faces for which there is a different closest patch
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
labelList ownPatch(mesh.nFaces(), -1);
|
|
labelList neiPatch(mesh.nFaces(), -1);
|
|
|
|
const polyBoundaryMesh& patches = mesh.boundaryMesh();
|
|
|
|
forAll(patches, patchI)
|
|
{
|
|
const polyPatch& pp = patches[patchI];
|
|
|
|
forAll(pp, i)
|
|
{
|
|
ownPatch[pp.start()+i] = patchI;
|
|
neiPatch[pp.start()+i] = patchI;
|
|
}
|
|
}
|
|
|
|
label nChanged = 0;
|
|
forAll(closestPatch, i)
|
|
{
|
|
label faceI = pp.addressing()[i];
|
|
|
|
if (closestPatch[i] != -1 && closestPatch[i] != ownPatch[faceI])
|
|
{
|
|
ownPatch[faceI] = closestPatch[i];
|
|
neiPatch[faceI] = closestPatch[i];
|
|
nChanged++;
|
|
}
|
|
}
|
|
|
|
Info<< "Repatched " << returnReduce(nChanged, sumOp<label>())
|
|
<< " faces in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl
|
|
<< endl;
|
|
|
|
return meshRefiner_.createBaffles(ownPatch, neiPatch);
|
|
}
|
|
|
|
|
|
void Foam::snappySnapDriver::detectWarpedFaces
|
|
(
|
|
const scalar featureCos,
|
|
const indirectPrimitivePatch& pp,
|
|
|
|
DynamicList<label>& splitFaces,
|
|
DynamicList<labelPair>& splits
|
|
) const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
const faceList& localFaces = pp.localFaces();
|
|
const pointField& localPoints = pp.localPoints();
|
|
const labelList& bFaces = pp.addressing();
|
|
|
|
splitFaces.clear();
|
|
splitFaces.setCapacity(bFaces.size());
|
|
splits.clear();
|
|
splits.setCapacity(bFaces.size());
|
|
|
|
// Determine parallel consistent normals on points
|
|
const vectorField pointNormals(PatchTools::pointNormals(mesh, pp));
|
|
|
|
face f0(4);
|
|
face f1(4);
|
|
|
|
forAll(localFaces, faceI)
|
|
{
|
|
const face& f = localFaces[faceI];
|
|
|
|
if (f.size() >= 4)
|
|
{
|
|
// See if splitting face across diagonal would make two faces with
|
|
// biggish normal angle
|
|
|
|
labelPair minDiag(-1, -1);
|
|
scalar minCos(GREAT);
|
|
|
|
for (label startFp = 0; startFp < f.size()-2; startFp++)
|
|
{
|
|
label minFp = f.rcIndex(startFp);
|
|
|
|
for
|
|
(
|
|
label endFp = f.fcIndex(f.fcIndex(startFp));
|
|
endFp < f.size() && endFp != minFp;
|
|
endFp++
|
|
)
|
|
{
|
|
// Form two faces
|
|
f0.setSize(endFp-startFp+1);
|
|
label i0 = 0;
|
|
for (label fp = startFp; fp <= endFp; fp++)
|
|
{
|
|
f0[i0++] = f[fp];
|
|
}
|
|
f1.setSize(f.size()+2-f0.size());
|
|
label i1 = 0;
|
|
for (label fp = endFp; fp != startFp; fp = f.fcIndex(fp))
|
|
{
|
|
f1[i1++] = f[fp];
|
|
}
|
|
f1[i1++] = f[startFp];
|
|
|
|
//Info<< "Splitting face:" << f << " into f0:" << f0
|
|
// << " f1:" << f1 << endl;
|
|
|
|
vector n0 = f0.normal(localPoints);
|
|
scalar n0Mag = mag(n0);
|
|
vector n1 = f1.normal(localPoints);
|
|
scalar n1Mag = mag(n1);
|
|
|
|
if (n0Mag > ROOTVSMALL && n1Mag > ROOTVSMALL)
|
|
{
|
|
scalar cosAngle = (n0/n0Mag) & (n1/n1Mag);
|
|
if (cosAngle < minCos)
|
|
{
|
|
minCos = cosAngle;
|
|
minDiag = labelPair(startFp, endFp);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
if (minCos < featureCos)
|
|
{
|
|
splitFaces.append(bFaces[faceI]);
|
|
splits.append(minDiag);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::labelList Foam::snappySnapDriver::getInternalOrBaffleDuplicateFace() const
|
|
{
|
|
const fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
labelList internalOrBaffleFaceZones;
|
|
{
|
|
List<surfaceZonesInfo::faceZoneType> fzTypes(2);
|
|
fzTypes[0] = surfaceZonesInfo::INTERNAL;
|
|
fzTypes[1] = surfaceZonesInfo::BAFFLE;
|
|
internalOrBaffleFaceZones = meshRefiner_.getZones(fzTypes);
|
|
}
|
|
|
|
List<labelPair> baffles
|
|
(
|
|
meshRefiner_.subsetBaffles
|
|
(
|
|
mesh,
|
|
internalOrBaffleFaceZones,
|
|
localPointRegion::findDuplicateFacePairs(mesh)
|
|
)
|
|
);
|
|
|
|
labelList faceToDuplicate(mesh.nFaces(), -1);
|
|
forAll(baffles, i)
|
|
{
|
|
const labelPair& p = baffles[i];
|
|
faceToDuplicate[p[0]] = p[1];
|
|
faceToDuplicate[p[1]] = p[0];
|
|
}
|
|
|
|
return faceToDuplicate;
|
|
}
|
|
|
|
|
|
void Foam::snappySnapDriver::doSnap
|
|
(
|
|
const dictionary& snapDict,
|
|
const dictionary& motionDict,
|
|
const bool mergePatchFaces,
|
|
const scalar featureCos,
|
|
const scalar planarAngle,
|
|
const snapParameters& snapParams
|
|
)
|
|
{
|
|
fvMesh& mesh = meshRefiner_.mesh();
|
|
|
|
Info<< nl
|
|
<< "Morphing phase" << nl
|
|
<< "--------------" << nl
|
|
<< endl;
|
|
|
|
// faceZone handling
|
|
// ~~~~~~~~~~~~~~~~~
|
|
//
|
|
// We convert all faceZones into baffles during snapping so we can use
|
|
// a standard mesh motion (except for the mesh checking which for baffles
|
|
// created from internal faces should check across the baffles). The state
|
|
// is stored in two variables:
|
|
// baffles : pairs of boundary faces
|
|
// duplicateFace : from mesh face to its baffle colleague (or -1 for
|
|
// normal faces)
|
|
// There are three types of faceZones according to the faceType property:
|
|
//
|
|
// internal
|
|
// --------
|
|
// - baffles: need to be checked across
|
|
// - duplicateFace: from face to duplicate face. Contains
|
|
// all faces on faceZone to prevents merging patch faces.
|
|
//
|
|
// baffle
|
|
// ------
|
|
// - baffles: no need to be checked across
|
|
// - duplicateFace: contains all faces on faceZone to prevent
|
|
// merging patch faces.
|
|
//
|
|
// boundary
|
|
// --------
|
|
// - baffles: no need to be checked across. Also points get duplicated
|
|
// so will no longer be baffles
|
|
// - duplicateFace: contains no faces on faceZone since both sides can
|
|
// merge faces independently.
|
|
|
|
|
|
|
|
// faceZones of type internal
|
|
const labelList internalFaceZones
|
|
(
|
|
meshRefiner_.getZones
|
|
(
|
|
List<surfaceZonesInfo::faceZoneType>
|
|
(
|
|
1,
|
|
surfaceZonesInfo::INTERNAL
|
|
)
|
|
)
|
|
);
|
|
|
|
|
|
// Create baffles (pairs of faces that share the same points)
|
|
// Baffles stored as owner and neighbour face that have been created.
|
|
{
|
|
List<labelPair> baffles;
|
|
labelList originatingFaceZone;
|
|
meshRefiner_.createZoneBaffles
|
|
(
|
|
identity(mesh.faceZones().size()),
|
|
baffles,
|
|
originatingFaceZone
|
|
);
|
|
}
|
|
|
|
// Duplicate points on faceZones of type boundary
|
|
meshRefiner_.dupNonManifoldBoundaryPoints();
|
|
|
|
|
|
bool doFeatures = false;
|
|
label nFeatIter = 1;
|
|
if (snapParams.nFeatureSnap() > 0)
|
|
{
|
|
doFeatures = true;
|
|
nFeatIter = snapParams.nFeatureSnap();
|
|
|
|
Info<< "Snapping to features in " << nFeatIter
|
|
<< " iterations ..." << endl;
|
|
}
|
|
|
|
|
|
bool meshOk = false;
|
|
|
|
|
|
// Get the labels of added patches.
|
|
labelList adaptPatchIDs(meshRefiner_.meshedPatches());
|
|
|
|
|
|
|
|
{
|
|
autoPtr<indirectPrimitivePatch> ppPtr
|
|
(
|
|
meshRefinement::makePatch
|
|
(
|
|
mesh,
|
|
adaptPatchIDs
|
|
)
|
|
);
|
|
|
|
|
|
// Distance to attract to nearest feature on surface
|
|
scalarField snapDist(calcSnapDistance(mesh, snapParams, ppPtr()));
|
|
|
|
|
|
// Construct iterative mesh mover.
|
|
Info<< "Constructing mesh displacer ..." << endl;
|
|
Info<< "Using mesh parameters " << motionDict << nl << endl;
|
|
|
|
autoPtr<motionSmoother> meshMoverPtr
|
|
(
|
|
new motionSmoother
|
|
(
|
|
mesh,
|
|
ppPtr(),
|
|
adaptPatchIDs,
|
|
meshRefinement::makeDisplacementField
|
|
(
|
|
pointMesh::New(mesh),
|
|
adaptPatchIDs
|
|
),
|
|
motionDict
|
|
)
|
|
);
|
|
|
|
|
|
// Check initial mesh
|
|
Info<< "Checking initial mesh ..." << endl;
|
|
labelHashSet wrongFaces(mesh.nFaces()/100);
|
|
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
|
|
const label nInitErrors = returnReduce
|
|
(
|
|
wrongFaces.size(),
|
|
sumOp<label>()
|
|
);
|
|
|
|
Info<< "Detected " << nInitErrors << " illegal faces"
|
|
<< " (concave, zero area or negative cell pyramid volume)"
|
|
<< endl;
|
|
|
|
|
|
Info<< "Checked initial mesh in = "
|
|
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
|
|
|
|
// Extract baffles across internal faceZones (for checking mesh quality
|
|
// across
|
|
labelPairList internalBaffles
|
|
(
|
|
meshRefiner_.subsetBaffles
|
|
(
|
|
mesh,
|
|
internalFaceZones,
|
|
localPointRegion::findDuplicateFacePairs(mesh)
|
|
)
|
|
);
|
|
|
|
|
|
|
|
// Pre-smooth patch vertices (so before determining nearest)
|
|
preSmoothPatch
|
|
(
|
|
meshRefiner_,
|
|
snapParams,
|
|
nInitErrors,
|
|
internalBaffles,
|
|
meshMoverPtr()
|
|
);
|
|
|
|
|
|
|
|
//- Only if in feature attraction mode:
|
|
// Nearest feature
|
|
vectorField patchAttraction;
|
|
// Constraints at feature
|
|
List<pointConstraint> patchConstraints;
|
|
|
|
|
|
//- Any faces to split
|
|
DynamicList<label> splitFaces;
|
|
//- Indices in face to split across
|
|
DynamicList<labelPair> splits;
|
|
|
|
|
|
for (label iter = 0; iter < nFeatIter; iter++)
|
|
{
|
|
Info<< nl
|
|
<< "Morph iteration " << iter << nl
|
|
<< "-----------------" << endl;
|
|
|
|
// Splitting iteration?
|
|
bool doSplit = false;
|
|
if
|
|
(
|
|
doFeatures
|
|
&& snapParams.nFaceSplitInterval() > 0
|
|
&& (
|
|
(iter == nFeatIter-1)
|
|
|| (iter > 0 && (iter%snapParams.nFaceSplitInterval()) == 0)
|
|
)
|
|
)
|
|
{
|
|
doSplit = true;
|
|
}
|
|
|
|
|
|
|
|
indirectPrimitivePatch& pp = ppPtr();
|
|
motionSmoother& meshMover = meshMoverPtr();
|
|
|
|
|
|
// Calculate displacement at every patch point. Insert into
|
|
// meshMover.
|
|
// Calculate displacement at every patch point
|
|
pointField nearestPoint;
|
|
vectorField nearestNormal;
|
|
|
|
if (snapParams.detectNearSurfacesSnap())
|
|
{
|
|
nearestPoint.setSize(pp.nPoints(), vector::max);
|
|
nearestNormal.setSize(pp.nPoints(), Zero);
|
|
}
|
|
|
|
vectorField disp = calcNearestSurface
|
|
(
|
|
snapParams.strictRegionSnap(), // attract points to region only
|
|
meshRefiner_,
|
|
globalToMasterPatch_, // for if strictRegionSnap
|
|
globalToSlavePatch_, // for if strictRegionSnap
|
|
snapDist,
|
|
pp,
|
|
|
|
nearestPoint,
|
|
nearestNormal
|
|
);
|
|
|
|
|
|
// Override displacement at thin gaps
|
|
if (snapParams.detectNearSurfacesSnap())
|
|
{
|
|
detectNearSurfaces
|
|
(
|
|
Foam::cos(degToRad(planarAngle)),// planar cos for gaps
|
|
pp,
|
|
nearestPoint, // surfacepoint from nearest test
|
|
nearestNormal, // surfacenormal from nearest test
|
|
|
|
disp
|
|
);
|
|
}
|
|
|
|
// Override displacement with feature edge attempt
|
|
if (doFeatures)
|
|
{
|
|
splitFaces.clear();
|
|
splits.clear();
|
|
disp = calcNearestSurfaceFeature
|
|
(
|
|
snapParams,
|
|
!doSplit, // alignMeshEdges
|
|
iter,
|
|
featureCos,
|
|
scalar(iter+1)/nFeatIter,
|
|
|
|
snapDist,
|
|
disp,
|
|
nearestNormal,
|
|
meshMover,
|
|
|
|
patchAttraction,
|
|
patchConstraints,
|
|
|
|
splitFaces,
|
|
splits
|
|
);
|
|
}
|
|
|
|
// Check for displacement being outwards.
|
|
outwardsDisplacement(pp, disp);
|
|
|
|
// Set initial distribution of displacement field (on patches)
|
|
// from patchDisp and make displacement consistent with b.c.
|
|
// on displacement pointVectorField.
|
|
meshMover.setDisplacement(disp);
|
|
|
|
|
|
if (debug&meshRefinement::ATTRACTION)
|
|
{
|
|
dumpMove
|
|
(
|
|
mesh.time().path()
|
|
/ "patchDisplacement_" + name(iter) + ".obj",
|
|
pp.localPoints(),
|
|
pp.localPoints() + disp
|
|
);
|
|
}
|
|
|
|
// Get smoothly varying internal displacement field.
|
|
smoothDisplacement(snapParams, meshMover);
|
|
|
|
// Apply internal displacement to mesh.
|
|
meshOk = scaleMesh
|
|
(
|
|
snapParams,
|
|
nInitErrors,
|
|
internalBaffles,
|
|
meshMover
|
|
);
|
|
|
|
if (!meshOk)
|
|
{
|
|
WarningInFunction
|
|
<< "Did not succesfully snap mesh."
|
|
<< " Continuing to snap to resolve easy" << nl
|
|
<< " surfaces but the"
|
|
<< " resulting mesh will not satisfy your quality"
|
|
<< " constraints" << nl << endl;
|
|
}
|
|
|
|
if (debug&meshRefinement::MESH)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing scaled mesh to time "
|
|
<< meshRefiner_.timeName() << endl;
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
mesh.time().path()/meshRefiner_.timeName()
|
|
);
|
|
Info<< "Writing displacement field ..." << endl;
|
|
meshMover.displacement().write();
|
|
tmp<pointScalarField> magDisp
|
|
(
|
|
mag(meshMover.displacement())
|
|
);
|
|
magDisp().write();
|
|
}
|
|
|
|
// Use current mesh as base mesh
|
|
meshMover.correct();
|
|
|
|
|
|
|
|
// See if any faces need splitting
|
|
label nTotalSplit = returnReduce(splitFaces.size(), sumOp<label>());
|
|
if (nTotalSplit && doSplit)
|
|
{
|
|
// Filter out baffle faces from faceZones of type
|
|
// internal/baffle
|
|
|
|
labelList duplicateFace(getInternalOrBaffleDuplicateFace());
|
|
|
|
{
|
|
labelList oldSplitFaces(splitFaces.xfer());
|
|
List<labelPair> oldSplits(splits.xfer());
|
|
forAll(oldSplitFaces, i)
|
|
{
|
|
if (duplicateFace[oldSplitFaces[i]] == -1)
|
|
{
|
|
splitFaces.append(oldSplitFaces[i]);
|
|
splits.append(oldSplits[i]);
|
|
}
|
|
}
|
|
nTotalSplit = returnReduce
|
|
(
|
|
splitFaces.size(),
|
|
sumOp<label>()
|
|
);
|
|
}
|
|
|
|
// Update mesh
|
|
meshRefiner_.splitFacesUndo
|
|
(
|
|
splitFaces,
|
|
splits,
|
|
motionDict,
|
|
|
|
duplicateFace,
|
|
internalBaffles
|
|
);
|
|
|
|
// Redo meshMover
|
|
meshMoverPtr.clear();
|
|
ppPtr.clear();
|
|
|
|
// Update mesh mover
|
|
ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
|
|
meshMoverPtr.reset
|
|
(
|
|
new motionSmoother
|
|
(
|
|
mesh,
|
|
ppPtr(),
|
|
adaptPatchIDs,
|
|
meshRefinement::makeDisplacementField
|
|
(
|
|
pointMesh::New(mesh),
|
|
adaptPatchIDs
|
|
),
|
|
motionDict
|
|
)
|
|
);
|
|
|
|
// Update snapping distance
|
|
snapDist = calcSnapDistance(mesh, snapParams, ppPtr());
|
|
|
|
|
|
if (debug&meshRefinement::MESH)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing split-faces mesh to time "
|
|
<< meshRefiner_.timeName() << endl;
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
mesh.time().path()/meshRefiner_.timeName()
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
if (debug&meshRefinement::MESH)
|
|
{
|
|
forAll(internalBaffles, i)
|
|
{
|
|
const labelPair& p = internalBaffles[i];
|
|
const point& fc0 = mesh.faceCentres()[p[0]];
|
|
const point& fc1 = mesh.faceCentres()[p[1]];
|
|
|
|
if (mag(fc0-fc1) > meshRefiner_.mergeDistance())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Separated baffles : f0:" << p[0]
|
|
<< " centre:" << fc0
|
|
<< " f1:" << p[1] << " centre:" << fc1
|
|
<< " distance:" << mag(fc0-fc1)
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Merge any introduced baffles (from faceZones of faceType 'internal')
|
|
{
|
|
autoPtr<mapPolyMesh> mapPtr = meshRefiner_.mergeZoneBaffles
|
|
(
|
|
true, // internal zones
|
|
false // baffle zones
|
|
);
|
|
|
|
if (mapPtr.valid())
|
|
{
|
|
if (debug & meshRefinement::MESH)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing baffle-merged mesh to time "
|
|
<< meshRefiner_.timeName() << endl;
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
meshRefiner_.timeName()
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Repatch faces according to nearest. Do not repatch baffle faces.
|
|
{
|
|
labelList duplicateFace(getInternalOrBaffleDuplicateFace());
|
|
|
|
repatchToSurface(snapParams, adaptPatchIDs, duplicateFace);
|
|
}
|
|
|
|
if (mergePatchFaces)
|
|
{
|
|
labelList duplicateFace(getInternalOrBaffleDuplicateFace());
|
|
|
|
// Repatching might have caused faces to be on same patch and hence
|
|
// mergeable so try again to merge coplanar faces. Do not merge baffle
|
|
// faces to ensure they both stay the same.
|
|
label nChanged = meshRefiner_.mergePatchFacesUndo
|
|
(
|
|
featureCos, // minCos
|
|
featureCos, // concaveCos
|
|
meshRefiner_.meshedPatches(),
|
|
motionDict,
|
|
duplicateFace // faces not to merge
|
|
);
|
|
|
|
nChanged += meshRefiner_.mergeEdgesUndo(featureCos, motionDict);
|
|
|
|
if (nChanged > 0 && debug & meshRefinement::MESH)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
Info<< "Writing patchFace merged mesh to time "
|
|
<< meshRefiner_.timeName() << endl;
|
|
meshRefiner_.write
|
|
(
|
|
meshRefinement::debugType(debug),
|
|
meshRefinement::writeType
|
|
(
|
|
meshRefinement::writeLevel()
|
|
| meshRefinement::WRITEMESH
|
|
),
|
|
meshRefiner_.timeName()
|
|
);
|
|
}
|
|
}
|
|
|
|
if (debug & meshRefinement::MESH)
|
|
{
|
|
const_cast<Time&>(mesh.time())++;
|
|
}
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|