ENH: snappyHexMesh: handle parallel cyclics. Fixes #1731.

This commit is contained in:
mattijs
2020-07-08 16:47:26 +01:00
parent 423bbcefa8
commit c06c63f49f
7 changed files with 532 additions and 111 deletions

View File

@ -17,6 +17,7 @@ meshRefinement/meshRefinementGapRefine.C
meshRefinement/meshRefinementBlock.C
meshRefinement/wallPoints.C
meshRefinement/patchFaceOrientation.C
meshRefinement/weightedPosition.C
refinementFeatures/refinementFeatures.C
refinementSurfaces/surfaceZonesInfo.C

View File

@ -1093,7 +1093,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
newPoints[meshPoints[i]] += disp[i];
}
syncTools::syncPointList
syncTools::syncPointPositions
(
mesh_,
newPoints,

View File

@ -0,0 +1,153 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "weightedPosition.H"
#include "vectorTensorTransform.H"
#include "coupledPolyPatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::weightedPosition Foam::pTraits<Foam::weightedPosition>::zero
(
scalar(0),
Foam::point::zero
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::weightedPosition::weightedPosition()
:
Tuple2<scalar, point>()
{}
Foam::weightedPosition::weightedPosition(const scalar s, const point& p)
:
Tuple2<scalar, point>(s, p)
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::weightedPosition::getPoints
(
const UList<weightedPosition>& in,
List<point>& out
)
{
out.setSize(in.size());
forAll(in, i)
{
out[i] = in[i].second();
if (mag(in[i].first()) > VSMALL)
{
out[i] /= in[i].first();
}
}
}
void Foam::weightedPosition::setPoints
(
const UList<point>& in,
List<weightedPosition>& out
)
{
out.setSize(in.size());
forAll(in, i)
{
out[i].second() = out[i].first()*in[i];
}
}
void Foam::weightedPosition::plusEqOp
(
weightedPosition& x,
const weightedPosition& y
)
{
x.first() += y.first();
x.second() += y.second();
}
void Foam::weightedPosition::operator()
(
const vectorTensorTransform& vt,
const bool forward,
List<weightedPosition>& fld
) const
{
pointField pfld;
getPoints(fld, pfld);
if (forward)
{
pfld = vt.transformPosition(pfld);
}
else
{
pfld = vt.invTransformPosition(pfld);
}
setPoints(pfld, fld);
}
void Foam::weightedPosition::operator()
(
const vectorTensorTransform& vt,
const bool forward,
List<List<weightedPosition>>& flds
) const
{
for (List<weightedPosition>& fld : flds)
{
operator()(vt, forward, fld);
}
}
void Foam::weightedPosition::operator()
(
const coupledPolyPatch& cpp,
Field<weightedPosition>& fld
) const
{
pointField pfld;
getPoints(fld, pfld);
cpp.transformPosition(pfld);
setPoints(pfld, fld);
}
// ************************************************************************* //

View File

@ -0,0 +1,180 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::weightedPosition
Description
Wrapper for position + weight to be used in e.g. averaging.
This avoids the problems when synchronising locations with e.g.
parallel cyclics. The separation vector only applies to locations
and not e.g. summed locations. Note that there is no corresponding
problem for rotational cyclics.
Typical use might be to e.g. average face centres to points on a patch
const labelListList& pointFaces = pp.pointFaces();
const pointField& faceCentres = pp.faceCentres();
Field<weightedPosition> avgBoundary(pointFaces.size());
forAll(pointFaces, pointi)
{
const labelList& pFaces = pointFaces[pointi];
avgBoundary[pointi].first() = pFaces.size();
avgBoundary[pointi].second() = sum(pointField(faceCentres, pFaces));
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
avgBoundary,
weightedPosition::plusEqOp, // combine op
pTraits<weightedPosition>::zero,// null value (not used)
pTraits<weightedPosition>::zero // transform class
);
SourceFiles
weightedPosition.C
\*---------------------------------------------------------------------------*/
#ifndef weightedPosition_H
#define weightedPosition_H
#include "Tuple2.H"
#include "point.H"
#include "Field.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declarations
class weightedPosition;
class vectorTensorTransform;
class coupledPolyPatch;
//- pTraits
template<>
class pTraits<weightedPosition>
{
public:
typedef weightedPosition cmptType;
static const weightedPosition zero;
};
/*---------------------------------------------------------------------------*\
Class weightedPosition Declaration
\*---------------------------------------------------------------------------*/
class weightedPosition
:
public Tuple2<scalar, point>
{
public:
// Constructors
//- Construct null
weightedPosition();
//- Construct from components
weightedPosition(const scalar s, const point& p);
// Member Functions
// Helpers
//- Get points
static void getPoints
(
const UList<weightedPosition>& in,
List<point>& out
);
//- Set points
static void setPoints
(
const UList<point>& in,
List<weightedPosition>& out
);
//- Summation operator
static void plusEqOp(weightedPosition& x, const weightedPosition& y);
// Transformation operators
void operator()
(
const vectorTensorTransform& vt,
const bool forward,
List<weightedPosition>& fld
) const;
void operator()
(
const vectorTensorTransform& vt,
const bool forward,
List<List<weightedPosition>>& flds
) const;
void operator()
(
const coupledPolyPatch& cpp,
Field<weightedPosition>& fld
) const;
template<template<class> class Container>
void operator()
(
const coupledPolyPatch& cpp,
Container<weightedPosition>& map
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "weightedPositionTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,61 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "vectorTensorTransform.H"
#include "coupledPolyPatch.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<template<class> class Container>
void Foam::weightedPosition::operator()
(
const coupledPolyPatch& cpp,
Container<weightedPosition>& map
)
const
{
Field<point> fld(map.size());
label i = 0;
forAllConstIters(map, iter)
{
point pt(iter->second());
if (mag(iter->first()) > VSMALL)
{
pt /= iter->first();
}
fld[i++] = pt;
}
cpp.transformPosition(fld);
i = 0;
forAllIters(map, iter)
{
iter->second() = fld[i++]*iter->first();
}
}
// ************************************************************************* //

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
Copyright (C) 2015-2019 OpenCFD Ltd.
Copyright (C) 2015-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -48,6 +48,7 @@ Description
#include "localPointRegion.H"
#include "PatchTools.H"
#include "refinementFeatures.H"
#include "weightedPosition.H"
#include "profiling.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -234,8 +235,11 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothInternalDisplacement
// Calculate average of connected cells
labelList nCells(mesh.nPoints(), Zero);
pointField sumLocation(mesh.nPoints(), Zero);
Field<weightedPosition> sumLocation
(
mesh.nPoints(),
pTraits<weightedPosition>::zero
);
forAll(isMovingPoint, pointi)
{
@ -243,22 +247,22 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothInternalDisplacement
{
const labelList& pCells = mesh.pointCells(pointi);
forAll(pCells, i)
sumLocation[pointi].first() = pCells.size();
for (const label celli : pCells)
{
sumLocation[pointi] += mesh.cellCentres()[pCells[i]];
nCells[pointi]++;
sumLocation[pointi].second() += mesh.cellCentres()[celli];
}
}
}
// Sum
syncTools::syncPointList(mesh, nCells, plusEqOp<label>(), label(0));
syncTools::syncPointList
(
mesh,
sumLocation,
plusEqOp<point>(),
vector::zero
weightedPosition::plusEqOp, // combine op
pTraits<weightedPosition>::zero,// null value (not used)
pTraits<weightedPosition>::zero // transform class
);
tmp<pointField> tdisplacement(new pointField(mesh.nPoints(), Zero));
@ -268,10 +272,12 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothInternalDisplacement
forAll(displacement, pointi)
{
if (nCells[pointi] > 0)
const weightedPosition& wp = sumLocation[pointi];
if (mag(wp.first()) > VSMALL)
{
displacement[pointi] =
sumLocation[pointi]/nCells[pointi]-mesh.points()[pointi];
wp.second()/wp.first()
- mesh.points()[pointi];
nAdapted++;
}
}
@ -355,56 +361,61 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
// Get average position of boundary face centres
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vectorField avgBoundary(pointFaces.size(), Zero);
labelList nBoundary(pointFaces.size(), Zero);
forAll(pointFaces, patchPointi)
Field<weightedPosition> avgBoundary
(
pointFaces.size(),
pTraits<weightedPosition>::zero
);
{
const labelList& pFaces = pointFaces[patchPointi];
forAll(pFaces, pfi)
forAll(pointFaces, patchPointi)
{
label facei = pFaces[pfi];
const labelList& pFaces = pointFaces[patchPointi];
if (isMasterFace.test(pp.addressing()[facei]))
forAll(pFaces, pfi)
{
avgBoundary[patchPointi] += pp[facei].centre(points);
nBoundary[patchPointi]++;
label facei = pFaces[pfi];
if (isMasterFace.test(pp.addressing()[facei]))
{
avgBoundary[patchPointi].first() += 1.0;
avgBoundary[patchPointi].second() +=
pp[facei].centre(points);
}
}
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
avgBoundary,
plusEqOp<point>(), // combine op
vector::zero // null value
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
nBoundary,
plusEqOp<label>(), // combine op
label(0) // null value
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
avgBoundary,
weightedPosition::plusEqOp, // combine op
pTraits<weightedPosition>::zero,// null value (not used)
pTraits<weightedPosition>::zero // transform class
);
forAll(avgBoundary, i)
{
avgBoundary[i] /= nBoundary[i];
// Normalise
forAll(avgBoundary, i)
{
// Note: what if there is no master boundary face?
if (mag(avgBoundary[i].first()) > VSMALL)
{
avgBoundary[i].second() /= avgBoundary[i].first();
}
}
}
// Get average position of internal face centres
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
vectorField avgInternal;
labelList nInternal;
Field<weightedPosition> avgInternal;
{
vectorField globalSum(mesh.nPoints(), Zero);
labelList globalNum(mesh.nPoints(), Zero);
Field<weightedPosition> globalSum
(
mesh.nPoints(),
pTraits<weightedPosition>::zero
);
// Note: no use of pointFaces
const faceList& faces = mesh.faces();
@ -416,8 +427,9 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
forAll(f, fp)
{
globalSum[f[fp]] += fc;
globalNum[f[fp]]++;
weightedPosition& wp = globalSum[f[fp]];
wp.first() += 1.0;
wp.second() += fc;
}
}
@ -444,8 +456,9 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
forAll(f, fp)
{
globalSum[f[fp]] += fc;
globalNum[f[fp]]++;
weightedPosition& wp = globalSum[f[fp]];
wp.first() += 1.0;
wp.second() += fc;
}
}
}
@ -455,35 +468,28 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
(
mesh,
globalSum,
plusEqOp<vector>(), // combine op
vector::zero // null value
);
syncTools::syncPointList
(
mesh,
globalNum,
plusEqOp<label>(), // combine op
label(0) // null value
weightedPosition::plusEqOp, // combine op
pTraits<weightedPosition>::zero,// null value (not used)
pTraits<weightedPosition>::zero // transform class
);
avgInternal.setSize(meshPoints.size());
nInternal.setSize(meshPoints.size());
forAll(avgInternal, patchPointi)
{
label meshPointi = meshPoints[patchPointi];
const weightedPosition& wp = globalSum[meshPointi];
nInternal[patchPointi] = globalNum[meshPointi];
if (nInternal[patchPointi] == 0)
avgInternal[patchPointi].first() = wp.first();
if (mag(wp.first()) < VSMALL)
{
avgInternal[patchPointi] = globalSum[meshPointi];
// Set to zero?
avgInternal[patchPointi].second() = wp.second();
}
else
{
avgInternal[patchPointi] =
globalSum[meshPointi]
/ nInternal[patchPointi];
avgInternal[patchPointi].second() = wp.second()/wp.first();
}
}
}
@ -514,13 +520,15 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
label meshPointi = meshPoints[i];
const point& currentPos = pp.points()[meshPointi];
// Now we have the two average points: avgBoundary and avgInternal
// and how many boundary/internal faces connect to the point
// (nBoundary, nInternal)
// Now we have the two average points and their counts:
// avgBoundary and avgInternal
// Do some blending between the two.
// Note: the following section has some reasoning behind it but the
// blending factors can be experimented with.
const weightedPosition& internal = avgInternal[i];
const weightedPosition& boundary = avgBoundary[i];
point newPos;
if (!nonManifoldPoint.test(i))
@ -532,14 +540,17 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
point avgPos =
(
internalBlend*nInternal[i]*avgInternal[i]
+(1-internalBlend)*nBoundary[i]*avgBoundary[i]
internalBlend*internal.first()*internal.second()
+(1-internalBlend)*boundary.first()*boundary.second()
)
/ (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
/ (
internalBlend*internal.first()
+(1-internalBlend)*boundary.first()
);
newPos = (1-blend)*avgPos + blend*currentPos;
}
else if (nInternal[i] == 0)
else if (internal.first() == 0)
{
// Non-manifold without internal faces. Use any connected cell
// as internal point instead. Use precalculated any cell to avoid
@ -550,7 +561,7 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
scalar cellCBlend = 0.8;
scalar blend = 0.1;
point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
point avgPos = (1-cellCBlend)*boundary.second() + cellCBlend*cc;
newPos = (1-blend)*avgPos + blend*currentPos;
}
@ -561,8 +572,8 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::smoothPatchDisplacement
scalar blend = 0.1;
point avgPos =
internalBlend*avgInternal[i]
+ (1-internalBlend)*avgBoundary[i];
internalBlend*internal.second()
+ (1-internalBlend)*boundary.second();
newPos = (1-blend)*avgPos + blend*currentPos;
}
@ -980,26 +991,22 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::avgCellCentres
{
const labelListList& pointFaces = pp.pointFaces();
tmp<pointField> tavgBoundary
Field<weightedPosition> avgBoundary
(
new pointField(pointFaces.size(), Zero)
pointFaces.size(),
pTraits<weightedPosition>::zero
);
pointField& avgBoundary = tavgBoundary.ref();
labelList nBoundary(pointFaces.size(), Zero);
forAll(pointFaces, pointi)
{
const labelList& pFaces = pointFaces[pointi];
avgBoundary[pointi].first() = pFaces.size();
forAll(pFaces, pfi)
{
label facei = pFaces[pfi];
label meshFacei = pp.addressing()[facei];
label own = mesh.faceOwner()[meshFacei];
avgBoundary[pointi] += mesh.cellCentres()[own];
nBoundary[pointi]++;
label own = mesh.faceOwner()[pp.addressing()[facei]];
avgBoundary[pointi].second() += mesh.cellCentres()[own];
}
}
@ -1008,22 +1015,14 @@ Foam::tmp<Foam::pointField> Foam::snappySnapDriver::avgCellCentres
mesh,
pp.meshPoints(),
avgBoundary,
plusEqOp<point>(), // combine op
vector::zero // null value
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
nBoundary,
plusEqOp<label>(), // combine op
label(0) // null value
weightedPosition::plusEqOp, // combine op
pTraits<weightedPosition>::zero,// null value (not used)
pTraits<weightedPosition>::zero // transform class
);
forAll(avgBoundary, i)
{
avgBoundary[i] /= nBoundary[i];
}
tmp<pointField> tavgBoundary(new pointField(avgBoundary.size()));
weightedPosition::getPoints(avgBoundary, tavgBoundary.ref());
return tavgBoundary;
}

View File

@ -676,15 +676,42 @@ void Foam::snappySnapDriver::calcNearestFacePointProperties
List<point>(),
mapDistribute::transform()
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceCentres,
listPlusEqOp<point>(),
List<point>(),
mapDistribute::transformPosition()
);
{
// Make into displacement before synchronising to avoid any problems
// with parallel cyclics
pointField localPoints(pp.points(), pp.meshPoints());
forAll(pointFaceCentres, pointi)
{
const point& pt = pp.points()[pp.meshPoints()[pointi]];
List<point>& pFc = pointFaceCentres[pointi];
for (point& p : pFc)
{
p -= pt;
}
}
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
pointFaceCentres,
listPlusEqOp<point>(),
List<point>(),
mapDistribute::transform()
);
forAll(pointFaceCentres, pointi)
{
const point& pt = pp.points()[pp.meshPoints()[pointi]];
List<point>& pFc = pointFaceCentres[pointi];
for (point& p : pFc)
{
p += pt;
}
}
}
syncTools::syncPointList
(
mesh,