Commit a001c0d2 authored by mattijs's avatar mattijs

BUG: volPointInterpolation: fixes #1831

parent ab551573
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -291,9 +291,11 @@ void Foam::volPointInterpolation::interpolateBoundaryField
// Do points on 'normal' patches from the surrounding patch faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(boundary.meshPoints(), i)
const labelList& mp = boundary.meshPoints();
forAll(mp, i)
{
label pointi = boundary.meshPoints()[i];
label pointi = mp[i];
if (isPatchPoint_[pointi])
{
......@@ -319,6 +321,17 @@ void Foam::volPointInterpolation::interpolateBoundaryField
// And add separated contributions
addSeparated(pf);
// Optionally normalise
if (normalisationPtr_)
{
const scalarField& normalisation = normalisationPtr_();
forAll(mp, i)
{
pfi[mp[i]] *= normalisation[i];
}
}
// Push master data to slaves. It is possible (not sure how often) for
// a coupled point to have its master on a different patch so
// to make sure just push master data to slaves.
......
......@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -32,6 +33,7 @@ License
#include "demandDrivenData.H"
#include "pointConstraints.H"
#include "surfaceFields.H"
#include "processorPointPatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -43,6 +45,25 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::volPointInterpolation::hasSeparated(const pointMesh& pMesh)
{
const pointBoundaryMesh& pbm = pMesh.boundary();
bool hasSpecial = false;
for (const auto& pp : pbm)
{
if (isA<coupledFacePointPatch>(pp) && !isType<processorPointPatch>(pp))
{
hasSpecial = true;
break;
}
}
reduce(hasSpecial, orOp<bool>());
return hasSpecial;
}
void Foam::volPointInterpolation::calcBoundaryAddressing()
{
if (debug)
......@@ -70,8 +91,10 @@ void Foam::volPointInterpolation::calcBoundaryAddressing()
boundaryIsPatchFace_.setSize(boundary.size());
boundaryIsPatchFace_ = false;
isPatchPoint_.setSize(mesh().nPoints());
isPatchPoint_ = false;
// Store per mesh point whether it is on any 'real' patch. Currently
// boolList just so we can use syncUntransformedData (does not take
// bitSet. Tbd)
boolList isPatchPoint(mesh().nPoints(), false);
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
......@@ -99,7 +122,7 @@ void Foam::volPointInterpolation::calcBoundaryAddressing()
forAll(f, fp)
{
isPatchPoint_[f[fp]] = true;
isPatchPoint[f[fp]] = true;
}
}
}
......@@ -107,28 +130,19 @@ void Foam::volPointInterpolation::calcBoundaryAddressing()
// Make sure point status is synchronised so even processor that holds
// no face of a certain patch still can have boundary points marked.
if (debug)
{
boolList oldData(isPatchPoint_);
pointConstraints::syncUntransformedData
(
mesh(),
isPatchPoint_,
orEqOp<bool>()
);
pointConstraints::syncUntransformedData
(
mesh(),
isPatchPoint,
orEqOp<bool>()
);
forAll(isPatchPoint_, pointi)
{
if (isPatchPoint_[pointi] != oldData[pointi])
{
Pout<< "volPointInterpolation::calcBoundaryAddressing():"
<< " added dangling mesh point:" << pointi
<< " at:" << mesh().points()[pointi]
<< endl;
}
}
// Convert to bitSet
isPatchPoint_.setSize(mesh().nPoints());
isPatchPoint_.assign(isPatchPoint);
if (debug)
{
label nPatchFace = 0;
forAll(boundaryIsPatchFace_, i)
{
......@@ -242,6 +256,98 @@ void Foam::volPointInterpolation::makeBoundaryWeights(scalarField& sumWeights)
}
void Foam::volPointInterpolation::interpolateOne
(
const tmp<scalarField>& tnormalisation,
pointScalarField& pf
) const
{
if (debug)
{
Pout<< "volPointInterpolation::interpolateOne("
<< "pointScalarField&) : "
<< "interpolating oneField"
<< " from cells to BOUNDARY points "
<< pf.name() << endl;
}
const primitivePatch& boundary = boundaryPtr_();
const labelList& mp = boundary.meshPoints();
Field<scalar>& pfi = pf.primitiveFieldRef();
// 1. Interpolate coupled boundary points from cells
{
forAll(mp, i)
{
const label pointi = mp[i];
if (!isPatchPoint_[pointi])
{
const scalarList& pw = pointWeights_[pointi];
scalar& val = pfi[pointi];
val = Zero;
forAll(pw, pointCelli)
{
val += pw[pointCelli];
}
}
}
}
// 2. Interpolate to the patches preserving fixed value BCs
{
forAll(mp, i)
{
const label pointi = mp[i];
if (isPatchPoint_[pointi])
{
const labelList& pFaces = boundary.pointFaces()[i];
const scalarList& pWeights = boundaryPointWeights_[i];
scalar& val = pfi[pointi];
val = Zero;
forAll(pFaces, j)
{
if (boundaryIsPatchFace_[pFaces[j]])
{
val += pWeights[j];
}
}
}
}
// Sum collocated contributions
pointConstraints::syncUntransformedData
(
mesh(),
pfi,
plusEqOp<scalar>()
);
// And add separated contributions
addSeparated(pf);
// Optionally normalise
if (tnormalisation)
{
const scalarField& normalisation = tnormalisation();
forAll(mp, i)
{
pfi[mp[i]] *= normalisation[i];
}
}
}
// Apply constraints
pointConstraints::New(pf.mesh()).constrain(pf, false);
}
void Foam::volPointInterpolation::makeWeights()
{
if (debug)
......@@ -251,22 +357,28 @@ void Foam::volPointInterpolation::makeWeights()
<< endl;
}
const pointMesh& pMesh = pointMesh::New(mesh());
// Update addressing over all boundary faces
calcBoundaryAddressing();
// Running sum of weights
pointScalarField sumWeights
tmp<pointScalarField> tsumWeights
(
IOobject
new pointScalarField
(
"volPointSumWeights",
mesh().polyMesh::instance(),
mesh()
),
pointMesh::New(mesh()),
dimensionedScalar(dimless, Zero)
IOobject
(
"volPointSumWeights",
mesh().polyMesh::instance(),
mesh()
),
pMesh,
dimensionedScalar(dimless, Zero)
)
);
pointScalarField& sumWeights = tsumWeights.ref();
// Create internal weights; add to sumWeights
......@@ -277,19 +389,8 @@ void Foam::volPointInterpolation::makeWeights()
makeBoundaryWeights(sumWeights);
//forAll(boundary.meshPoints(), i)
//{
// label pointi = boundary.meshPoints()[i];
//
// if (isPatchPoint_[pointi])
// {
// Pout<< "Calculated Weight at boundary point:" << i
// << " at:" << mesh().points()[pointi]
// << " sumWeight:" << sumWeights[pointi]
// << " from:" << boundaryPointWeights_[i]
// << endl;
// }
//}
const primitivePatch& boundary = boundaryPtr_();
const labelList& mp = boundary.meshPoints();
// Sum collocated contributions
......@@ -300,9 +401,11 @@ void Foam::volPointInterpolation::makeWeights()
plusEqOp<scalar>()
);
// And add separated contributions
addSeparated(sumWeights);
// Push master data to slaves. It is possible (not sure how often) for
// a coupled point to have its master on a different patch so
// to make sure just push master data to slaves. Reuse the syncPointData
......@@ -322,11 +425,9 @@ void Foam::volPointInterpolation::makeWeights()
}
// Normalise boundary weights
const primitivePatch& boundary = boundaryPtr_();
forAll(boundary.meshPoints(), i)
forAll(mp, i)
{
label pointi = boundary.meshPoints()[i];
const label pointi = mp[i];
scalarList& pw = boundaryPointWeights_[i];
// Note:pw only sized for isPatchPoint
......@@ -337,6 +438,30 @@ void Foam::volPointInterpolation::makeWeights()
}
// Normalise separated contributions
if (hasSeparated_)
{
if (debug)
{
Pout<< "volPointInterpolation::makeWeights() : "
<< "detected separated coupled patches"
<< " - allocating normalisation" << endl;
}
// Sum up effect of interpolating one (on boundary points only)
interpolateOne(tmp<scalarField>(), sumWeights);
// Store as normalisation factor (on boundary points only)
normalisationPtr_ = new scalarField(mp.size());
normalisationPtr_.ref() = scalar(1.0);
normalisationPtr_.ref() /= scalarField(sumWeights, mp);
}
else
{
normalisationPtr_.clear();
}
if (debug)
{
Pout<< "volPointInterpolation::makeWeights() : "
......@@ -350,7 +475,8 @@ void Foam::volPointInterpolation::makeWeights()
Foam::volPointInterpolation::volPointInterpolation(const fvMesh& vm)
:
MeshObject<fvMesh, Foam::UpdateableMeshObject, volPointInterpolation>(vm)
MeshObject<fvMesh, Foam::UpdateableMeshObject, volPointInterpolation>(vm),
hasSeparated_(hasSeparated(pointMesh::New(vm)))
{
makeWeights();
}
......@@ -366,6 +492,9 @@ Foam::volPointInterpolation::~volPointInterpolation()
void Foam::volPointInterpolation::updateMesh(const mapPolyMesh&)
{
// Recheck whether has coupled patches
hasSeparated_ = hasSeparated(pointMesh::New(mesh()));
makeWeights();
}
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016 OpenCFD Ltd.
Copyright (C) 2016-2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -72,19 +72,30 @@ class volPointInterpolation
//- Boundary addressing
autoPtr<primitivePatch> boundaryPtr_;
//- Per boundary face whether is on non-coupled patch
boolList boundaryIsPatchFace_;
//- Per boundary face whether is on non-coupled, non-empty patch
bitSet boundaryIsPatchFace_;
//- Per mesh(!) point whether is on non-coupled patch (on any
// processor)
boolList isPatchPoint_;
//- Per mesh(!) point whether is on non-coupled, non-empty patch (on
// any processor)
bitSet isPatchPoint_;
//- Per boundary point the weights per pointFaces.
scalarListList boundaryPointWeights_;
//- Whether mesh has any non-processor coupled patches
bool hasSeparated_;
//- Optional scaling for coupled point patch fields using the
// swapAddSeparated to add additional contributions (e.g. cyclicAMI).
// (note: one-to-one matching (processor) uses a different mechanism)
tmp<scalarField> normalisationPtr_;
// Private Member Functions
//- Constructor helper: check if any non-processor coupled patches
static bool hasSeparated(const pointMesh& pMesh);
//- Construct addressing over all boundary faces
void calcBoundaryAddressing();
......@@ -97,6 +108,14 @@ class volPointInterpolation
//- Construct all point weighting factors
void makeWeights();
//- Helper: interpolate oneField onto boundary points only. Used to
// correct normalisation
void interpolateOne
(
const tmp<scalarField>& tnormalisation,
pointScalarField& pf
) const;
//- Helper: push master point data to collocated points
template<class Type>
void pushUntransformedData(List<Type>&) const;
......@@ -209,70 +228,71 @@ public:
// Low level
//- Interpolate internal field from volField to pointField
// using inverse distance weighting
template<class Type>
void interpolateInternalField
(
const GeometricField<Type, fvPatchField, volMesh>&,
GeometricField<Type, pointPatchField, pointMesh>&
) const;
//- Interpolate boundary field without applying constraints/boundary
// conditions
template<class Type>
void interpolateBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
GeometricField<Type, pointPatchField, pointMesh>& pf
) const;
//- Interpolate boundary with constraints/boundary conditions
template<class Type>
void interpolateBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
GeometricField<Type, pointPatchField, pointMesh>& pf,
const bool overrideFixedValue
) const;
//- Interpolate from volField to pointField
// using inverse distance weighting
template<class Type>
void interpolate
(
const GeometricField<Type, fvPatchField, volMesh>&,
GeometricField<Type, pointPatchField, pointMesh>&
) const;
//- Interpolate volField using inverse distance weighting
// returning pointField with name. Optionally caches
template<class Type>
tmp<GeometricField<Type, pointPatchField, pointMesh>> interpolate
(
const GeometricField<Type, fvPatchField, volMesh>&,
const word& name,
const bool cache
) const;
//- Interpolate dimensioned internal field from cells to points
// using inverse distance weighting
template<class Type>
void interpolateDimensionedInternalField
(
const DimensionedField<Type, volMesh>& vf,
DimensionedField<Type, pointMesh>& pf
) const;
//- Interpolate dimensionedField using inverse distance weighting
// returning pointField with name. Optionally caches
template<class Type>
tmp<DimensionedField<Type, pointMesh>> interpolate
(
const DimensionedField<Type, volMesh>&,
const word& name,
const bool cache
) const;
//- Interpolate internal field from volField to pointField
// using inverse distance weighting
template<class Type>
void interpolateInternalField
(
const GeometricField<Type, fvPatchField, volMesh>&,
GeometricField<Type, pointPatchField, pointMesh>&
) const;
//- Interpolate boundary field without applying constraints/boundary
// conditions
template<class Type>
void interpolateBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
GeometricField<Type, pointPatchField, pointMesh>& pf
) const;
//- Interpolate boundary with constraints/boundary conditions
template<class Type>
void interpolateBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
GeometricField<Type, pointPatchField, pointMesh>& pf,
const bool overrideFixedValue
) const;
//- Interpolate from volField to pointField
// using inverse distance weighting
template<class Type>
void interpolate
(
const GeometricField<Type, fvPatchField, volMesh>&,
GeometricField<Type, pointPatchField, pointMesh>&
) const;
//- Interpolate volField using inverse distance weighting
// returning pointField with name. Optionally caches
template<class Type>
tmp<GeometricField<Type, pointPatchField, pointMesh>> interpolate
(
const GeometricField<Type, fvPatchField, volMesh>&,
const word& name,
const bool cache
) const;
//- Interpolate dimensioned internal field from cells to points
// using inverse distance weighting
template<class Type>
void interpolateDimensionedInternalField
(
const DimensionedField<Type, volMesh>& vf,
DimensionedField<Type, pointMesh>& pf
) const;
//- Interpolate dimensionedField using inverse distance weighting
// returning pointField with name. Optionally caches
template<class Type>
tmp<DimensionedField<Type, pointMesh>> interpolate
(
const DimensionedField<Type, volMesh>&,
const word& name,
const bool cache
) const;
// Interpolation for displacement (applies 2D correction)
......@@ -283,7 +303,6 @@ public:
const volVectorField&,
pointVectorField&
) const;
};
......
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