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: autoHexMesh: parallel consistency, fix curvature refinement

Browse Source
This commit is contained in:
mattijs
2013-10-18 11:52:38 +01:00
parent a76199f37c
commit 71c0a5d1d7
5 changed files with 652 additions and 272 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -910,6 +910,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
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();
@ -1220,6 +1221,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
}
const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
label nOverride = 0;
// 1. All points to non-interface surfaces
@ -1332,7 +1334,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
}
}
if (override)
if (override && isMasterPoint[meshPoints[pointI]])
{
nOverride++;
}
@ -1399,8 +1401,6 @@ void Foam::autoSnapDriver::detectNearSurfaces
);
label nOverride = 0;
forAll(hit1, i)
{
label pointI = zonePointIndices[i];
@ -1472,7 +1472,7 @@ void Foam::autoSnapDriver::detectNearSurfaces
}
}
if (override)
if (override && isMasterPoint[meshPoints[pointI]])
{
nOverride++;
}
@ -1690,7 +1690,13 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
scalarField magDisp(mag(patchDisp));
Info<< "Wanted displacement : average:"
<< gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
<< meshRefinement::gAverage
(
mesh,
syncTools::getMasterPoints(mesh),
pp.meshPoints(),
magDisp
)
<< " min:" << gMin(magDisp)
<< " max:" << gMax(magDisp) << endl;
}
@ -2548,25 +2554,30 @@ void Foam::autoSnapDriver::doSnap
adaptPatchIDs
)
);
indirectPrimitivePatch& pp = ppPtr();
// Distance to attract to nearest feature on surface
const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp));
const scalarField snapDist(calcSnapDistance(mesh, snapParams, ppPtr()));
// Construct iterative mesh mover.
Info<< "Constructing mesh displacer ..." << endl;
Info<< "Using mesh parameters " << motionDict << nl << endl;
const pointMesh& pMesh = pointMesh::New(mesh);
motionSmoother meshMover
autoPtr<motionSmoother> meshMoverPtr
(
mesh,
pp,
adaptPatchIDs,
meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs),
motionDict
new motionSmoother
(
mesh,
ppPtr(),
adaptPatchIDs,
meshRefinement::makeDisplacementField
(
pointMesh::New(mesh),
adaptPatchIDs
),
motionDict
)
);
@ -2595,16 +2606,95 @@ void Foam::autoSnapDriver::doSnap
snapParams,
nInitErrors,
baffles,
meshMover
meshMoverPtr()
);
//- Only if in feature attraction mode:
// Nearest feature
vectorField patchAttraction;
// Constraints at feature
List<pointConstraint> patchConstraints;
for (label iter = 0; iter < nFeatIter; iter++)
{
Info<< nl
<< "Morph iteration " << iter << nl
<< "-----------------" << endl;
//if (iter > 0 && iter == nFeatIter/2)
//{
// Info<< "Splitting diagonal attractions" << endl;
// const labelList& bFaces = ppPtr().addressing();
//
// DynamicList<label> splitFaces(bFaces.size());
// DynamicList<labelPair> splits(bFaces.size());
//
// forAll(bFaces, faceI)
// {
// const labelPair split
// (
// findDiagonalAttraction
// (
// ppPtr(),
// patchAttraction,
// patchConstraints,
// faceI
// )
// );
//
// if (split != labelPair(-1, -1))
// {
// splitFaces.append(bFaces[faceI]);
// splits.append(split);
// }
// }
//
// autoPtr<mapPolyMesh> mapPtr = meshRefiner_.splitFaces
// (
// splitFaces,
// splits
// );
//
// const labelList& faceMap = mapPtr().faceMap();
// meshRefinement::updateList(faceMap, -1, duplicateFace);
// const labelList& reverseFaceMap = mapPtr().reverseFaceMap();
// forAll(baffles, i)
// {
// labelPair& baffle = baffles[i];
// baffle.first() = reverseFaceMap[baffle.first()];
// baffle.second() = reverseFaceMap[baffle.second()];
// }
//
// meshMoverPtr.clear();
// ppPtr.clear();
//
// ppPtr = meshRefinement::makePatch(mesh, adaptPatchIDs);
// meshMoverPtr.reset
// (
// new motionSmoother
// (
// mesh,
// ppPtr(),
// adaptPatchIDs,
// meshRefinement::makeDisplacementField
// (
// pointMesh::New(mesh),
// adaptPatchIDs
// ),
// motionDict
// )
// );
//}
indirectPrimitivePatch& pp = ppPtr();
motionSmoother& meshMover = meshMoverPtr();
// Calculate displacement at every patch point. Insert into
// meshMover.
// Calculate displacement at every patch point
@ -2652,7 +2742,9 @@ void Foam::autoSnapDriver::doSnap
scalar(iter+1)/nFeatIter,
snapDist,
disp,
meshMover
meshMover,
patchAttraction,
patchConstraints
);
}

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

