Skip to content

GitLab

Commit 868a0902 authored by mattijs's avatar mattijs
Browse files

BUG: snappyHexMesh: synchronisation when patch edges align with processor edges

parent 376f60b8
Hide whitespace changes
Inline Side-by-side
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View file @ 868a0902
......@@ -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>());
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
View file @ 868a0902
......@@ -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)
// {
......
src/meshTools/algorithms/PointEdgeWave/PointEdgeWave.C
View file @ 868a0902
......@@ -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)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment