From 71c0a5d1d7c8205bc948e66800db8743e72910a9 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 18 Oct 2013 11:52:38 +0100
Subject: [PATCH] ENH: autoHexMesh: parallel consistency, fix curvature
refinement
---
.../autoHexMeshDriver/autoSnapDriver.C | 126 ++++-
.../autoHexMeshDriver/autoSnapDriver.H | 16 +-
.../autoHexMeshDriver/autoSnapDriverFeature.C | 252 ++++++---
.../meshRefinement/meshRefinement.H | 19 +-
.../meshRefinement/meshRefinementRefine.C | 495 ++++++++++++------
5 files changed, 644 insertions(+), 264 deletions(-)
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 07ff697df6a..adb7c385903 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -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
);
}
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
index 77327c0689b..b8282bedfcd 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H
@@ -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;
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
index 887674ea7ae..32133ad5f84 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
@@ -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)
{
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 8a5f327b410..873167bc4f0 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -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)
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
index 7444493e298..ce84ceec38e 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
@@ -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
+ const vector& n1 = normals1[i];
+
+ forAll(normals2, j)
{
- // 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& n2 = normals2[j];
- // Set normal to that of highest surface. Not really necessary
- // over here but we reuse cellMax info when doing coupled faces.
- if (surfaceLevel > cellMaxLevel)
+ 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],
+ if (fLevels[hitI] > cellLevel[own])
+ {
+ cellSurfLevels[own].append(fLevels[hitI]);
+ cellSurfNormals[own].append(fNormals[hitI]);
+ }
- 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]);
+ }
+ }
}
+ }
+ }
- if (mesh_.isInternalFace(faceI))
- {
- label nei = mesh_.faceNeighbour()[faceI];
-
- forAll(fLevels, hitI)
- {
- checkCurvature
- (
- curvature,
- nAllowRefine,
- fLevels[hitI],
- fNormals[hitI],
- nei,
- cellMaxLevel[nei],
- cellMaxNormal[nei],
-
- refineCell,
- nRefine
- );
- }
+ // Bit of statistics
+ {
+ 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;
}
- // 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++)
- {
- label bFaceI = faceI-mesh_.nInternalFaces();
- label own = mesh_.faceOwner()[faceI];
+ bool reachedLimit = false;
- neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
- neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
- }
- syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
- syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
- // Loop over all faces. Could only be checkFaces.. except if they're coupled
+ // 1. Check normals per cell
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~
- // Internal faces
- for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+ for
+ (
+ label cellI = 0;
+ !reachedLimit && cellI < cellSurfNormals.size();
+ cellI++
+ )
{
- label own = mesh_.faceOwner()[faceI];
- label nei = mesh_.faceNeighbour()[faceI];
+ const vectorList& normals = cellSurfNormals[cellI];
+ const labelList& levels = cellSurfLevels[cellI];
- if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
+ // n^2 comparison of all normals in a cell
+ for (label i = 0; !reachedLimit && i < normals.size(); i++)
{
- // Have valid data on both sides. Check curvature.
- if ((cellMaxNormal[own] & cellMaxNormal[nei]) < curvature)
+ for (label j = i+1; !reachedLimit && j < normals.size(); j++)
{
- // See which side to refine
- if (cellLevel[own] < cellMaxLevel[own])
+ if ((normals[i] & normals[j]) < curvature)
{
- if
- (
- !markForRefine
+ label maxLevel = max(levels[i], levels[j]);
+
+ if (cellLevel[cellI] < maxLevel)
+ {
+ if
(
- cellMaxLevel[own],
- nAllowRefine,
- refineCell[own],
- nRefine
+ !markForRefine
+ (
+ maxLevel,
+ nAllowRefine,
+ refineCell[cellI],
+ nRefine
+ )
)
- )
- {
- if (debug)
{
- Pout<< "Stopped refining since reaching my cell"
- << " limit of " << mesh_.nCells()+7*nRefine
- << endl;
+ if (debug)
+ {
+ Pout<< "Stopped refining since reaching my cell"
+ << " limit of " << mesh_.nCells()+7*nRefine
+ << endl;
+ }
+ reachedLimit = true;
+ break;
}
- break;
}
}
+ }
+ }
+ }
- if (cellLevel[nei] < cellMaxLevel[nei])
+
+
+ // 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;
+ !reachedLimit && faceI < mesh_.nInternalFaces();
+ faceI++
+ )
+ {
+ label own = mesh_.faceOwner()[faceI];
+ label nei = mesh_.faceNeighbour()[faceI];
+
+ const vectorList& ownNormals = cellSurfNormals[own];
+ const labelList& ownLevels = cellSurfLevels[own];
+ const vectorList& neiNormals = cellSurfNormals[nei];
+ const labelList& neiLevels = cellSurfLevels[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>())
--
GitLab