@ -125,7 +125,9 @@ class autoSnapDriver
void smoothAndConstrain
(
const PackedBoolList& isMasterEdge,
const indirectPrimitivePatch& pp,
const labelList& meshEdges,
const List<pointConstraint>& constraints,
vectorField& disp
) const;
@ -196,6 +198,16 @@ class autoSnapDriver
List<pointConstraint>& patchConstraints
) const;
//- Detect any diagonal attraction. Returns indices in face
// or (-1, -1) if none
labelPair findDiagonalAttraction
(
const indirectPrimitivePatch& pp,
const vectorField& patchAttraction,
const List<pointConstraint>& patchConstraints,
const label faceI
) const;
//- Return hit if on multiple points
pointIndexHit findMultiPatchPoint
(
@ -371,7 +383,9 @@ class autoSnapDriver
const scalar featureAttract,
const scalarField& snapDist,
const vectorField& nearestDisp,
motionSmoother& meshMover
motionSmoother& meshMover,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;

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

@ -38,50 +38,18 @@ License
#include "treeDataPoint.H"
#include "indexedOctree.H"
#include "snapParameters.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
class listTransform
{
public:
void operator()
(
const vectorTensorTransform& vt,
const bool forward,
List<List<point> >& fld
) const
{
const tensor T
(
forward
? vt.R()
: vt.R().T()
);
forAll(fld, i)
{
List<point>& elems = fld[i];
forAll(elems, elemI)
{
elems[elemI] = transform(T, elems[elemI]);
}
}
}
};
template<class T>
class listPlusEqOp
{
public:
void operator()
(
List<T>& x,
const List<T>& y
) const
void operator()(List<T>& x, const List<T>& y) const
{
label sz = x.size();
x.setSize(sz+y.size());
@ -159,7 +127,9 @@ bool Foam::autoSnapDriver::isFeaturePoint
void Foam::autoSnapDriver::smoothAndConstrain
(
const PackedBoolList& isMasterEdge,
const indirectPrimitivePatch& pp,
const labelList& meshEdges,
const List<pointConstraint>& constraints,
vectorField& disp
) const
@ -197,11 +167,16 @@ void Foam::autoSnapDriver::smoothAndConstrain
{
forAll(pEdges, i)
{
label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
if (constraints[nbrPointI].first() >= nConstraints)
label edgeI = pEdges[i];
if (isMasterEdge[meshEdges[edgeI]])
{
dispSum[pointI] += disp[nbrPointI];
dispCount[pointI]++;
label nbrPointI = edges[pEdges[i]].otherVertex(pointI);
if (constraints[nbrPointI].first() >= nConstraints)
{
dispSum[pointI] += disp[nbrPointI];
dispCount[pointI]++;
}
}
}
}
@ -564,6 +539,9 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
{
const fvMesh& mesh = meshRefiner_.mesh();
const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
// For now just get all surrounding face data. Expensive - should just
// store and sync data on coupled points only
// (see e.g PatchToolsNormals.C)
@ -583,7 +561,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
forAll(pFaces, i)
{
label faceI = pFaces[i];
if (faceSurfaceGlobalRegion[faceI] != -1)
if (isMasterFace[faceI] && faceSurfaceGlobalRegion[faceI] != -1)
{
nFaces++;
}
@ -605,7 +583,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
label faceI = pFaces[i];
label globalRegionI = faceSurfaceGlobalRegion[faceI];
if (globalRegionI != -1)
if (isMasterFace[faceI] && globalRegionI != -1)
{
pNormals[nFaces] = faceSurfaceNormal[faceI];
pDisp[nFaces] = faceDisp[faceI];
@ -694,7 +672,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
pointFaceSurfNormals,
listPlusEqOp<point>(),
List<point>(),
listTransform()
mapDistribute::transform()
);
syncTools::syncPointList
(
@ -703,7 +681,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
pointFaceDisp,
listPlusEqOp<point>(),
List<point>(),
listTransform()
mapDistribute::transform()
);
syncTools::syncPointList
(
@ -712,7 +690,7 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
pointFaceCentres,
listPlusEqOp<point>(),
List<point>(),
listTransform()
mapDistribute::transformPosition()
);
syncTools::syncPointList
(
@ -722,6 +700,25 @@ void Foam::autoSnapDriver::calcNearestFacePointProperties
listPlusEqOp<label>(),
List<label>()
);
// Sort the data according to the face centres. This is only so we get
// consistent behaviour serial and parallel.
labelList visitOrder;
forAll(pointFaceDisp, pointI)
{
List<point>& pNormals = pointFaceSurfNormals[pointI];
List<point>& pDisp = pointFaceDisp[pointI];
List<point>& pFc = pointFaceCentres[pointI];
labelList& pFid = pointFacePatchID[pointI];
sortedOrder(mag(pFc)(), visitOrder);
pNormals = List<point>(pNormals, visitOrder);
pDisp = List<point>(pDisp, visitOrder);
pFc = List<point>(pFc, visitOrder);
pFid = UIndirectList<label>(pFid, visitOrder);
}
}
@ -1360,6 +1357,70 @@ void Foam::autoSnapDriver::stringFeatureEdges
}
// If only two attractions and across face return the face indices
Foam::labelPair Foam::autoSnapDriver::findDiagonalAttraction
(
const indirectPrimitivePatch& pp,
const vectorField& patchAttraction,
const List<pointConstraint>& patchConstraints,
const label faceI
) const
{
const face& f = pp.localFaces()[faceI];
// For now just detect any attraction. Improve this to look at
// actual attraction position and orientation
labelPair attractIndices(-1, -1);
if (f.size() >= 4)
{
forAll(f, fp)
{
label pointI = f[fp];
if (patchConstraints[pointI].first() >= 2)
{
// Attract to feature edge or point
if (attractIndices[0] == -1)
{
// First attraction. Store
attractIndices[0] = fp;
}
else if (attractIndices[1] == -1)
{
// Second attraction. Check if not consecutive to first
// attraction
label fp0 = attractIndices[0];
if (f.fcIndex(fp0) == fp || f.fcIndex(fp) == fp0)
{
// Consecutive. Skip.
attractIndices = labelPair(-1, -1);
break;
}
else
{
attractIndices[1] = fp;
}
}
else
{
// More than two attractions. Skip.
attractIndices = labelPair(-1, -1);
break;
}
}
}
if (attractIndices[1] == -1)
{
// Found only one attraction. Skip.
attractIndices = labelPair(-1, -1);
}
}
return attractIndices;
}
Foam::pointIndexHit Foam::autoSnapDriver::findNearFeatureEdge
(
const indirectPrimitivePatch& pp,
@ -1929,7 +1990,7 @@ void Foam::autoSnapDriver::featureAttractionUsingFeatureEdges
edgeFaceNormals,
listPlusEqOp<point>(),
List<point>(),
listTransform()
mapDistribute::transform()
);
}
@ -2778,7 +2839,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
const scalar featureAttract,
const scalarField& snapDist,
const vectorField& nearestDisp,
motionSmoother& meshMover
motionSmoother& meshMover,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const
{
const Switch implicitFeatureAttraction = snapParams.implicitFeatureSnap();
@ -2851,7 +2914,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// - faceSurfaceNormal
// - faceDisp
// - faceCentres&faceNormal
// - faceCentres
List<List<point> > pointFaceSurfNormals(pp.nPoints());
List<List<point> > pointFaceDisp(pp.nPoints());
List<List<point> > pointFaceCentres(pp.nPoints());
@ -2884,10 +2947,11 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// here.
// Nearest feature
vectorField patchAttraction(localPoints.size(), vector::zero);
patchAttraction.setSize(localPoints.size());
patchAttraction = vector::zero;
// Constraints at feature
List<pointConstraint> patchConstraints(localPoints.size());
patchConstraints.setSize(localPoints.size());
patchConstraints = pointConstraint();
if (implicitFeatureAttraction)
{
@ -2951,15 +3015,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
patchConstraints
);
const PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
{
vector avgPatchDisp = meshRefinement::gAverage
(
mesh,
isMasterPoint,
pp.meshPoints(),
patchDisp
);
vector avgPatchAttr = meshRefinement::gAverage
(
mesh,
isMasterPoint,
pp.meshPoints(),
patchAttraction
);
Info<< "Attraction:" << endl
<< " linear : max:" << gMax(patchDisp)
<< " avg:" << gAverage(patchDisp)
<< endl
<< " feature : max:" << gMax(patchAttraction)
<< " avg:" << gAverage(patchAttraction)
<< endl;
Info<< "Attraction:" << endl
<< " linear : max:" << gMaxMagSqr(patchDisp)
<< " avg:" << avgPatchDisp << endl
<< " feature : max:" << gMaxMagSqr(patchAttraction)
<< " avg:" << avgPatchAttr << endl;
}
// So now we have:
// - patchDisp : point movement to go to nearest point on surface
@ -2986,37 +3064,45 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// Count
{
const labelList& meshPoints = pp.meshPoints();
label nMasterPoints = 0;
label nPlanar = 0;
label nEdge = 0;
label nPoint = 0;
forAll(patchConstraints, pointI)
{
if (patchConstraints[pointI].first() == 1)
if (isMasterPoint[meshPoints[pointI]])
{
nPlanar++;
}
else if (patchConstraints[pointI].first() == 2)
{
nEdge++;
}
else if (patchConstraints[pointI].first() == 3)
{
nPoint++;
nMasterPoints++;
if (patchConstraints[pointI].first() == 1)
{
nPlanar++;
}
else if (patchConstraints[pointI].first() == 2)
{
nEdge++;
}
else if (patchConstraints[pointI].first() == 3)
{
nPoint++;
}
}
}
label nTotPoints = returnReduce(pp.nPoints(), sumOp<label>());
reduce(nMasterPoints, sumOp<label>());
reduce(nPlanar, sumOp<label>());
reduce(nEdge, sumOp<label>());
reduce(nPoint, sumOp<label>());
Info<< "Feature analysis : total points:"
<< nTotPoints
Info<< "Feature analysis : total master points:"
<< nMasterPoints
<< " attraction to :" << nl
<< " feature point : " << nPoint << nl
<< " feature edge : " << nEdge << nl
<< " nearest surface : " << nPlanar << nl
<< " rest : " << nTotPoints-nPoint-nEdge-nPlanar
<< " rest : " << nMasterPoints-nPoint-nEdge-nPlanar
<< nl
<< endl;
}
@ -3035,11 +3121,29 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
if (featureAttract < 1-0.001)
{
const PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
const vectorField pointNormals
(
PatchTools::pointNormals
(
mesh,
pp
)
);
const labelList meshEdges
(
pp.meshEdges(mesh.edges(), mesh.pointEdges())
);
// 1. Smoothed all displacement
vectorField smoothedPatchDisp = patchDisp;
smoothAndConstrain
(
isMasterEdge,
pp,
meshEdges,
patchConstraints,
smoothedPatchDisp
);
@ -3047,16 +3151,18 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
// 2. Smoothed tangential component
vectorField tangPatchDisp = patchDisp;
tangPatchDisp -= (pp.pointNormals() & patchDisp) * pp.pointNormals();
tangPatchDisp -= (pointNormals & patchDisp) * pointNormals;
smoothAndConstrain
(
isMasterEdge,
pp,
meshEdges,
patchConstraints,
tangPatchDisp
);
// Re-add normal component
tangPatchDisp += (pp.pointNormals() & patchDisp) * pp.pointNormals();
tangPatchDisp += (pointNormals & patchDisp) * pointNormals;
if (debug&meshRefinement::OBJINTERSECTIONS)
{

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

@ -259,23 +259,14 @@ private:
label& nRefine
) const;
//- Mark cell if local curvature > curvature or
// markDifferingRegions = true and intersections with different
// regions.
bool checkCurvature
//- Helper: count number of normals1 that are in normals2
label countMatches
(
const scalar curvature,
const label nAllowRefine,
const label surfaceLevel,
const vector& surfaceNormal,
const label cellI,
label& cellMaxLevel,
vector& cellMaxNormal,
labelList& refineCell,
label& nRefine
const List<point>& normals1,
const List<point>& normals2,
const scalar tol = 1e-6
) const;
//- Mark cells for surface curvature based refinement. Marks if
// local curvature > curvature or if on different regions
// (markDifferingRegions)

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

@ -38,7 +38,74 @@ License
#include "featureEdgeMesh.H"
#include "Cloud.H"
//#include "globalIndex.H"
//#include "OBJstream.H"
#include "OBJstream.H"
#include "cellSet.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
//- To compare normals
static bool less(const vector& x, const vector& y)
{
for (direction i = 0; i < vector::nComponents; i++)
{
if (x[i] < y[i])
{
return true;
}
else if (x[i] > y[i])
{
return false;
}
}
// All components the same
return false;
}
//- To compare normals
class normalLess
{
const vectorList& values_;
public:
normalLess(const vectorList& values)
:
values_(values)
{}
bool operator()(const label a, const label b) const
{
return less(values_[a], values_[b]);
}
};
//- template specialization for pTraits<labelList> so we can have fields
template<>
class pTraits<labelList>
{
public:
//- Component type
typedef labelList cmptType;
};
//- template specialization for pTraits<labelList> so we can have fields
template<>
class pTraits<vectorList>
{
public:
//- Component type
typedef vectorList cmptType;
};
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -954,64 +1021,33 @@ Foam::label Foam::meshRefinement::markSurfaceRefinement
}
// Checks if multiple intersections of a cell (by a surface with a higher
// max than the cell level) and if so if the normals at these intersections
// make a large angle.
// Returns false if the nRefine limit has been reached, true otherwise.
bool Foam::meshRefinement::checkCurvature
// Count number of matches of first argument in second argument
Foam::label Foam::meshRefinement::countMatches
(
const scalar curvature,
const label nAllowRefine,
const label surfaceLevel, // current intersection max level
const vector& surfaceNormal,// current intersection normal
const label cellI,
label& cellMaxLevel, // cached max surface level for this cell
vector& cellMaxNormal, // cached surface normal for this cell
labelList& refineCell,
label& nRefine
const List<point>& normals1,
const List<point>& normals2,
const scalar tol
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
label nMatches = 0;
// Test if surface applicable
if (surfaceLevel > cellLevel[cellI])
forAll(normals1, i)
{
if (cellMaxLevel == -1)
{
// First visit of cell. Store
cellMaxLevel = surfaceLevel;
cellMaxNormal = surfaceNormal;
}
else
{
// Second or more visit. Check curvature.
if ((cellMaxNormal & surfaceNormal) < curvature)
{
return markForRefine
(
surfaceLevel, // mark with any non-neg number.
nAllowRefine,
refineCell[cellI],
nRefine
);
}
const vector& n1 = normals1[i];
// Set normal to that of highest surface. Not really necessary
// over here but we reuse cellMax info when doing coupled faces.
if (surfaceLevel > cellMaxLevel)
forAll(normals2, j)
{
const vector& n2 = normals2[j];
if (magSqr(n1-n2) < tol)
{
cellMaxLevel = surfaceLevel;
cellMaxNormal = surfaceNormal;
nMatches++;
break;
}
}
}
// Did not reach refinement limit.
return true;
return nMatches;
}
@ -1039,6 +1075,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
// on a different surface gets refined (if its current level etc.)
const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
// Collect candidate faces (i.e. intersecting any surface and
// owner/neighbour not yet refined.
labelList testFaces(getRefineCandidateFaces(refineCell));
@ -1068,6 +1107,12 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
start[i] = cellCentres[own];
end[i] = neiCc[bFaceI];
if (!isMasterFace[faceI])
{
Swap(start[i], end[i]);
}
minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
}
}
@ -1081,10 +1126,9 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
// Test for all intersections (with surfaces of higher max level than
// minLevel) and cache per cell the max surface level and the local normal
// on that surface.
labelList cellMaxLevel(mesh_.nCells(), -1);
vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
// minLevel) and cache per cell the interesting inter
labelListList cellSurfLevels(mesh_.nCells());
List<vectorList> cellSurfNormals(mesh_.nCells());
{
// Per segment the normals of the surfaces
@ -1104,12 +1148,29 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
surfaceNormal,
surfaceLevel
);
// Sort the data according to surface location. This will guarantee
// that on coupled faces both sides visit the intersections in
// the same order so will decide the same
labelList visitOrder;
forAll(surfaceNormal, pointI)
{
vectorList& pNormals = surfaceNormal[pointI];
labelList& pLevel = surfaceLevel[pointI];
sortedOrder(pNormals, visitOrder, normalLess(pNormals));
pNormals = List<point>(pNormals, visitOrder);
pLevel = UIndirectList<label>(pLevel, visitOrder);
}
// Clear out unnecessary data
start.clear();
end.clear();
minLevel.clear();
// Extract per cell information on the surface with the highest max
// Convert face-wise data to cell.
forAll(testFaces, i)
{
label faceI = testFaces[i];
@ -1118,164 +1179,280 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
const vectorList& fNormals = surfaceNormal[i];
const labelList& fLevels = surfaceLevel[i];
forAll(fLevels, hitI)
forAll(fNormals, hitI)
{
checkCurvature
(
curvature,
nAllowRefine,
fLevels[hitI],
fNormals[hitI],
own,
cellMaxLevel[own],
cellMaxNormal[own],
refineCell,
nRefine
);
}
if (mesh_.isInternalFace(faceI))
{
label nei = mesh_.faceNeighbour()[faceI];
forAll(fLevels, hitI)
if (fLevels[hitI] > cellLevel[own])
{
checkCurvature
(
curvature,
nAllowRefine,
cellSurfLevels[own].append(fLevels[hitI]);
cellSurfNormals[own].append(fNormals[hitI]);
}
fLevels[hitI],
fNormals[hitI],
nei,
cellMaxLevel[nei],
cellMaxNormal[nei],
refineCell,
nRefine
);
if (mesh_.isInternalFace(faceI))
{
label nei = mesh_.faceNeighbour()[faceI];
if (fLevels[hitI] > cellLevel[nei])
{
cellSurfLevels[nei].append(fLevels[hitI]);
cellSurfNormals[nei].append(fNormals[hitI]);
}
}
}
}
}
// 2. Find out a measure of surface curvature and region edges.
// Send over surface region and surface normal to neighbour cell.
labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
// Bit of statistics
{
label bFaceI = faceI-mesh_.nInternalFaces();
label own = mesh_.faceOwner()[faceI];
neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
label nSet = 0;
label nNormals = 9;
forAll(cellSurfNormals, cellI)
{
const vectorList& normals = cellSurfNormals[cellI];
if (normals.size())
{
nSet++;
nNormals += normals.size();
}
}
reduce(nSet, sumOp<label>());
reduce(nNormals, sumOp<label>());
Info<< "markSurfaceCurvatureRefinement :"
<< " cells:" << mesh_.globalData().nTotalCells()
<< " of which with normals:" << nSet
<< " ; total normals stored:" << nNormals
<< endl;
}
syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
// Loop over all faces. Could only be checkFaces.. except if they're coupled
bool reachedLimit = false;
// 1. Check normals per cell
// ~~~~~~~~~~~~~~~~~~~~~~~~~
for
(
label cellI = 0;
!reachedLimit && cellI < cellSurfNormals.size();
cellI++
)
{
const vectorList& normals = cellSurfNormals[cellI];
const labelList& levels = cellSurfLevels[cellI];
// n^2 comparison of all normals in a cell
for (label i = 0; !reachedLimit && i < normals.size(); i++)
{
for (label j = i+1; !reachedLimit && j < normals.size(); j++)
{
if ((normals[i] & normals[j]) < curvature)
{
label maxLevel = max(levels[i], levels[j]);
if (cellLevel[cellI] < maxLevel)
{
if
(
!markForRefine
(
maxLevel,
nAllowRefine,
refineCell[cellI],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching my cell"
<< " limit of " << mesh_.nCells()+7*nRefine
<< endl;
}
reachedLimit = true;
break;
}
}
}
}
}
}
// 2. Find out a measure of surface curvature
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Look at normals between neighbouring surfaces
// Loop over all faces. Could only be checkFaces, except if they're coupled
// Internal faces
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
for
(
label faceI = 0;
!reachedLimit && faceI < mesh_.nInternalFaces();
faceI++
)
{
label own = mesh_.faceOwner()[faceI];
label nei = mesh_.faceNeighbour()[faceI];
if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
{
// Have valid data on both sides. Check curvature.
if ((cellMaxNormal[own] & cellMaxNormal[nei]) < curvature)
{
// See which side to refine
if (cellLevel[own] < cellMaxLevel[own])
{
if
(
!markForRefine
(
cellMaxLevel[own],
nAllowRefine,
refineCell[own],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching my cell"
<< " limit of " << mesh_.nCells()+7*nRefine
<< endl;
}
break;
}
}
const vectorList& ownNormals = cellSurfNormals[own];
const labelList& ownLevels = cellSurfLevels[own];
const vectorList& neiNormals = cellSurfNormals[nei];
const labelList& neiLevels = cellSurfLevels[nei];
if (cellLevel[nei] < cellMaxLevel[nei])
// Special case: owner normals set is a subset of the neighbour
// normals. Do not do curvature refinement since other cell's normals
// do not add any information. Avoids weird corner extensions of
// refinement regions.
bool ownIsSubset =
countMatches(ownNormals, neiNormals)
== ownNormals.size();
bool neiIsSubset =
countMatches(neiNormals, ownNormals)
== neiNormals.size();
if (!ownIsSubset && !neiIsSubset)
{
// n^2 comparison of between ownNormals and neiNormals
for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
{
for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
{
if
(
!markForRefine
(
cellMaxLevel[nei],
nAllowRefine,
refineCell[nei],
nRefine
)
)
// Have valid data on both sides. Check curvature.
if ((ownNormals[i] & neiNormals[j]) < curvature)
{
if (debug)
// See which side to refine.
if (cellLevel[own] < ownLevels[i])
{
Pout<< "Stopped refining since reaching my cell"
<< " limit of " << mesh_.nCells()+7*nRefine
<< endl;
if
(
!markForRefine
(
ownLevels[i],
nAllowRefine,
refineCell[own],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching"
<< " my cell limit of "
<< mesh_.nCells()+7*nRefine << endl;
}
reachedLimit = true;
break;
}
}
if (cellLevel[nei] < neiLevels[j])
{
if
(
!markForRefine
(
neiLevels[j],
nAllowRefine,
refineCell[nei],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching"
<< " my cell limit of "
<< mesh_.nCells()+7*nRefine << endl;
}
reachedLimit = true;
break;
}
}
break;
}
}
}
}
}
// Send over surface normal to neighbour cell.
List<vectorList> neiSurfaceNormals;
syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
// Boundary faces
for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
for
(
label faceI = mesh_.nInternalFaces();
!reachedLimit && faceI < mesh_.nFaces();
faceI++
)
{
label own = mesh_.faceOwner()[faceI];
label bFaceI = faceI - mesh_.nInternalFaces();
if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
const vectorList& ownNormals = cellSurfNormals[own];
const labelList& ownLevels = cellSurfLevels[own];
const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
// Special case: owner normals set is a subset of the neighbour
// normals. Do not do curvature refinement since other cell's normals
// do not add any information. Avoids weird corner extensions of
// refinement regions.
bool ownIsSubset =
countMatches(ownNormals, neiNormals)
== ownNormals.size();
bool neiIsSubset =
countMatches(neiNormals, ownNormals)
== neiNormals.size();
if (!ownIsSubset && !neiIsSubset)
{
// Have valid data on both sides. Check curvature.
if ((cellMaxNormal[own] & neiBndMaxNormal[bFaceI]) < curvature)
// n^2 comparison of between ownNormals and neiNormals
for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
{
if
(
!markForRefine
(
cellMaxLevel[own],
nAllowRefine,
refineCell[own],
nRefine
)
)
for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
{
if (debug)
// Have valid data on both sides. Check curvature.
if ((ownNormals[i] & neiNormals[j]) < curvature)
{
Pout<< "Stopped refining since reaching my cell"
<< " limit of " << mesh_.nCells()+7*nRefine
<< endl;
if (cellLevel[own] < ownLevels[i])
{
if
(
!markForRefine
(
ownLevels[i],
nAllowRefine,
refineCell[own],
nRefine
)
)
{
if (debug)
{
Pout<< "Stopped refining since reaching"
<< " my cell limit of "
<< mesh_.nCells()+7*nRefine
<< endl;
}
reachedLimit = true;
break;
}
}
}
break;
}
}
}
}
if
(
returnReduce(nRefine, sumOp<label>())