zhyang-dev/OpenFOAM-5.x
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

refineMesh: Correct parallel operation

Browse Source 
Patch contributed by Mattijs Janssens
Resolves bug-report https://bugs.openfoam.org/view.php?id=2621
This commit is contained in:
Henry Weller
2017-07-22 22:39:40 +01:00
parent 0c4175bec7
commit c6f365f6cc
6 changed files with 309 additions and 89 deletions
Download Patch File Download Diff File Expand all files Collapse all files

68
applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
View File

@ -40,18 +40,11 @@ Description
#include "argList.H"
#include "polyMesh.H"
#include "Time.H"
#include "undoableMeshCutter.H"
#include "hexCellLooper.H"
#include "cellSet.H"
#include "twoDPointCorrector.H"
#include "directions.H"
#include "OFstream.H"
#include "multiDirRefinement.H"
#include "labelIOList.H"
#include "wedgePolyPatch.H"
#include "plane.H"
#include "SubField.H"
#include "IOdictionary.H"
#include "syncTools.H"
using namespace Foam;
@ -62,7 +55,7 @@ static const scalar edgeTol = 1e-3;
// Print edge statistics on mesh.
void printEdgeStats(const primitiveMesh& mesh)
void printEdgeStats(const polyMesh& mesh)
{
label nX = 0;
label nY = 0;
@ -70,22 +63,26 @@ void printEdgeStats(const primitiveMesh& mesh)
scalar minX = GREAT;
scalar maxX = -GREAT;
vector x(1, 0, 0);
static const vector x(1, 0, 0);
scalar minY = GREAT;
scalar maxY = -GREAT;
vector y(0, 1, 0);
static const vector y(0, 1, 0);
scalar minZ = GREAT;
scalar maxZ = -GREAT;
vector z(0, 0, 1);
static const vector z(0, 0, 1);
scalar minOther = GREAT;
scalar maxOther = -GREAT;
PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
const edgeList& edges = mesh.edges();
forAll(edges, edgeI)
{
if (isMasterEdge[edgeI])
{
const edge& e = edges[edgeI];
@ -119,17 +116,37 @@ void printEdgeStats(const primitiveMesh& mesh)
maxOther = max(maxOther, eMag);
}
}
}
Pout<< "Mesh edge statistics:" << endl
label nEdges = mesh.nEdges();
reduce(nEdges, sumOp<label>());
reduce(nX, sumOp<label>());
reduce(nY, sumOp<label>());
reduce(nZ, sumOp<label>());
reduce(minX, minOp<scalar>());
reduce(maxX, maxOp<scalar>());
reduce(minY, minOp<scalar>());
reduce(maxY, maxOp<scalar>());
reduce(minZ, minOp<scalar>());
reduce(maxZ, maxOp<scalar>());
reduce(minOther, minOp<scalar>());
reduce(maxOther, maxOp<scalar>());
Info<< "Mesh edge statistics:" << nl
<< " x aligned : number:" << nX << "\tminLen:" << minX
<< "\tmaxLen:" << maxX << endl
<< "\tmaxLen:" << maxX << nl
<< " y aligned : number:" << nY << "\tminLen:" << minY
<< "\tmaxLen:" << maxY << endl
<< "\tmaxLen:" << maxY << nl
<< " z aligned : number:" << nZ << "\tminLen:" << minZ
<< "\tmaxLen:" << maxZ << endl
<< " other : number:" << mesh.nEdges() - nX - nY - nZ
<< "\tmaxLen:" << maxZ << nl
<< " other : number:" << nEdges - nX - nY - nZ
<< "\tminLen:" << minOther
<< "\tmaxLen:" << maxOther << endl << endl;
<< "\tmaxLen:" << maxOther << nl << endl;
}
@ -228,7 +245,8 @@ int main(int argc, char *argv[])
cellSet cells(mesh, setName);
Pout<< "Read " << cells.size() << " cells from cellSet "
Info<< "Read " << returnReduce(cells.size(), sumOp<label>())
<< " cells from cellSet "
<< cells.instance()/cells.local()/cells.name()
<< endl << endl;
@ -239,12 +257,7 @@ int main(int argc, char *argv[])
Info<< "Refining all cells" << nl << endl;
// Select all cells
refCells.setSize(mesh.nCells());
forAll(mesh.cells(), celli)
{
refCells[celli] = celli;
}
refCells = identity(mesh.nCells());
if (mesh.nGeometricD() == 3)
{
@ -340,7 +353,9 @@ int main(int argc, char *argv[])
}
}
Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
Info<< "Writing refined cells ("
<< returnReduce(newCells.size(), sumOp<label>())
<< ") to cellSet "
<< newCells.instance()/newCells.local()/newCells.name()
<< endl << endl;
@ -348,7 +363,6 @@ int main(int argc, char *argv[])
//
// Invert cell split to construct map from new to old
//

198
src/dynamicMesh/meshCut/cellCuts/cellCuts.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -33,12 +33,25 @@ License
#include "geomCellLooper.H"
#include "OFstream.H"
#include "plane.H"
#include "syncTools.H"
#include "dummyTransform.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(cellCuts, 0);
//- Template specialization for pTraits<edge> so we can use syncTools
// functionality
template<>
class pTraits<edge>
{
public:
//- Component type
typedef edge cmptType;
};
}
@ -114,6 +127,148 @@ Foam::label Foam::cellCuts::firstUnique
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::cellCuts::syncProc()
{
if (!Pstream::parRun())
{
return;
}
syncTools::syncPointList(mesh(), pointIsCut_, orEqOp<bool>(), false);
syncTools::syncEdgeList(mesh(), edgeIsCut_, orEqOp<bool>(), false);
syncTools::syncEdgeList(mesh(), edgeWeight_, maxEqOp<scalar>(), -GREAT);
{
const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
// Convert faceSplitCut into face-local data: vertex and edge w.r.t.
// vertex 0: (since this is same on both sides)
//
// Sending-side vertex Receiving-side vertex
// 0 0
// 1 3
// 2 2
// 3 1
//
// Sending-side edge Receiving side edge
// 0-1 3-0
// 1-2 2-3
// 2-3 1-2
// 3-0 0-1
//
// Encoding is as index:
// 0 : not set
// >0 : vertex, vertex is index-1
// <0 : edge, edge is -index-1
edgeList relCuts(nBnd, edge(0, 0));
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start()+i;
label bFacei = facei-mesh().nInternalFaces();
const Map<edge>::const_iterator iter =
faceSplitCut_.find(facei);
if (iter != faceSplitCut_.end())
{
const face& f = mesh().faces()[facei];
const labelList& fEdges = mesh().faceEdges()[facei];
const edge& cuts = iter();
forAll(cuts, i)
{
if (isEdge(cuts[i]))
{
label edgei = getEdge(cuts[i]);
label index = findIndex(fEdges, edgei);
relCuts[bFacei][i] = -index-1;
}
else
{
label index = findIndex(f, getVertex(cuts[i]));
relCuts[bFacei][i] = index+1;
}
}
}
}
}
}
// Exchange
syncTools::syncBoundaryFaceList
(
mesh(),
relCuts,
eqOp<edge>(),
dummyTransform()
);
// Convert relCuts back into mesh based data
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start()+i;
label bFacei = facei-mesh().nInternalFaces();
const edge& relCut = relCuts[bFacei];
if (relCut != edge(0, 0))
{
const face& f = mesh().faces()[facei];
const labelList& fEdges = mesh().faceEdges()[facei];
edge absoluteCut(0, 0);
forAll(relCut, i)
{
if (relCut[i] < 0)
{
label oppFp = -relCut[i]-1;
label fp = f.size()-1-oppFp;
absoluteCut[i] = edgeToEVert(fEdges[fp]);
}
else
{
label oppFp = relCut[i]-1;
label fp = f.size()-1-oppFp;
absoluteCut[i] = vertToEVert(f[fp]);
}
}
if
(
!faceSplitCut_.insert(facei, absoluteCut)
&& faceSplitCut_[facei] != absoluteCut
)
{
FatalErrorInFunction
<< "Cut " << faceSplitCut_[facei]
<< " on face " << mesh().faceCentres()[facei]
<< " of coupled patch " << pp.name()
<< " is not consistent with coupled cut "
<< absoluteCut
<< exit(FatalError);
}
}
}
}
}
}
}
void Foam::cellCuts::writeUncutOBJ
(
const fileName& dir,
@ -337,7 +492,7 @@ Foam::label Foam::cellCuts::vertexVertexToFace
void Foam::cellCuts::calcFaceCuts() const
{
if (faceCutsPtr_)
if (faceCutsPtr_.valid())
{
FatalErrorInFunction
<< "faceCuts already calculated" << abort(FatalError);
@ -345,9 +500,8 @@ void Foam::cellCuts::calcFaceCuts() const
const faceList& faces = mesh().faces();
faceCutsPtr_ = new labelListList(mesh().nFaces());
labelListList& faceCuts = *faceCutsPtr_;
faceCutsPtr_.reset(new labelListList(mesh().nFaces()));
labelListList& faceCuts = faceCutsPtr_();
for (label facei = 0; facei < mesh().nFaces(); facei++)
{
@ -2624,6 +2778,16 @@ void Foam::cellCuts::check() const
// Check that cut faces have a neighbour that is cut.
boolList nbrCellIsCut;
{
boolList cellIsCut(mesh().nCells(), false);
forAll(cellLoops_, celli)
{
cellIsCut[celli] = cellLoops_[celli].size();
}
syncTools::swapBoundaryCellList(mesh(), cellIsCut, nbrCellIsCut);
}
forAllConstIter(Map<edge>, faceSplitCut_, iter)
{
label facei = iter.key();
@ -2645,9 +2809,11 @@ void Foam::cellCuts::check() const
}
else
{
label bFacei = facei - mesh().nInternalFaces();
label own = mesh().faceOwner()[facei];
if (cellLoops_[own].empty())
if (cellLoops_[own].empty() && !nbrCellIsCut[bFacei])
{
FatalErrorInFunction
<< "Boundary face:" << facei << " cut by " << iter()
@ -2675,7 +2841,6 @@ Foam::cellCuts::cellCuts
pointIsCut_(expand(mesh.nPoints(), meshVerts)),
edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
faceCutsPtr_(nullptr),
faceSplitCut_(cutCells.size()),
cellLoops_(mesh.nCells()),
nLoops_(-1),
@ -2718,7 +2883,6 @@ Foam::cellCuts::cellCuts
pointIsCut_(expand(mesh.nPoints(), meshVerts)),
edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
faceCutsPtr_(nullptr),
faceSplitCut_(mesh.nFaces()/10 + 1),
cellLoops_(mesh.nCells()),
nLoops_(-1),
@ -2734,6 +2898,9 @@ Foam::cellCuts::cellCuts
calcLoopsAndAddressing(identity(mesh.nCells()));
// Adds cuts on other side of coupled boundaries
syncProc();
// Calculate planes and flip cellLoops if necessary
orientPlanesAndLoops();
@ -2763,7 +2930,6 @@ Foam::cellCuts::cellCuts
pointIsCut_(mesh.nPoints(), false),
edgeIsCut_(mesh.nEdges(), false),
edgeWeight_(mesh.nEdges(), -GREAT),
faceCutsPtr_(nullptr),
faceSplitCut_(cellLabels.size()),
cellLoops_(mesh.nCells()),
nLoops_(-1),
@ -2778,6 +2944,9 @@ Foam::cellCuts::cellCuts
// Makes sure cuts are consistent
setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);
// Adds cuts on other side of coupled boundaries
syncProc();
// Calculate planes and flip cellLoops if necessary
orientPlanesAndLoops();
@ -2806,7 +2975,6 @@ Foam::cellCuts::cellCuts
pointIsCut_(mesh.nPoints(), false),
edgeIsCut_(mesh.nEdges(), false),
edgeWeight_(mesh.nEdges(), -GREAT),
faceCutsPtr_(nullptr),
faceSplitCut_(refCells.size()),
cellLoops_(mesh.nCells()),
nLoops_(-1),
@ -2821,6 +2989,9 @@ Foam::cellCuts::cellCuts
// Makes sure cuts are consistent
setFromCellCutter(cellCutter, refCells);
// Adds cuts on other side of coupled boundaries
syncProc();
// Calculate planes and flip cellLoops if necessary
orientPlanesAndLoops();
@ -2850,7 +3021,6 @@ Foam::cellCuts::cellCuts
pointIsCut_(mesh.nPoints(), false),
edgeIsCut_(mesh.nEdges(), false),
edgeWeight_(mesh.nEdges(), -GREAT),
faceCutsPtr_(nullptr),
faceSplitCut_(cellLabels.size()),
cellLoops_(mesh.nCells()),
nLoops_(-1),
@ -2866,6 +3036,9 @@ Foam::cellCuts::cellCuts
// Makes sure cuts are consistent
setFromCellCutter(cellCutter, cellLabels, cutPlanes);
// Adds cuts on other side of coupled boundaries
syncProc();
// Calculate planes and flip cellLoops if necessary
orientPlanesAndLoops();
@ -2900,7 +3073,6 @@ Foam::cellCuts::cellCuts
pointIsCut_(pointIsCut),
edgeIsCut_(edgeIsCut),
edgeWeight_(edgeWeight),
faceCutsPtr_(nullptr),
faceSplitCut_(faceSplitCut),
cellLoops_(cellLoops),
nLoops_(nLoops),
@ -2924,7 +3096,7 @@ Foam::cellCuts::~cellCuts()
void Foam::cellCuts::clearOut()
{
deleteDemandDrivenData(faceCutsPtr_);
faceCutsPtr_.clear();
}

12
src/dynamicMesh/meshCut/cellCuts/cellCuts.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -128,7 +128,7 @@ class cellCuts
//- Cuts per existing face (includes those along edge of face)
// Cuts in no particular order.
mutable labelListList* faceCutsPtr_;
mutable autoPtr<labelListList> faceCutsPtr_;
//- Per face : cut across edge (so not along existing edge)
// (can only be one per face)
@ -418,6 +418,10 @@ class cellCuts
const List<scalarField>& cellLoopWeights
);
//- Update basic cut information to be consistent across
// coupled boundaries.
void syncProc();
//- Cut cells and update basic cut information from cellLoops.
// Checks each loop for consistency with existing cut pattern.
void setFromCellCutter
@ -557,11 +561,11 @@ public:
// Cuts in no particular order
const labelListList& faceCuts() const
{
if (!faceCutsPtr_)
if (!faceCutsPtr_.valid())
{
calcFaceCuts();
}
return *faceCutsPtr_;
return faceCutsPtr_();
}
//- Gives for split face the two cuts that split the face into two.

10
src/dynamicMesh/meshCut/directions/directions.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -257,7 +257,11 @@ Foam::vectorField Foam::directions::propagateDirection
}
}
Pout<< "Calculated local coords for " << defaultDir
reduce(nGeom, sumOp<label>());
reduce(nTopo, sumOp<label>());
reduce(nUnset, sumOp<label>());
Info<< "Calculated local coords for " << defaultDir
<< endl
<< " Geometric cut cells : " << nGeom << endl
<< " Topological cut cells : " << nTopo << endl
@ -322,7 +326,7 @@ Foam::directions::directions
vector normal = tan1 ^ tan2;
normal /= mag(normal);
Pout<< "Global Coordinate system:" << endl
Info<< "Global Coordinate system:" << endl
<< " normal : " << normal << endl
<< " tan1 : " << tan1 << endl
<< " tan2 : " << tan2

48
src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -33,6 +33,7 @@ License
#include "polyAddPoint.H"
#include "polyAddFace.H"
#include "polyAddCell.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -536,7 +537,7 @@ void Foam::meshCutter::setRefinement
addedPoints_.clear();
addedPoints_.resize(cuts.nLoops());
if (cuts.nLoops() == 0)
if (returnReduce(cuts.nLoops(), sumOp<label>()) == 0)
{
return;
}
@ -544,6 +545,40 @@ void Foam::meshCutter::setRefinement
const labelListList& anchorPts = cuts.cellAnchorPoints();
const labelListList& cellLoops = cuts.cellLoops();
if (debug)
{
// Check that any edge is cut only if any cell using it is cut
boolList edgeOnCutCell(mesh().nEdges(), false);
forAll(cuts.cellLoops(), celli)
{
if (cuts.cellLoops()[celli].size())
{
const labelList& cEdges = mesh().cellEdges(celli);
forAll(cEdges, i)
{
edgeOnCutCell[cEdges[i]] = true;
}
}
}
syncTools::syncEdgeList(mesh(), edgeOnCutCell, orEqOp<bool>(), false);
forAll(cuts.edgeIsCut(), edgeI)
{
if (cuts.edgeIsCut()[edgeI] && !edgeOnCutCell[edgeI])
{
const edge& e = mesh().edges()[edgeI];
WarningInFunction
<< "Problem: cut edge but none of the cells using"
<< " it is cut\n"
<< "edge:" << edgeI << " verts:" << e
<< " at:" << e.line(mesh().points())
<< endl; //abort(FatalError);
}
}
}
//
// Add new points along cut edges.
//
@ -554,15 +589,6 @@ void Foam::meshCutter::setRefinement
{
const edge& e = mesh().edges()[edgeI];
// Check if there is any cell using this edge.
if (debug && findCutCell(cuts, mesh().edgeCells()[edgeI]) == -1)
{
FatalErrorInFunction
<< "Problem: cut edge but none of the cells using it is\n"
<< "edge:" << edgeI << " verts:" << e
<< abort(FatalError);
}
// One of the edge end points should be master point of nbCelli.
label masterPointi = e.start();

4
src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -270,7 +270,7 @@ Foam::Map<Foam::label> Foam::refinementIterator::setRefinement
while (!stop);
if (nRefCells == oldRefCells)
if (returnReduce((nRefCells == oldRefCells), andOp<bool>()))
{
WarningInFunction
<< "stopped refining."