From 868a090243a011484384203f6038db4f4ce3d861 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 16 Jan 2013 17:08:26 +0000
Subject: [PATCH] BUG: snappyHexMesh: synchronisation when patch edges align
with processor edges
---
.../autoHexMeshDriver/autoLayerDriver.C | 42 ++++-
.../autoHexMeshDriver/autoLayerDriverShrink.C | 110 +++++++++++++-
.../algorithms/PointEdgeWave/PointEdgeWave.C | 143 ++++++++++++------
3 files changed, 245 insertions(+), 50 deletions(-)
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index d0e43dbb10e..b224ef609d8 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -293,6 +293,46 @@ void Foam::autoLayerDriver::handleNonManifolds
}
}
+ // 3. Remote check for end of layer across coupled boundaries
+ {
+ PackedBoolList isCoupledEdge(mesh.nEdges());
+
+ const labelList& cpEdges = mesh.globalData().coupledPatchMeshEdges();
+ forAll(cpEdges, i)
+ {
+ isCoupledEdge[cpEdges[i]] = true;
+ }
+ syncTools::syncEdgeList
+ (
+ mesh,
+ isCoupledEdge,
+ orEqOp<unsigned int>(),
+ 0
+ );
+
+ forAll(edgeGlobalFaces, edgeI)
+ {
+ label meshEdgeI = meshEdges[edgeI];
+
+ if
+ (
+ pp.edgeFaces()[edgeI].size() == 1
+ && edgeGlobalFaces[edgeI].size() == 1
+ && isCoupledEdge[meshEdgeI]
+ )
+ {
+ // Edge of patch but no continuation across processor.
+ const edge& e = pp.edges()[edgeI];
+ //Pout<< "** Stopping extrusion on edge "
+ // << pp.localPoints()[e[0]]
+ // << pp.localPoints()[e[1]] << endl;
+ nonManifoldPoints.insert(pp.meshPoints()[e[0]]);
+ nonManifoldPoints.insert(pp.meshPoints()[e[1]]);
+ }
+ }
+ }
+
+
label nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
index 8ffcaa5f578..bf25526141c 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -351,6 +351,11 @@ void Foam::autoLayerDriver::smoothNormals
isFixedPoint.set(meshPointI, 1);
}
+ // Make sure that points that are coupled to meshPoints but not on a patch
+ // are fixed as well
+ syncTools::syncPointList(mesh, isFixedPoint, maxEqOp<unsigned int>(), 0);
+
+
// Correspondence between local edges/points and mesh edges/points
const labelList meshEdges(identity(mesh.nEdges()));
const labelList meshPoints(identity(mesh.nPoints()));
@@ -831,6 +836,19 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
)
);
+ // pointNormals
+ if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
+ {
+ pointField meshPointNormals(mesh.nPoints(), point(1, 0, 0));
+ UIndirectList<point>(meshPointNormals, pp.meshPoints()) = pointNormals;
+ meshRefinement::testSyncPointList
+ (
+ "pointNormals",
+ mesh,
+ meshPointNormals
+ );
+ }
+
// Smooth patch normal vectors
smoothPatchNormals
(
@@ -841,6 +859,18 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
pointNormals
);
+ // smoothed pointNormals
+ if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
+ {
+ pointField meshPointNormals(mesh.nPoints(), point(1, 0, 0));
+ UIndirectList<point>(meshPointNormals, pp.meshPoints()) = pointNormals;
+ meshRefinement::testSyncPointList
+ (
+ "smoothed pointNormals",
+ mesh,
+ meshPointNormals
+ );
+ }
// Calculate distance to pp points
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -883,6 +913,28 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
);
}
+
+ // Check sync of wall distance
+ if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
+ {
+ pointField origin(pointWallDist.size());
+ scalarField distSqr(pointWallDist.size());
+ scalarField passiveS(pointWallDist.size());
+ pointField passiveV(pointWallDist.size());
+ forAll(pointWallDist, pointI)
+ {
+ origin[pointI] = pointWallDist[pointI].origin();
+ distSqr[pointI] = pointWallDist[pointI].distSqr();
+ passiveS[pointI] = pointWallDist[pointI].s();
+ passiveV[pointI] = pointWallDist[pointI].v();
+ }
+ meshRefinement::testSyncPointList("origin", mesh, origin);
+ meshRefinement::testSyncPointList("distSqr", mesh, distSqr);
+ meshRefinement::testSyncPointList("passiveS", mesh, passiveS);
+ meshRefinement::testSyncPointList("passiveV", mesh, passiveV);
+ }
+
+
// 2. Find points with max distance and transport information back to
// wall.
{
@@ -1095,6 +1147,13 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr());
medialVec[pointI] = pointMedialDist[pointI].origin();
}
+
+ // Check
+ if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
+ {
+ meshRefinement::testSyncPointList("medialDist", mesh, medialDist);
+ meshRefinement::testSyncPointList("medialVec", mesh, medialVec);
+ }
}
// Extract transported surface normals as pointVectorField
@@ -1106,6 +1165,11 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
// Smooth normal vectors. Do not change normals on pp.meshPoints
smoothNormals(nSmoothNormals, isMasterEdge, meshPoints, dispVec);
+ if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
+ {
+ meshRefinement::testSyncPointList("smoothed dispVec", mesh, dispVec);
+ }
+
// Calculate ratio point medial distance to point wall distance
forAll(medialRatio, pointI)
{
@@ -1122,6 +1186,14 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
}
}
+
+ if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
+ {
+ // medialRatio
+ meshRefinement::testSyncPointList("medialRatio", mesh, medialRatio);
+ }
+
+
if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
{
Info<< "medialAxisSmoothingInfo :"
@@ -1417,6 +1489,42 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
+
+ // Check a bit the sync of displacements
+ if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
+ {
+ // initial mesh
+ meshRefinement::testSyncPointList("mesh.points()", mesh, mesh.points());
+
+ // pointWallDist
+ scalarField pWallDist(pointWallDist.size());
+ forAll(pointWallDist, pointI)
+ {
+ pWallDist[pointI] = pointWallDist[pointI].s();
+ }
+ meshRefinement::testSyncPointList("pointWallDist", mesh, pWallDist);
+
+ // dispVec
+ meshRefinement::testSyncPointList("dispVec", mesh, dispVec);
+
+ // displacement before and after correction
+ meshRefinement::testSyncPointList
+ (
+ "displacement BEFORE",
+ mesh,
+ displacement
+ );
+
+ meshMover.correctBoundaryConditions(displacement);
+ meshRefinement::testSyncPointList
+ (
+ "displacement AFTER",
+ mesh,
+ displacement
+ );
+ }
+
+
//XXXXX
// // Smear displacement away from fixed values (medialRatio=0 or 1)
// {
diff --git a/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C b/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C
index db337fc098f..ad57079dbcc 100644
--- a/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C
+++ b/src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C
@@ -43,6 +43,30 @@ Foam::scalar Foam::PointEdgeWave<Type, TrackingData>::propagationTol_ = 0.01;
template <class Type, class TrackingData>
int Foam::PointEdgeWave<Type, TrackingData>::dummyTrackData_ = 12345;
+namespace Foam
+{
+ //- Reduction class. If x and y are not equal assign value.
+ template<class Type, class TrackingData>
+ class combineEqOp
+ {
+ TrackingData& td_;
+
+ public:
+ combineEqOp(TrackingData& td)
+ :
+ td_(td)
+ {}
+
+ void operator()(Type& x, const Type& y) const
+ {
+ if (!x.valid(td_) && y.valid(td_))
+ {
+ x = y;
+ }
+ }
+ };
+}
+
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -325,12 +349,12 @@ void Foam::PointEdgeWave<Type, TrackingData>::handleProcPatches()
// Adapt for leaving domain
leaveDomain(procPatch, thisPoints, patchInfo);
- if (debug)
- {
- Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
- << " communicating with " << procPatch.neighbProcNo()
- << " Sending:" << patchInfo.size() << endl;
- }
+ //if (debug)
+ //{
+ // Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
+ // << " communicating with " << procPatch.neighbProcNo()
+ // << " Sending:" << patchInfo.size() << endl;
+ //}
UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
toNeighbour << nbrPoints << patchInfo;
@@ -357,12 +381,12 @@ void Foam::PointEdgeWave<Type, TrackingData>::handleProcPatches()
fromNeighbour >> patchPoints >> patchInfo;
}
- if (debug)
- {
- Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
- << " communicating with " << procPatch.neighbProcNo()
- << " Received:" << patchInfo.size() << endl;
- }
+ //if (debug)
+ //{
+ // Pout<< "Processor patch " << patchI << ' ' << procPatch.name()
+ // << " communicating with " << procPatch.neighbProcNo()
+ // << " Received:" << patchInfo.size() << endl;
+ //}
// Apply transform to received data for non-parallel planes
if (!procPatch.parallel())
@@ -492,12 +516,12 @@ void Foam::PointEdgeWave<Type, TrackingData>::handleCyclicPatches()
transform(cycPatch, cycPatch.forwardT(), nbrInfo);
}
- if (debug)
- {
- Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
- << " Changed : " << nbrInfo.size()
- << endl;
- }
+ //if (debug)
+ //{
+ // Pout<< "Cyclic patch " << patchI << ' ' << patch.name()
+ // << " Changed : " << nbrInfo.size()
+ // << endl;
+ //}
// Adapt for entering domain
enterDomain(cycPatch, thisPoints, nbrInfo);
@@ -523,7 +547,8 @@ void Foam::PointEdgeWave<Type, TrackingData>::handleCyclicPatches()
}
-// Propagate information from edge to point. Return number of points changed.
+// Guarantee collocated points have same information.
+// Return number of points changed.
template <class Type, class TrackingData>
Foam::label Foam::PointEdgeWave<Type, TrackingData>::handleCollocatedPoints()
{
@@ -541,30 +566,22 @@ Foam::label Foam::PointEdgeWave<Type, TrackingData>::handleCollocatedPoints()
elems[pointI] = allPointInfo_[meshPoints[pointI]];
}
- // Reset changed points counter.
- nChangedPoints_ = 0;
-
// Pull slave data onto master. No need to update transformed slots.
slavesMap.distribute(elems, false);
// Combine master data with slave data
+ combineEqOp<Type, TrackingData> cop(td_);
+
forAll(slaves, pointI)
{
Type& elem = elems[pointI];
const labelList& slavePoints = slaves[pointI];
- label meshPointI = meshPoints[pointI];
-
// Combine master with untransformed slave data
forAll(slavePoints, j)
{
- updatePoint
- (
- meshPointI,
- elems[slavePoints[j]],
- elem
- );
+ cop(elem, elems[slavePoints[j]]);
}
// Copy result back to slave slots
@@ -580,7 +597,28 @@ Foam::label Foam::PointEdgeWave<Type, TrackingData>::handleCollocatedPoints()
// Extract back onto mesh
forAll(meshPoints, pointI)
{
- allPointInfo_[meshPoints[pointI]] = elems[pointI];
+ Type& elem = allPointInfo_[meshPoints[pointI]];
+
+ // Like updatePoint but bypass Type::updatePoint with its tolerance
+ // checking
+ if (!elem.valid(td_) || !elem.equal(elems[pointI], td_))
+ {
+ nEvals_++;
+ elem = elems[pointI];
+
+ // See if element now valid
+ if (elem.valid(td_))
+ {
+ --nUnvisitedPoints_;
+ }
+
+ // Update database of changed points
+ if (!changedPoint_[pointI])
+ {
+ changedPoint_[pointI] = true;
+ changedPoints_[nChangedPoints_++] = pointI;
+ }
+ }
}
// Sum nChangedPoints over all procs
@@ -658,7 +696,8 @@ Foam::PointEdgeWave<Type, TrackingData>::PointEdgeWave
if (debug)
{
- Pout<< "Seed points : " << nChangedPoints_ << endl;
+ Info<< typeName << ": Seed points : "
+ << returnReduce(nChangedPoints_, sumOp<label>()) << endl;
}
// Iterate until nothing changes
@@ -761,6 +800,9 @@ void Foam::PointEdgeWave<Type, TrackingData>::setPointInfo
changedPoints_[nChangedPoints_++] = pointI;
}
}
+
+ // Sync
+ handleCollocatedPoints();
}
@@ -826,10 +868,10 @@ Foam::label Foam::PointEdgeWave<Type, TrackingData>::edgeToPoint()
handleProcPatches();
}
- if (debug)
- {
- Pout<< "Changed points : " << nChangedPoints_ << endl;
- }
+ //if (debug)
+ //{
+ // Pout<< "Changed points : " << nChangedPoints_ << endl;
+ //}
// Sum nChangedPoints over all procs
label totNChanged = nChangedPoints_;
@@ -894,10 +936,10 @@ Foam::label Foam::PointEdgeWave<Type, TrackingData>::pointToEdge()
// Handled all changed points by now
nChangedPoints_ = 0;
- if (debug)
- {
- Pout<< "Changed edges : " << nChangedEdges_ << endl;
- }
+ //if (debug)
+ //{
+ // Pout<< "Changed edges : " << nChangedEdges_ << endl;
+ //}
// Sum nChangedPoints over all procs
label totNChanged = nChangedEdges_;
@@ -936,14 +978,15 @@ Foam::label Foam::PointEdgeWave<Type, TrackingData>::iterate
{
if (debug)
{
- Pout<< "Iteration " << iter << endl;
+ Info<< typeName << ": Iteration " << iter << endl;
}
label nEdges = pointToEdge();
if (debug)
{
- Pout<< "Total changed edges : " << nEdges << endl;
+ Info<< typeName << ": Total changed edges : "
+ << nEdges << endl;
}
if (nEdges == 0)
@@ -955,10 +998,14 @@ Foam::label Foam::PointEdgeWave<Type, TrackingData>::iterate
if (debug)
{
- Pout<< "Total changed points : " << nPoints << nl
- << "Total evaluations : " << nEvals_ << nl
- << "Remaining unvisited points: " << nUnvisitedPoints_ << nl
- << "Remaining unvisited edges : " << nUnvisitedEdges_ << nl
+ Info<< typeName << ": Total changed points : "
+ << nPoints << nl
+ << typeName << ": Total evaluations : "
+ << returnReduce(nEvals_, sumOp<label>()) << nl
+ << typeName << ": Remaining unvisited points: "
+ << returnReduce(nUnvisitedPoints_, sumOp<label>()) << nl
+ << typeName << ": Remaining unvisited edges : "
+ << returnReduce(nUnvisitedEdges_, sumOp<label>()) << nl
<< endl;
}
@@ -976,8 +1023,8 @@ Foam::label Foam::PointEdgeWave<Type, TrackingData>::iterate
label nPoints = handleCollocatedPoints();
if (debug)
{
- Pout<< "Collocated point sync : " << nPoints << nl
- << endl;
+ Info<< typeName << ": Collocated point sync : "
+ << nPoints << nl << endl;
}
if (nPoints == 0)
--
GitLab