zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: snappyHexMesh: detect small distance snapping

Browse Source
This commit is contained in:
mattijs
2013-09-13 16:18:04 +01:00
parent 1348eabe3b
commit 1a0a6b51dd
8 changed files with 987 additions and 64 deletions
Download Patch File Download Diff File Expand all files Collapse all files

750
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
View File

@ -792,11 +792,636 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
} }
Foam::tmp<Foam::pointField> Foam::autoSnapDriver::avgCellCentres
(
const fvMesh& mesh,
const indirectPrimitivePatch& pp
)
{
const labelListList& pointFaces = pp.pointFaces();
tmp<pointField> tavgBoundary
(
new pointField(pointFaces.size(), vector::zero)
);
pointField& avgBoundary = tavgBoundary();
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::autoSnapDriver::calcEdgeLen
//(
// const indirectPrimitivePatch& pp
//) const
//{
// // Get local edge length based on refinement level
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// // (Ripped from autoLayerDriver)
//
// 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::autoSnapDriver::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 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 from localPoints to beyond cell centre
pointField end(pp.nPoints());
forAll(localPoints, pointI)
{
const point& pt = localPoints[pointI];
end[pointI] = pt + 2*(avgCc[pointI]-pt);
}
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,
localPoints,
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 (override)
{
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(localPoints, zonePointIndices),
pointField(end, zonePointIndices),
surface1,
hit1,
region1,
normal1,
surface2,
hit2,
region2,
normal2
);
label nOverride = 0;
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 (override)
{
nOverride++;
}
}
}
}
Info<< "Overriding nearest with intersection of close gaps at "
<< returnReduce(nOverride, sumOp<label>())
<< " out of " << returnReduce(pp.nPoints(), sumOp<label>())
<< " points." << endl;
}
Foam::vectorField Foam::autoSnapDriver::calcNearestSurface Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
( (
const meshRefinement& meshRefiner, const meshRefinement& meshRefiner,
const scalarField& snapDist, const scalarField& snapDist,
const indirectPrimitivePatch& pp const indirectPrimitivePatch& pp,
pointField& nearestPoint,
vectorField& nearestNormal
) )
{ {
Info<< "Calculating patchDisplacement as distance to nearest surface" Info<< "Calculating patchDisplacement as distance to nearest surface"
@ -815,12 +1440,12 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
labelList snapSurf(localPoints.size(), -1); labelList snapSurf(localPoints.size(), -1);
// Divide surfaces into zoned and unzoned // Divide surfaces into zoned and unzoned
labelList zonedSurfaces = const labelList zonedSurfaces =
surfaceZonesInfo::getNamedSurfaces surfaceZonesInfo::getNamedSurfaces
( (
meshRefiner.surfaces().surfZones() meshRefiner.surfaces().surfZones()
); );
labelList unzonedSurfaces = const labelList unzonedSurfaces =
surfaceZonesInfo::getUnnamedSurfaces surfaceZonesInfo::getUnnamedSurfaces
( (
meshRefiner.surfaces().surfZones() meshRefiner.surfaces().surfZones()
@ -833,14 +1458,42 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
{ {
List<pointIndexHit> hitInfo; List<pointIndexHit> hitInfo;
labelList hitSurface; labelList hitSurface;
surfaces.findNearest
( if (nearestNormal.size() == localPoints.size())
unzonedSurfaces, {
localPoints, labelList hitRegion;
sqr(snapDist), // sqr of attract distance vectorField hitNormal;
hitSurface, surfaces.findNearestRegion
hitInfo (
); 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) forAll(hitInfo, pointI)
{ {
@ -888,14 +1541,43 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
// Find nearest for points both on faceZone and pp. // Find nearest for points both on faceZone and pp.
List<pointIndexHit> hitInfo; List<pointIndexHit> hitInfo;
labelList hitSurface; labelList hitSurface;
surfaces.findNearest
( if (nearestNormal.size() == localPoints.size())
labelList(1, zoneSurfI), {
pointField(localPoints, zonePointIndices), labelList hitRegion;
sqr(scalarField(minSnapDist, zonePointIndices)), vectorField hitNormal;
hitSurface, surfaces.findNearestRegion
hitInfo (
); surfacesToTest,
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,
pointField(localPoints, zonePointIndices),
sqr(scalarField(minSnapDist, zonePointIndices)),
hitSurface,
hitInfo
);
}
forAll(hitInfo, i) forAll(hitInfo, i)
{ {
@ -1529,6 +2211,7 @@ void Foam::autoSnapDriver::doSnap
const dictionary& snapDict, const dictionary& snapDict,
const dictionary& motionDict, const dictionary& motionDict,
const scalar featureCos, const scalar featureCos,
const scalar planarAngle,
const snapParameters& snapParams const snapParameters& snapParams
) )
{ {
@ -1791,13 +2474,40 @@ void Foam::autoSnapDriver::doSnap
// Calculate displacement at every patch point. Insert into // Calculate displacement at every patch point. Insert into
// meshMover. // meshMover.
// Calculate displacement at every patch point
pointField nearestPoint;
vectorField nearestNormal;
if (snapParams.detectNearSurfacesSnap())
{
nearestPoint.setSize(pp.nPoints(), vector::max);
nearestNormal.setSize(pp.nPoints(), vector::zero);
}
vectorField disp = calcNearestSurface vectorField disp = calcNearestSurface
( (
meshRefiner_, meshRefiner_,
snapDist, snapDist,
meshMover.patch() 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 // Override displacement with feature edge attempt
if (doFeatures) if (doFeatures)
{ {

23
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
View File

@ -148,6 +148,7 @@ class autoSnapDriver
( (
const label iter, const label iter,
const indirectPrimitivePatch& pp, const indirectPrimitivePatch& pp,
const scalarField& faceSnapDist,
vectorField& faceDisp, vectorField& faceDisp,
vectorField& faceSurfaceNormal, vectorField& faceSurfaceNormal,
labelList& faceSurfaceRegion, labelList& faceSurfaceRegion,
@ -432,6 +433,23 @@ public:
const word& zoneName const word& zoneName
); );
//- Helper: calculate average cell centre per point
static tmp<pointField> avgCellCentres
(
const fvMesh& mesh,
const indirectPrimitivePatch&
);
//- Per patch point override displacement if in gap situation
void detectNearSurfaces
(
const scalar planarCos,
const indirectPrimitivePatch&,
const pointField& nearestPoint,
const vectorField& nearestNormal,
vectorField& disp
) const;
//- Per patch point calculate point on nearest surface. Set as //- Per patch point calculate point on nearest surface. Set as
// boundary conditions of motionSmoother displacement field. Return // boundary conditions of motionSmoother displacement field. Return
// displacement of patch points. // displacement of patch points.
@ -439,7 +457,9 @@ public:
( (
const meshRefinement& meshRefiner, const meshRefinement& meshRefiner,
const scalarField& snapDist, const scalarField& snapDist,
const indirectPrimitivePatch& const indirectPrimitivePatch&,
pointField& nearestPoint,
vectorField& nearestNormal
); );
////- Per patch point calculate point on nearest surface. Set as ////- Per patch point calculate point on nearest surface. Set as
@ -483,6 +503,7 @@ public:
const dictionary& snapDict, const dictionary& snapDict,
const dictionary& motionDict, const dictionary& motionDict,
const scalar featureCos, const scalar featureCos,
const scalar planarAngle,
const snapParameters& snapParams const snapParameters& snapParams
); );

168
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
View File

@ -319,6 +319,7 @@ void Foam::autoSnapDriver::calcNearestFace
( (
const label iter, const label iter,
const indirectPrimitivePatch& pp, const indirectPrimitivePatch& pp,
const scalarField& faceSnapDist,
vectorField& faceDisp, vectorField& faceDisp,
vectorField& faceSurfaceNormal, vectorField& faceSurfaceNormal,
labelList& faceSurfaceGlobalRegion, labelList& faceSurfaceGlobalRegion,
@ -409,7 +410,7 @@ void Foam::autoSnapDriver::calcNearestFace
( (
labelList(1, zoneSurfI), labelList(1, zoneSurfI),
fc, fc,
sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist sqr(faceSnapDist),// sqr of attract dist
hitSurface, hitSurface,
hitInfo, hitInfo,
hitRegion, hitRegion,
@ -429,6 +430,14 @@ void Foam::autoSnapDriver::calcNearestFace
hitRegion[hitI] hitRegion[hitI]
); );
} }
else
{
WarningIn
(
"autoSnapDriver::calcNearestFace(..)"
) << "Did not find surface near face centre " << fc[hitI]
<< endl;
}
} }
} }
@ -467,7 +476,7 @@ void Foam::autoSnapDriver::calcNearestFace
( (
unzonedSurfaces, unzonedSurfaces,
fc, fc,
sqr(scalarField(fc.size(), GREAT)),// sqr of attract dist sqr(faceSnapDist),// sqr of attract dist
hitSurface, hitSurface,
hitInfo, hitInfo,
hitRegion, hitRegion,
@ -487,6 +496,14 @@ void Foam::autoSnapDriver::calcNearestFace
hitRegion[hitI] hitRegion[hitI]
); );
} }
else
{
WarningIn
(
"autoSnapDriver::calcNearestFace(..)"
) << "Did not find surface near face centre " << fc[hitI]
<< endl;
}
} }
@ -560,24 +577,44 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
forAll(pp.pointFaces(), pointI) forAll(pp.pointFaces(), pointI)
{ {
const labelList& pFaces = pp.pointFaces()[pointI]; const labelList& pFaces = pp.pointFaces()[pointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
pNormals.setSize(pFaces.size());
List<point>& pDisp = pointFaceDisp[pointI];
pDisp.setSize(pFaces.size());
List<point>& pFc = pointFaceCentres[pointI];
pFc.setSize(pFaces.size());
labelList& pFid = pointFacePatchID[pointI];
pFid.setSize(pFaces.size());
// Count valid face normals
label nFaces = 0;
forAll(pFaces, i) forAll(pFaces, i)
{ {
label faceI = pFaces[i]; label faceI = pFaces[i];
pNormals[i] = faceSurfaceNormal[faceI]; if (faceSurfaceGlobalRegion[faceI] != -1)
pDisp[i] = faceDisp[faceI]; {
pFc[i] = pp.faceCentres()[faceI]; nFaces++;
//label meshFaceI = pp.addressing()[faceI]; }
//pFid[i] = mesh.boundaryMesh().whichPatch(meshFaceI); }
pFid[i] = globalToMasterPatch_[faceSurfaceGlobalRegion[faceI]];
List<point>& pNormals = pointFaceSurfNormals[pointI];
pNormals.setSize(nFaces);
List<point>& pDisp = pointFaceDisp[pointI];
pDisp.setSize(nFaces);
List<point>& pFc = pointFaceCentres[pointI];
pFc.setSize(nFaces);
labelList& pFid = pointFacePatchID[pointI];
pFid.setSize(nFaces);
nFaces = 0;
forAll(pFaces, i)
{
label faceI = pFaces[i];
label globalRegionI = faceSurfaceGlobalRegion[faceI];
if (globalRegionI != -1)
{
pNormals[nFaces] = faceSurfaceNormal[faceI];
pDisp[nFaces] = faceDisp[faceI];
pFc[nFaces] = pp.faceCentres()[faceI];
//label meshFaceI = pp.addressing()[faceI];
//pFid[nFaces] = mesh.boundaryMesh().whichPatch(meshFaceI);
pFid[nFaces] = globalToMasterPatch_[globalRegionI];
nFaces++;
}
} }
} }
@ -2123,9 +2160,85 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
} }
// Collect candidate points for attraction
DynamicList<label> attractPoints(pp.nPoints());
{
const fvMesh& mesh = meshRefiner_.mesh();
boolList isFeatureEdgeOrPoint(pp.nPoints(), false);
label nFeats = 0;
forAll(rawPatchConstraints, pointI)
{
if (rawPatchConstraints[pointI].first() >= 2)
{
isFeatureEdgeOrPoint[pointI] = true;
nFeats++;
}
}
Info<< "Intially selected " << returnReduce(nFeats, sumOp<label>())
<< " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
<< " for reverse attraction." << endl;
// Make sure is synchronised (note: check if constraint is already
// synced in which case this is not needed here)
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
isFeatureEdgeOrPoint,
orEqOp<bool>(), // combine op
false
);
for (label nGrow = 0; nGrow < 1; nGrow++)
{
forAll(pp.localFaces(), faceI)
{
const face& f = pp.localFaces()[faceI];
forAll(f, fp)
{
if (isFeatureEdgeOrPoint[f[fp]])
{
// Mark all points on face
forAll(f, fp)
{
isFeatureEdgeOrPoint[f[fp]] = true;
}
break;
}
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
isFeatureEdgeOrPoint,
orEqOp<bool>(), // combine op
false
);
}
// Collect attractPoints
forAll(isFeatureEdgeOrPoint, pointI)
{
if (isFeatureEdgeOrPoint[pointI])
{
attractPoints.append(pointI);
}
}
Info<< "Selected " << returnReduce(attractPoints.size(), sumOp<label>())
<< " points out of " << returnReduce(pp.nPoints(), sumOp<label>())
<< " for reverse attraction." << endl;
}
indexedOctree<treeDataPoint> ppTree indexedOctree<treeDataPoint> ppTree
( (
treeDataPoint(pp.localPoints()), treeDataPoint(pp.localPoints(), attractPoints),
bb, // overall search domain bb, // overall search domain
8, // maxLevel 8, // maxLevel
10, // leafsize 10, // leafsize
@ -2159,7 +2272,8 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
if (nearInfo.hit()) if (nearInfo.hit())
{ {
label pointI = nearInfo.index(); label pointI =
ppTree.shapes().pointLabels()[nearInfo.index()];
const point attraction = featPt-pp.localPoints()[pointI]; const point attraction = featPt-pp.localPoints()[pointI];
// Check if this point is already being attracted. If so // Check if this point is already being attracted. If so
@ -2220,7 +2334,9 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
if (nearInfo.hit()) if (nearInfo.hit())
{ {
label pointI = nearInfo.index(); label pointI =
ppTree.shapes().pointLabels()[nearInfo.index()];
const point& pt = pp.localPoints()[pointI]; const point& pt = pp.localPoints()[pointI];
const point attraction = featPt-pt; const point attraction = featPt-pt;
@ -2681,6 +2797,19 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
const fvMesh& mesh = meshRefiner_.mesh(); const fvMesh& mesh = meshRefiner_.mesh();
// Calculate attraction distance per face (from the attraction distance
// per point)
scalarField faceSnapDist(pp.size(), -GREAT);
forAll(pp.localFaces(), faceI)
{
const face& f = pp.localFaces()[faceI];
forAll(f, fp)
{
faceSnapDist[faceI] = max(faceSnapDist[faceI], snapDist[f[fp]]);
}
}
// Displacement and orientation per pp face // Displacement and orientation per pp face
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -2695,6 +2824,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
( (
iter, iter,
pp, pp,
faceSnapDist,
faceDisp, faceDisp,
faceSurfaceNormal, faceSurfaceNormal,
faceSurfaceGlobalRegion, faceSurfaceGlobalRegion,

9
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -40,11 +40,12 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
multiRegionFeatureSnap_ multiRegionFeatureSnap_
( (
dict.lookupOrDefault("multiRegionFeatureSnap", false) dict.lookupOrDefault("multiRegionFeatureSnap", false)
),
detectNearSurfacesSnap_
(
dict.lookupOrDefault("detectNearSurfacesSnap", true)
) )
{} {}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* // // ************************************************************************* //

9
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -70,6 +70,8 @@ class snapParameters
const Switch multiRegionFeatureSnap_; const Switch multiRegionFeatureSnap_;
const Switch detectNearSurfacesSnap_;
// Private Member Functions // Private Member Functions
@ -141,6 +143,11 @@ public:
return multiRegionFeatureSnap_; return multiRegionFeatureSnap_;
} }
Switch detectNearSurfacesSnap() const
{
return detectNearSurfacesSnap_;
}
}; };

30
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
View File

@ -289,16 +289,6 @@ private:
label& nRefine label& nRefine
) const; ) const;
//- Is local topology a small gap?
bool isGap
(
const scalar,
const vector&,
const vector&,
const vector&,
const vector&
) const;
//- Mark cell if local a gap topology or //- Mark cell if local a gap topology or
bool checkProximity bool checkProximity
( (
@ -664,6 +654,26 @@ public:
// Refinement // Refinement
//- Is local topology a small gap?
bool isGap
(
const scalar,
const vector&,
const vector&,
const vector&,
const vector&
) const;
//- Is local topology a small gap normal to the test vector
bool isNormalGap
(
const scalar,
const vector&,
const vector&,
const vector&,
const vector&
) const;
//- Calculate list of cells to refine. //- Calculate list of cells to refine.
labelList refineCandidates labelList refineCandidates
( (

6
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
View File

@ -1061,13 +1061,17 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
meshMover meshMover
); );
pointField nearestPoint;
vectorField nearestNormal;
const vectorField disp const vectorField disp
( (
autoSnapDriver::calcNearestSurface autoSnapDriver::calcNearestSurface
( (
*this, *this,
snapDist, // attraction snapDist, // attraction
pp pp,
nearestPoint,
nearestNormal
) )
); );

56
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -1214,7 +1214,6 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
} }
// Mark small gaps
bool Foam::meshRefinement::isGap bool Foam::meshRefinement::isGap
( (
const scalar planarCos, const scalar planarCos,
@ -1225,12 +1224,53 @@ bool Foam::meshRefinement::isGap
const vector& normal1 const vector& normal1
) const ) const
{ {
////- hits differ and angles are oppositeish //- hits differ and angles are oppositeish and
//return // hits have a normal distance
// (mag(point0-point1) > mergeDistance()) vector d = point1-point0;
// && ((normal0 & normal1) < (-1+planarCos)); scalar magD = mag(d);
if (magD > mergeDistance())
{
scalar cosAngle = (normal0 & normal1);
vector avg = vector::zero;
if (cosAngle < (-1+planarCos))
{
// Opposite normals
avg = 0.5*(normal0-normal1);
}
else if (cosAngle > (1-planarCos))
{
avg = 0.5*(normal0+normal1);
}
if (avg != vector::zero)
{
return true;
}
else
{
return false;
}
}
else
{
return false;
}
}
// Mark small gaps
bool Foam::meshRefinement::isNormalGap
(
const scalar planarCos,
const vector& point0,
const vector& normal0,
const vector& point1,
const vector& normal1
) const
{
//- hits differ and angles are oppositeish and //- hits differ and angles are oppositeish and
// hits have a normal distance // hits have a normal distance
vector d = point1-point0; vector d = point1-point0;
@ -1316,7 +1356,7 @@ bool Foam::meshRefinement::checkProximity
// - different location // - different location
// - opposite surface // - opposite surface
bool closeSurfaces = isGap bool closeSurfaces = isNormalGap
( (
planarCos, planarCos,
cellMaxLocation, cellMaxLocation,
@ -1565,7 +1605,7 @@ Foam::label Foam::meshRefinement::markProximityRefinement
// Have valid data on both sides. Check planarCos. // Have valid data on both sides. Check planarCos.
if if
( (
isGap isNormalGap
( (
planarCos, planarCos,
cellMaxLocation[own], cellMaxLocation[own],
@ -1635,7 +1675,7 @@ Foam::label Foam::meshRefinement::markProximityRefinement
// Have valid data on both sides. Check planarCos. // Have valid data on both sides. Check planarCos.
if if
( (
isGap isNormalGap
( (
planarCos, planarCos,
cellMaxLocation[own], cellMaxLocation[own],