mappedPatchBase: Fixed mapping with interpolation

Mapping with interpolation now behaves correctly when a single cell maps
to multiple faces. In addition, the mapping structure is more compact
and no longer results in copies being made of entire internal fields.

This has been achieved by rewriting the mapping strategy in
mappedPatchBase so that it maps from a reduced set of sampling locations
to the patch faces, rather than directly from the cells/faces to the
patch faces. This is more efficient, but it also permits multiple
sampling locations to be sent to a single cell/face. Values can then be
interpolated to these points before being sent back to the patch faces.

Previously a single cell/face could only be sent a single location onto
which to interpolate; typically that of the first associated patch face.
The resulting interpolated value was then sent back to all associated
patch faces. This meant that some (potentially most) patch faces did
receive a value interpolated to the correct position.
This commit is contained in:
Will Bainbridge
2022-08-18 16:23:16 +01:00
parent 7e798e86d0
commit c7ccc2dee9
28 changed files with 1011 additions and 980 deletions

View File

@ -846,16 +846,32 @@ void Foam::distributionMapBase::compact
// from the submap and do the same to the constructMap locally
// (and in same order).
// Send elemIsUsed field to neighbour. Use nonblocking code from
// distributionMapBase but in reverse order.
// Send elemIsUsed field to neighbours
List<boolList> recvFields(Pstream::nProcs());
// "Send" to myself (i.e., write directly into recvFields)
{
const labelList& map = constructMap_[Pstream::myProcNo()];
recvFields[Pstream::myProcNo()].setSize(map.size());
forAll(map, i)
{
recvFields[Pstream::myProcNo()][i] = accessAndFlip
(
elemIsUsed,
map[i],
constructHasFlip_,
noOp() // do not flip elemIsUsed value
);
}
}
// Send to others. Use nonblocking code from distributionMapBase but in
// reverse order.
if (Pstream::parRun())
{
label startOfRequests = Pstream::nRequests();
// Set up receives from neighbours
List<boolList> recvFields(Pstream::nProcs());
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& map = subMap_[domain];
@ -874,7 +890,6 @@ void Foam::distributionMapBase::compact
}
}
List<boolList> sendFields(Pstream::nProcs());
for (label domain = 0; domain < Pstream::nProcs(); domain++)
@ -907,62 +922,35 @@ void Foam::distributionMapBase::compact
}
}
Pstream::waitRequests(startOfRequests);
}
// Compact out all submap entries that are referring to unused elements
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& map = subMap_[domain];
// Set up 'send' to myself - write directly into recvFields
labelList newMap(map.size());
label newI = 0;
forAll(map, i)
{
const labelList& map = constructMap_[Pstream::myProcNo()];
recvFields[Pstream::myProcNo()].setSize(map.size());
forAll(map, i)
if (recvFields[domain][i])
{
recvFields[Pstream::myProcNo()][i] = accessAndFlip
(
elemIsUsed,
map[i],
constructHasFlip_,
noOp() // do not flip elemIsUsed value
);
// So element is used on destination side
newMap[newI++] = map[i];
}
}
// Wait for all to finish
Pstream::waitRequests(startOfRequests);
// Compact out all submap entries that are referring to unused elements
for (label domain = 0; domain < Pstream::nProcs(); domain++)
if (newI < map.size())
{
const labelList& map = subMap_[domain];
labelList newMap(map.size());
label newI = 0;
forAll(map, i)
{
if (recvFields[domain][i])
{
// So element is used on destination side
newMap[newI++] = map[i];
}
}
if (newI < map.size())
{
newMap.setSize(newI);
subMap_[domain].transfer(newMap);
}
newMap.setSize(newI);
subMap_[domain].transfer(newMap);
}
}
// 2. remove from construct map - since end-result (element in elemIsUsed)
// not used.
label maxConstructIndex = -1;
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& map = constructMap_[domain];
@ -1013,16 +1001,32 @@ void Foam::distributionMapBase::compact
// from the submap and do the same to the constructMap locally
// (and in same order).
// Send elemIsUsed field to neighbour. Use nonblocking code from
// distributionMapBase but in reverse order.
// Send elemIsUsed field to neighbours
List<boolList> recvFields(Pstream::nProcs());
// "Send" to myself (i.e., write directly into recvFields)
{
const labelList& map = constructMap_[Pstream::myProcNo()];
recvFields[Pstream::myProcNo()].setSize(map.size());
forAll(map, i)
{
recvFields[Pstream::myProcNo()][i] = accessAndFlip
(
elemIsUsed,
map[i],
constructHasFlip_,
noOp() // do not flip elemIsUsed value
);
}
}
// Send to others. Use nonblocking code from distributionMapBase but in
// reverse order.
if (Pstream::parRun())
{
label startOfRequests = Pstream::nRequests();
// Set up receives from neighbours
List<boolList> recvFields(Pstream::nProcs());
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& map = subMap_[domain];
@ -1041,7 +1045,6 @@ void Foam::distributionMapBase::compact
}
}
List<boolList> sendFields(Pstream::nProcs());
for (label domain = 0; domain < Pstream::nProcs(); domain++)
@ -1054,12 +1057,13 @@ void Foam::distributionMapBase::compact
subField.setSize(map.size());
forAll(map, i)
{
label index = map[i];
if (constructHasFlip_)
{
index = mag(index)-1;
}
subField[i] = elemIsUsed[index];
subField[i] = accessAndFlip
(
elemIsUsed,
map[i],
constructHasFlip_,
noOp() // do not flip elemIsUsed value
);
}
OPstream::write
@ -1073,108 +1077,79 @@ void Foam::distributionMapBase::compact
}
}
// Set up 'send' to myself - write directly into recvFields
{
const labelList& map = constructMap_[Pstream::myProcNo()];
recvFields[Pstream::myProcNo()].setSize(map.size());
forAll(map, i)
{
label index = map[i];
if (constructHasFlip_)
{
index = mag(index)-1;
}
recvFields[Pstream::myProcNo()][i] = elemIsUsed[index];
}
}
// Wait for all to finish
Pstream::waitRequests(startOfRequests);
}
// Work out which elements on the sending side are needed
{
oldToNewSub.setSize(localSize, -1);
boolList sendElemIsUsed(localSize, false);
// Work out which elements on the sending side are needed
{
oldToNewSub.setSize(localSize, -1);
boolList sendElemIsUsed(localSize, false);
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& map = subMap_[domain];
forAll(map, i)
{
if (recvFields[domain][i])
{
label index = map[i];
if (subHasFlip_)
{
index = mag(index)-1;
}
sendElemIsUsed[index] = true;
}
}
}
label newI = 0;
forAll(sendElemIsUsed, i)
{
if (sendElemIsUsed[i])
{
oldToNewSub[i] = newI++;
}
}
}
// Compact out all submap entries that are referring to unused elements
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& map = subMap_[domain];
labelList newMap(map.size());
label newI = 0;
forAll(map, i)
{
if (recvFields[domain][i])
{
// So element is used on destination side
label index = map[i];
label sign = 1;
if (subHasFlip_)
{
if (index < 0)
{
sign = -1;
}
index = mag(index)-1;
}
label newIndex = oldToNewSub[index];
if (subHasFlip_)
{
newIndex = sign*(newIndex+1);
}
newMap[newI++] = newIndex;
sendElemIsUsed[index] = true;
}
}
newMap.setSize(newI);
subMap_[domain].transfer(newMap);
}
label newI = 0;
forAll(sendElemIsUsed, i)
{
if (sendElemIsUsed[i])
{
oldToNewSub[i] = newI++;
}
}
}
// Compact out all submap entries that are referring to unused elements
for (label domain = 0; domain < Pstream::nProcs(); domain++)
{
const labelList& map = subMap_[domain];
labelList newMap(map.size());
label newI = 0;
forAll(map, i)
{
if (recvFields[domain][i])
{
// So element is used on destination side
label index = map[i];
label sign = 1;
if (subHasFlip_)
{
if (index < 0)
{
sign = -1;
}
index = mag(index)-1;
}
label newIndex = oldToNewSub[index];
if (subHasFlip_)
{
newIndex = sign*(newIndex+1);
}
newMap[newI++] = newIndex;
}
}
newMap.setSize(newI);
subMap_[domain].transfer(newMap);
}
// 2. remove from construct map - since end-result (element in elemIsUsed)
// not used.
oldToNewConstruct.setSize(elemIsUsed.size(), -1);
constructSize_ = 0;
forAll(elemIsUsed, i)

View File

@ -0,0 +1,169 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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::RemoteData
Description
Struct for keeping processor, element (cell, face, point) and a piece of
data. Used for finding minimum values across multiple processes.
SourceFiles
RemoteDataI.H
\*---------------------------------------------------------------------------*/
#ifndef RemoteData_H
#define RemoteData_H
#include "Istream.H"
#include "Ostream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of friend functions and operators
template<class Type>
class RemoteData;
template<class Type>
inline bool operator==(const RemoteData<Type>&, const RemoteData<Type>&);
template<class Type>
inline bool operator!=(const RemoteData<Type>&, const RemoteData<Type>&);
template<class Type>
inline Istream& operator>>(Istream&, RemoteData<Type>&);
template<class Type>
inline Ostream& operator<<(Ostream&, const RemoteData<Type>&);
/*---------------------------------------------------------------------------*\
Class RemoteData Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class RemoteData
{
public:
// Public Data
//- Processor index
label proci;
//- Element index
label elementi;
//- Data
Type data;
// Constructors
//- Construct null
inline RemoteData();
//- Construct from components
inline RemoteData(const label, const label, const Type&);
//- Construct from stream
inline RemoteData(Istream& is);
// Public Classes
//- Operator to take smallest valid value
struct smallestEqOp
{
inline void operator()
(
RemoteData<Type>& x,
const RemoteData<Type>& y
) const;
};
//- Operator to take smallest first valid value
struct smallestFirstEqOp
{
inline void operator()
(
RemoteData<Type>& x,
const RemoteData<Type>& y
) const;
};
// Friend Operators
//- Equality comparison
friend bool operator== <Type>
(
const RemoteData<Type>& a,
const RemoteData<Type>& b
);
//- Inequality comparison
friend bool operator!= <Type>
(
const RemoteData<Type>& a,
const RemoteData<Type>& b
);
// IOstream Operators
//- Write to stream
friend Ostream& operator<< <Type>
(
Ostream& os,
const RemoteData<Type>& p
);
//- Read from stream
friend Istream& operator>> <Type>
(
Istream& is,
RemoteData<Type>& p
);
};
template<class Type>
struct pTraits<RemoteData<Type>>
{
typedef RemoteData<Type> cmptType;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "RemoteDataI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,150 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "RemoteData.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
inline Foam::RemoteData<Type>::RemoteData()
:
proci(-1),
elementi(-1),
data()
{}
template<class Type>
inline Foam::RemoteData<Type>::RemoteData
(
const label p,
const label e,
const Type& d
)
:
proci(p),
elementi(e),
data(d)
{}
template<class Type>
inline Foam::RemoteData<Type>::RemoteData(Istream& is)
{
is >> *this;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
template<class Type>
inline void Foam::RemoteData<Type>::smallestEqOp::operator()
(
RemoteData<Type>& x,
const RemoteData<Type>& y
) const
{
if (y.proci != -1)
{
if (x.proci == -1)
{
x = y;
}
else if (y.data < x.data)
{
x = y;
}
}
}
template<class Type>
inline void Foam::RemoteData<Type>::smallestFirstEqOp::operator()
(
RemoteData<Type>& x,
const RemoteData<Type>& y
) const
{
if (y.proci != -1)
{
if (x.proci == -1)
{
x = y;
}
else if (y.data.first() < x.data.first())
{
x = y;
}
}
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template<class Type>
inline bool Foam::operator==
(
const RemoteData<Type>& a,
const RemoteData<Type>& b
)
{
return
a.proci == b.proci
&& a.elementi == b.elementi
&& a.data == b.data;
}
template<class Type>
inline bool Foam::operator!=
(
const RemoteData<Type>& a,
const RemoteData<Type>& b
)
{
return !(a == b);
}
// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * //
template<class Type>
inline Foam::Ostream& Foam::operator<<(Ostream& os, const RemoteData<Type>& p)
{
return os
<< p.proci << token::SPACE
<< p.elementi << token::SPACE
<< p.data;
}
template<class Type>
inline Foam::Istream& Foam::operator>>(Istream& is, RemoteData<Type>& p)
{
return is >> p.proci >> p.elementi >> p.data;
}
// ************************************************************************* //

View File

@ -252,13 +252,7 @@ baffleThickness() const
{
const mappedPatchBase& mpp =
refCast<const mappedPatchBase>(patch().patch());
tmp<scalarField> nbrBaffleThickness
(
new scalarField(nbrField().baffleThickness())
);
mpp.distribute(nbrBaffleThickness.ref());
return nbrBaffleThickness;
return mpp.distribute(nbrField().baffleThickness());
}
}
@ -274,13 +268,7 @@ tmp<scalarField> thermalBaffle1DFvPatchScalarField<solidType>::qs() const
{
const mappedPatchBase& mpp =
refCast<const mappedPatchBase>(patch().patch());
tmp<scalarField> nbrQs
(
new scalarField(nbrField().qs())
);
mpp.distribute(nbrQs.ref());
return nbrQs;
return mpp.distribute(nbrField().qs());
}
}
@ -391,8 +379,7 @@ void thermalBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
scalarField kappaDelta(kappap*patch().deltaCoeffs());
// Neighbour properties
scalarField nbrTp(nbrField());
mpp.distribute(nbrTp);
const scalarField nbrTp(mpp.distribute(nbrField()));
// Solid properties
scalarField kappas(patch().size(), 0.0);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -216,27 +216,22 @@ void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
<< "currently of type " << nbrTp.type() << exit(FatalError);
}
const scalarField myKDelta(kappa(*this)*patch().deltaCoeffs());
const thisType& nbrField = refCast<const thisType>(nbrTp);
// Swap to obtain full local values of neighbour internal field
tmp<scalarField> nbrIntFld(new scalarField(nbrField.size(), 0.0));
tmp<scalarField> nbrKDelta(new scalarField(nbrField.size(), 0.0));
if (contactRes_ == 0.0)
{
nbrIntFld.ref() = nbrField.patchInternalField();
nbrKDelta.ref() = nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs();
}
else
{
nbrIntFld.ref() = nbrField;
nbrKDelta.ref() = contactRes_;
}
mpp.distribute(nbrIntFld.ref());
mpp.distribute(nbrKDelta.ref());
tmp<scalarField> myKDelta = kappa(*this)*patch().deltaCoeffs();
const scalarField nbrIntFld
(
contactRes_ == 0
? mpp.distribute(nbrField.patchInternalField())
: mpp.distribute(nbrField)
);
const scalarField nbrKDelta
(
contactRes_ == 0
? mpp.distribute(nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs())
: tmp<scalarField>(new scalarField(contactRes_, nbrField.size()))
);
// Both sides agree on
// - temperature : (myKDelta*fld + nbrKDelta*nbrFld)/(myKDelta+nbrKDelta)
@ -251,11 +246,11 @@ void turbulentTemperatureCoupledBaffleMixedFvPatchScalarField::updateCoeffs()
// same on both sides. This leads to the choice of
// - refGradient = qs_/kappa;
// - refValue = neighbour value
// - mixFraction = nbrKDelta / (nbrKDelta + myKDelta())
// - mixFraction = nbrKDelta / (nbrKDelta + myKDelta)
this->refValue() = nbrIntFld();
this->refValue() = nbrIntFld;
this->refGrad() = qs_/kappa(*this);
this->valueFraction() = nbrKDelta()/(nbrKDelta() + myKDelta());
this->valueFraction() = nbrKDelta/(nbrKDelta + myKDelta);
mixedFvPatchScalarField::updateCoeffs();

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -182,7 +182,6 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
const fvPatch& nbrPatch =
refCast<const fvMesh>(nbrMesh).boundary()[samplePatchi];
scalarField Tc(patchInternalField());
scalarField& Tp = *this;
@ -204,36 +203,37 @@ void turbulentTemperatureRadCoupledMixedFvPatchScalarField::updateCoeffs()
const thisType& nbrField = refCast<const thisType>(nbrTp);
// Swap to obtain full local values of neighbour internal field
scalarField TcNbr(nbrField.patchInternalField());
mpp.distribute(TcNbr);
const scalarField TcNbr(mpp.distribute(nbrField.patchInternalField()));
const scalarField KDelta(kappa(*this)*patch().deltaCoeffs());
// Swap to obtain full local values of neighbour K*delta
scalarField KDeltaNbr;
if (contactRes_ == 0.0)
{
KDeltaNbr = nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs();
}
else
{
KDeltaNbr.setSize(nbrField.size(), contactRes_);
}
mpp.distribute(KDeltaNbr);
const scalarField KDeltaNbr
(
contactRes_ == 0
? mpp.distribute(nbrField.kappa(nbrField)*nbrPatch.deltaCoeffs())
: tmp<scalarField>(new scalarField(nbrField.size(), contactRes_))
);
scalarField KDelta(kappa(*this)*patch().deltaCoeffs());
const scalarField qr
(
qrName_ != "none"
? static_cast<const scalarField&>
(
patch().lookupPatchField<volScalarField, scalar>(qrName_)
)
: scalarField(Tp.size(), 0)
);
scalarField qr(Tp.size(), 0.0);
if (qrName_ != "none")
{
qr = patch().lookupPatchField<volScalarField, scalar>(qrName_);
}
scalarField qrNbr(Tp.size(), 0.0);
if (qrNbrName_ != "none")
{
qrNbr = nbrPatch.lookupPatchField<volScalarField, scalar>(qrNbrName_);
mpp.distribute(qrNbr);
}
const scalarField qrNbr
(
qrNbrName_ != "none"
? mpp.distribute
(
nbrPatch.lookupPatchField<volScalarField, scalar>(qrNbrName_)
)
: tmp<scalarField>(new scalarField(Tp.size(), 0))
);
valueFraction() = KDeltaNbr/(KDeltaNbr + KDelta);
refValue() = TcNbr;

View File

@ -146,11 +146,7 @@ mappedPatchFieldBase<Type>::sampleField() const
if (fieldName_ == patchField_.internalField().name())
{
// Optimisation: bypass field lookup
return
dynamic_cast<const fieldType&>
(
patchField_.internalField()
);
return dynamic_cast<const fieldType&>(patchField_.internalField());
}
else
{
@ -170,14 +166,13 @@ tmp<Field<Type>> mappedPatchFieldBase<Type>::mappedField() const
{
typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
const fvMesh& nbrMesh = refCast<const fvMesh>(mapper_.sampleMesh());
// Since we're inside initEvaluate/evaluate there might be processor
// comms underway. Change the tag we use.
int oldTag = UPstream::msgType();
UPstream::msgType() = oldTag + 1;
const fvMesh& thisMesh = patchField_.patch().boundaryMesh().mesh();
const fvMesh& nbrMesh = refCast<const fvMesh>(mapper_.sampleMesh());
// Result of obtaining remote values
tmp<Field<Type>> tnewValues(new Field<Type>(0));
Field<Type>& newValues = tnewValues.ref();
@ -186,24 +181,10 @@ tmp<Field<Type>> mappedPatchFieldBase<Type>::mappedField() const
{
case mappedPatchBase::NEARESTCELL:
{
const distributionMap& distMap = mapper_.map();
if (interpolationScheme_ != interpolationCell<Type>::typeName)
{
// Send back sample points to the processor that holds the cell
vectorField samples(mapper_.samplePoints());
distMap.reverseDistribute
(
(
mapper_.sameRegion()
? thisMesh.nCells()
: nbrMesh.nCells()
),
point::max,
samples
);
autoPtr<interpolation<Type>> interpolator
// Create an interpolation
autoPtr<interpolation<Type>> interpolatorPtr
(
interpolation<Type>::New
(
@ -211,28 +192,43 @@ tmp<Field<Type>> mappedPatchFieldBase<Type>::mappedField() const
sampleField()
)
);
const interpolation<Type>& interp = interpolator();
const interpolation<Type>& interpolator = interpolatorPtr();
newValues.setSize(samples.size(), pTraits<Type>::max);
forAll(samples, celli)
// Cells on which samples are generated
const labelList& sampleCells = mapper_.mapIndices();
// Send the patch points to the cells
pointField samplePoints(mapper_.samplePoints());
mapper_.map().reverseDistribute
(
sampleCells.size(),
samplePoints
);
// Interpolate values
Field<Type> sampleValues(sampleCells.size());
forAll(sampleCells, i)
{
if (samples[celli] != point::max)
if (sampleCells[i] != -1)
{
newValues[celli] = interp.interpolate
(
samples[celli],
celli
);
sampleValues[i] =
interpolator.interpolate
(
samplePoints[i],
sampleCells[i]
);
}
}
// Send the values back to the patch
newValues = sampleValues;
mapper_.map().distribute(newValues);
}
else
{
newValues = sampleField();
newValues = mapper_.distribute(sampleField());
}
distMap.distribute(newValues);
break;
}
case mappedPatchBase::NEARESTPATCHFACE:
@ -252,8 +248,8 @@ tmp<Field<Type>> mappedPatchFieldBase<Type>::mappedField() const
const fieldType& nbrField = sampleField();
newValues = nbrField.boundaryField()[nbrPatchID];
mapper_.distribute(newValues);
newValues =
mapper_.distribute(nbrField.boundaryField()[nbrPatchID]);
break;
}
@ -275,8 +271,7 @@ tmp<Field<Type>> mappedPatchFieldBase<Type>::mappedField() const
}
}
mapper_.distribute(allValues);
newValues.transfer(allValues);
newValues = mapper_.distribute(allValues);
break;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -127,8 +127,8 @@ void Foam::mappedFixedInternalValueFvPatchField<Type>::updateCoeffs()
const label samplePatchi = mpp.samplePolyPatch().index();
const fvPatchField<Type>& nbrPatchField =
this->sampleField().boundaryField()[samplePatchi];
nbrIntFld = nbrPatchField.patchInternalField();
mpp.distribute(nbrIntFld);
nbrIntFld = mpp.distribute(nbrPatchField.patchInternalField());
break;
}
@ -151,8 +151,7 @@ void Foam::mappedFixedInternalValueFvPatchField<Type>::updateCoeffs()
}
}
mpp.distribute(allValues);
nbrIntFld.transfer(allValues);
nbrIntFld = mpp.distribute(allValues);
break;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -112,11 +112,13 @@ void Foam::mappedFlowRateFvPatchVectorField::updateCoeffs()
nbrMesh
).boundary()[mpp.samplePolyPatch().index()];
scalarList phi =
nbrPatch.lookupPatchField<surfaceScalarField, scalar>(nbrPhiName_);
mpp.distribute(phi);
const scalarField phi
(
mpp.distribute
(
nbrPatch.lookupPatchField<surfaceScalarField, scalar>(nbrPhiName_)
)
);
const surfaceScalarField& phiName =
db().lookupObject<surfaceScalarField>(phiName_);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -172,11 +172,8 @@ void Foam::mappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
}
}
mpp.distribute(allUValues);
newUValues.transfer(allUValues);
mpp.distribute(allPhiValues);
newPhiValues.transfer(allPhiValues);
newUValues = mpp.distribute(allUValues);
newPhiValues = mpp.distribute(allPhiValues);
break;
}
@ -186,11 +183,8 @@ void Foam::mappedVelocityFluxFixedValueFvPatchField::updateCoeffs()
const label nbrPatchID =
nbrMesh.boundaryMesh().findPatchID(mpp.samplePatch());
newUValues = UField.boundaryField()[nbrPatchID];
mpp.distribute(newUValues);
newPhiValues = phiField.boundaryField()[nbrPatchID];
mpp.distribute(newPhiValues);
newUValues = mpp.distribute(UField.boundaryField()[nbrPatchID]);
newPhiValues = mpp.distribute(phiField.boundaryField()[nbrPatchID]);
break;
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -200,7 +200,7 @@ void Foam::SurfaceFilmModel<CloudType>::cacheFilmFields
diameterParcelPatch_ =
filmModel.cloudDiameterTrans().boundaryField()[filmPatchi];
filmModel.toPrimary(filmPatchi, diameterParcelPatch_, maxEqOp<scalar>());
filmModel.toPrimary(filmPatchi, diameterParcelPatch_);
UFilmPatch_ = filmModel.U().boundaryField()[filmPatchi];
filmModel.toPrimary(filmPatchi, UFilmPatch_);

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -83,19 +83,19 @@ protected:
// Cached injector fields per film patch
//- Parcel mass / patch face
scalarList massParcelPatch_;
scalarField massParcelPatch_;
//- Parcel diameter / patch face
scalarList diameterParcelPatch_;
scalarField diameterParcelPatch_;
//- Film velocity / patch face
List<vector> UFilmPatch_;
vectorField UFilmPatch_;
//- Film density / patch face
scalarList rhoFilmPatch_;
scalarField rhoFilmPatch_;
//- Film height of all film patches / patch face
scalarListList deltaFilmPatch_;
List<scalarField> deltaFilmPatch_;
// Counters

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -104,10 +104,10 @@ protected:
// Cached injector fields per film patch
//- Film temperature / patch face
scalarList TFilmPatch_;
scalarField TFilmPatch_;
//- Film specific heat capacity / patch face
scalarList CpFilmPatch_;
scalarField CpFilmPatch_;
// Interaction model data

View File

@ -24,24 +24,20 @@ License
\*---------------------------------------------------------------------------*/
#include "mappedPatchBase.H"
#include "addToRunTimeSelectionTable.H"
#include "ListListOps.H"
#include "meshSearchMeshObject.H"
#include "meshTools.H"
#include "OFstream.H"
#include "Random.H"
#include "OBJstream.H"
#include "treeDataFace.H"
#include "treeDataPoint.H"
#include "treeDataCell.H"
#include "indexedOctree.H"
#include "polyMesh.H"
#include "polyPatch.H"
#include "SubField.H"
#include "Time.H"
#include "distributionMap.H"
#include "SubField.H"
#include "triPointRef.H"
#include "syncTools.H"
#include "treeDataCell.H"
#include "DynamicField.H"
#include "RemoteData.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -53,14 +49,13 @@ namespace Foam
const char* Foam::NamedEnum
<
Foam::mappedPatchBase::sampleMode,
5
4
>::names[] =
{
"nearestCell",
"nearestPatchFace",
"nearestPatchFaceAMI",
"nearestFace",
"nearestOnlyCell"
"nearestFace"
};
template<>
@ -77,7 +72,7 @@ namespace Foam
}
const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 5>
const Foam::NamedEnum<Foam::mappedPatchBase::sampleMode, 4>
Foam::mappedPatchBase::sampleModeNames_;
const Foam::NamedEnum<Foam::mappedPatchBase::offsetMode, 3>
@ -86,76 +81,40 @@ const Foam::NamedEnum<Foam::mappedPatchBase::offsetMode, 3>
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::mappedPatchBase::collectSamples
(
pointField& samples,
labelList& patchFaceProcs,
labelList& patchFaces
) const
{
// Collect all sample points and the faces they come from.
{
List<pointField> globalSamples(Pstream::nProcs());
globalSamples[Pstream::myProcNo()] = samplePoints();
Pstream::gatherList(globalSamples);
Pstream::scatterList(globalSamples);
samples = ListListOps::combine<pointField>
(
globalSamples,
accessOp<pointField>()
);
}
{
labelListList globalFaces(Pstream::nProcs());
globalFaces[Pstream::myProcNo()] = identity(patch_.size());
Pstream::gatherList(globalFaces);
Pstream::scatterList(globalFaces);
patchFaces = ListListOps::combine<labelList>
(
globalFaces,
accessOp<labelList>()
);
}
{
labelList nPerProc(Pstream::nProcs());
nPerProc[Pstream::myProcNo()] = patch_.size();
Pstream::gatherList(nPerProc);
Pstream::scatterList(nPerProc);
patchFaceProcs.setSize(patchFaces.size());
label sampleI = 0;
forAll(nPerProc, proci)
{
for (label i = 0; i < nPerProc[proci]; i++)
{
patchFaceProcs[sampleI++] = proci;
}
}
}
}
void Foam::mappedPatchBase::findSamples
(
const sampleMode mode,
const pointField& samples,
labelList& sampleProcs,
const globalIndex& patchGlobalIndex,
labelList& sampleGlobalPatchFaces,
labelList& sampleIndices
) const
{
// Lookup the correct region
const polyMesh& mesh = sampleMesh();
// All the info for nearest. Construct to miss
List<nearInfo> nearest(samples.size());
// Gather the sample points into a single globally indexed list
List<point> allPoints(patchGlobalIndex.size());
{
List<pointField> procSamplePoints(Pstream::nProcs());
procSamplePoints[Pstream::myProcNo()] = this->samplePoints();
Pstream::gatherList(procSamplePoints);
Pstream::scatterList(procSamplePoints);
switch (mode)
forAll(procSamplePoints, proci)
{
forAll(procSamplePoints[proci], procSamplei)
{
allPoints
[
patchGlobalIndex.toGlobal(proci, procSamplei)
] = procSamplePoints[proci][procSamplei];
}
}
}
// Nearest info
List<RemoteData<scalar>> allNearest(patchGlobalIndex.size());
switch (mode_)
{
case NEARESTCELL:
{
@ -163,82 +122,45 @@ void Foam::mappedPatchBase::findSamples
{
FatalErrorInFunction
<< "No need to supply a patch name when in "
<< sampleModeNames_[mode] << " mode." << exit(FatalError);
<< sampleModeNames_[mode_] << " mode." << exit(FatalError);
}
//- Note: face-diagonal decomposition
const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
forAll(samples, sampleI)
forAll(allPoints, alli)
{
const point& sample = samples[sampleI];
const point& p = allPoints[alli];
label celli = tree.findInside(sample);
const label celli = tree.findInside(p);
if (celli == -1)
{
nearest[sampleI].second().first() = Foam::sqr(great);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
else
if (celli != -1)
{
const point& cc = mesh.cellCentres()[celli];
nearest[sampleI].first() = pointIndexHit
(
true,
cc,
celli
);
nearest[sampleI].second().first() = magSqr(cc-sample);
nearest[sampleI].second().second() = Pstream::myProcNo();
allNearest[alli].proci = Pstream::myProcNo();
allNearest[alli].elementi = celli;
allNearest[alli].data = magSqr(cc - p);
}
}
break;
}
case NEARESTONLYCELL:
{
if (samplePatch_.size() && samplePatch_ != "none")
{
FatalErrorInFunction
<< "No need to supply a patch name when in "
<< sampleModeNames_[mode] << " mode." << exit(FatalError);
}
//- Note: face-diagonal decomposition
const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
forAll(samples, sampleI)
{
const point& sample = samples[sampleI];
nearest[sampleI].first() = tree.findNearest(sample, sqr(great));
nearest[sampleI].second().first() = magSqr
(
nearest[sampleI].first().hitPoint()
-sample
);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
break;
}
case NEARESTPATCHFACE:
{
const polyPatch& pp = samplePolyPatch();
if (pp.empty())
{
forAll(samples, sampleI)
forAll(allPoints, alli)
{
nearest[sampleI].second().first() = Foam::sqr(great);
nearest[sampleI].second().second() = Pstream::myProcNo();
allNearest[alli].proci = -1;
allNearest[alli].elementi = -1;
allNearest[alli].data = Foam::sqr(great);
}
}
else
{
// patch faces
const labelList patchFaces(identity(pp.size()) + pp.start());
treeBoundBox patchBb
@ -260,30 +182,20 @@ void Foam::mappedPatchBase::findSamples
3.0 // duplicity
);
forAll(samples, sampleI)
forAll(allPoints, alli)
{
const point& sample = samples[sampleI];
const point& p = allPoints[alli];
pointIndexHit& nearInfo = nearest[sampleI].first();
nearInfo = boundaryTree.findNearest
(
sample,
magSqr(patchBb.span())
);
const pointIndexHit pih =
boundaryTree.findNearest(p, magSqr(patchBb.span()));
if (!nearInfo.hit())
if (pih.hit())
{
nearest[sampleI].second().first() = Foam::sqr(great);
nearest[sampleI].second().second() =
Pstream::myProcNo();
}
else
{
point fc(pp[nearInfo.index()].centre(pp.points()));
nearInfo.setPoint(fc);
nearest[sampleI].second().first() = magSqr(fc-sample);
nearest[sampleI].second().second() =
Pstream::myProcNo();
const point fc = pp[pih.index()].centre(pp.points());
allNearest[alli].proci = Pstream::myProcNo();
allNearest[alli].elementi = pih.index();
allNearest[alli].data = magSqr(fc - p);
}
}
}
@ -296,36 +208,26 @@ void Foam::mappedPatchBase::findSamples
{
FatalErrorInFunction
<< "No need to supply a patch name when in "
<< sampleModeNames_[mode] << " mode." << exit(FatalError);
<< sampleModeNames_[mode_] << " mode." << exit(FatalError);
}
//- Note: face-diagonal decomposition
const meshSearchMeshObject& meshSearchEngine =
meshSearchMeshObject::New(mesh);
forAll(samples, sampleI)
forAll(allPoints, alli)
{
const point& sample = samples[sampleI];
const point& p = allPoints[alli];
label facei = meshSearchEngine.findNearestFace(sample);
const label facei = meshSearchEngine.findNearestFace(p);
if (facei == -1)
{
nearest[sampleI].second().first() = Foam::sqr(great);
nearest[sampleI].second().second() = Pstream::myProcNo();
}
else
if (facei != -1)
{
const point& fc = mesh.faceCentres()[facei];
nearest[sampleI].first() = pointIndexHit
(
true,
fc,
facei
);
nearest[sampleI].second().first() = magSqr(fc-sample);
nearest[sampleI].second().second() = Pstream::myProcNo();
allNearest[alli].proci = Pstream::myProcNo();
allNearest[alli].elementi = facei;
allNearest[alli].data = magSqr(fc - p);
}
}
break;
@ -333,56 +235,96 @@ void Foam::mappedPatchBase::findSamples
case NEARESTPATCHFACEAMI:
{
// nothing to do here
return;
}
default:
{
FatalErrorInFunction
<< "problem." << abort(FatalError);
break;
}
}
// Find nearest. Combine on master.
Pstream::listCombineGather(nearest, nearestEqOp());
Pstream::listCombineScatter(nearest);
Pstream::listCombineGather
(
allNearest,
RemoteData<scalar>::smallestEqOp()
);
Pstream::listCombineScatter(allNearest);
if (debug)
// Determine the number of missing samples (no reduction necessary as this
// is computed from synchronised data)
label nNotFound = 0;
forAll(allPoints, alli)
{
InfoInFunction
<< "mesh " << sampleRegion() << " : " << endl;
forAll(nearest, sampleI)
if (allNearest[alli].proci == -1)
{
label proci = nearest[sampleI].second().second();
label localI = nearest[sampleI].first().index();
Info<< " " << sampleI << " coord:"<< samples[sampleI]
<< " found on processor:" << proci
<< " in local cell/face/point:" << localI
<< " with location:" << nearest[sampleI].first().rawPoint()
<< endl;
nNotFound ++;
}
}
// Convert back into proc+local index
sampleProcs.setSize(samples.size());
sampleIndices.setSize(samples.size());
forAll(nearest, sampleI)
// If any points were not found within cells then re-search for them using
// a nearest test, which should not fail. Warn that this is happening. If
// any points were not found for some other method, then fail.
if (nNotFound)
{
if (!nearest[sampleI].first().hit())
if (mode_ == NEARESTCELL)
{
sampleProcs[sampleI] = -1;
sampleIndices[sampleI] = -1;
WarningInFunction
<< "Did not find " << nNotFound
<< " out of " << returnReduce(patch_.size(), sumOp<label>())
<< " total samples. Sampling these on the nearest cell"
<< " centre instead." << endl << "On patch " << patch_.name()
<< " on region " << sampleRegion()
<< " with sample mode " << sampleModeNames_[mode_] << endl
<< "and offset mode " << offsetModeNames_[offsetMode_] << "."
<< endl;
//- Note: face-diagonal decomposition
const indexedOctree<Foam::treeDataCell>& tree = mesh.cellTree();
forAll(allPoints, alli)
{
const point& p = allPoints[alli];
if (allNearest[alli].proci == -1)
{
const pointIndexHit pih = tree.findNearest(p, sqr(great));
allNearest[alli].proci = Pstream::myProcNo();
allNearest[alli].elementi = pih.index();
allNearest[alli].data = magSqr(pih.hitPoint() - p);
}
}
Pstream::listCombineGather
(
allNearest,
RemoteData<scalar>::smallestEqOp()
);
Pstream::listCombineScatter(allNearest);
}
else
{
sampleProcs[sampleI] = nearest[sampleI].second().second();
sampleIndices[sampleI] = nearest[sampleI].first().index();
FatalErrorInFunction
<< "Mapping failed for " << nl
<< " patch: " << patch_.name() << nl
<< " sampleRegion: " << sampleRegion() << nl
<< " samplePatch: " << samplePatch() << nl
<< " sampleMode: " << sampleModeNames_[mode_] << nl
<< " offsetMode: " << offsetModeNames_[offsetMode_]
<< exit(FatalError);
}
}
// Build lists of samples
DynamicList<label> samplePatchGlobalFacesDyn;
DynamicList<label> sampleIndicesDyn;
forAll(allNearest, alli)
{
if (allNearest[alli].proci == Pstream::myProcNo())
{
samplePatchGlobalFacesDyn.append(alli);
sampleIndicesDyn.append(allNearest[alli].elementi);
}
}
sampleGlobalPatchFaces.transfer(samplePatchGlobalFacesDyn);
sampleIndices.transfer(sampleIndicesDyn);
}
@ -390,10 +332,6 @@ Foam::label Foam::mappedPatchBase::sampleSize() const
{
switch (mode_)
{
case NEARESTPATCHFACEAMI:
{
return samplePolyPatch().size();
}
case NEARESTCELL:
{
return sampleMesh().nCells();
@ -402,17 +340,17 @@ Foam::label Foam::mappedPatchBase::sampleSize() const
{
return samplePolyPatch().size();
}
case NEARESTPATCHFACEAMI:
{
return samplePolyPatch().size();
}
case NEARESTFACE:
{
return sampleMesh().nFaces() - sampleMesh().nInternalFaces();
}
default:
{
FatalErrorInFunction
<< "problem." << abort(FatalError);
return -1;
}
}
return -1;
}
@ -443,140 +381,73 @@ void Foam::mappedPatchBase::calcMapping() const
}
*/
// Get global list of all samples and the processor and face they come from.
pointField samples;
labelList patchFaceProcs;
labelList patchFaces;
collectSamples(samples, patchFaceProcs, patchFaces);
const globalIndex patchGlobalIndex(patch_.size());
// Find processor and cell/face samples are in and actual location.
labelList sampleProcs;
labelList sampleIndices;
findSamples(mode_, samples, sampleProcs, sampleIndices);
// Find processor and cell/face indices of samples
labelList sampleGlobalPatchFaces, sampleIndices;
findSamples(patchGlobalIndex, sampleGlobalPatchFaces, sampleIndices);
// Check for samples that were not found. This will only happen for
// NEARESTCELL since this finds a cell containing a location.
if (mode_ == NEARESTCELL)
{
const label nNotFound =
returnReduce(count(sampleProcs, -1), sumOp<label>());
if (nNotFound > 0)
{
WarningInFunction
<< "Did not find " << nNotFound
<< " out of " << sampleProcs.size() << " total samples."
<< " Sampling these on the nearest cell centre instead." << endl
<< "On patch " << patch_.name()
<< " on region " << sampleRegion()
<< " with sample mode " << sampleModeNames_[mode_] << endl
<< "and offset mode " << offsetModeNames_[offsetMode_] << "."
<< endl;
// Collect the samples that cannot be found
DynamicList<label> subMap;
DynamicField<point> subSamples;
forAll(sampleProcs, sampleI)
{
if (sampleProcs[sampleI] == -1)
{
subMap.append(sampleI);
subSamples.append(samples[sampleI]);
}
}
// And re-search for pure nearest (should not fail)
labelList subSampleProcs;
labelList subSampleIndices;
findSamples
(
NEARESTONLYCELL,
subSamples,
subSampleProcs,
subSampleIndices
);
// Insert
UIndirectList<label>(sampleProcs, subMap) = subSampleProcs;
UIndirectList<label>(sampleIndices, subMap) = subSampleIndices;
}
}
// Check for mapping failure
if (findIndex(sampleProcs, -1) != -1)
{
FatalErrorInFunction
<< "Mapping failed for " << nl
<< " patch: " << patch_.name() << nl
<< " sampleRegion: " << sampleRegion() << nl
<< " samplePatch: " << samplePatch() << nl
<< " sampleMode: " << sampleModeNames_[mode_] << nl
<< " offsetMode: " << offsetModeNames_[offsetMode_]
<< exit(FatalError);
}
// Determine schedule.
mapPtr_.reset(new distributionMap(sampleProcs, patchFaceProcs));
// Rework the schedule from indices into samples to cell data to send,
// face data to receive.
labelListList& subMap = mapPtr_().subMap();
labelListList& constructMap = mapPtr_().constructMap();
forAll(subMap, proci)
{
subMap[proci] = UIndirectList<label>
// Construct distribution schedule
List<Map<label>> compactMap;
mapPtr_.reset
(
new distributionMap
(
sampleIndices,
subMap[proci]
);
constructMap[proci] = UIndirectList<label>
patchGlobalIndex,
sampleGlobalPatchFaces,
compactMap
)
);
const labelList oldToNew(move(sampleGlobalPatchFaces));
const labelList oldSampleIndices(move(sampleIndices));
// Construct input mapping for data to be distributed
mapIndices_ = labelList(mapPtr_->constructSize(), -1);
UIndirectList<label>(mapIndices_, oldToNew) = oldSampleIndices;
// Reverse the map. This means the map is "forward" when going from samples
// to the patch, which is logical.
mapPtr_.reset
(
new distributionMap
(
patchFaces,
constructMap[proci]
);
}
patch_.size(),
move(mapPtr_->constructMap()),
move(mapPtr_->subMap())
)
);
// Redo constructSize
mapPtr_().constructSize() = patch_.size();
if (debug)
// Dump connecting lines for debugging
if (debug && mode_ == NEARESTCELL)
{
// Check that all elements get a value.
PackedBoolList used(patch_.size(), false);
OBJstream obj
(
patch_.name()
+ "_processor"
+ name(Pstream::myProcNo())
+ ".obj"
);
forAll(constructMap, proci)
// Samples -> Patch
pointField ccs(sampleMesh().cellCentres(), mapIndices_);
mapPtr_->distribute(ccs);
forAll(patch_, patchFacei)
{
const labelList& map = constructMap[proci];
forAll(map, i)
{
label facei = map[i];
if (!used[facei])
{
used[facei] = true;
}
else
{
FatalErrorInFunction
<< "On patch " << patch_.name()
<< " patchface " << facei
<< " is assigned to more than once."
<< abort(FatalError);
}
}
const point& fc = patch_.faceCentres()[patchFacei];
const point mid = 0.51*fc + 0.49*ccs[patchFacei];
obj.write(linePointRef(fc, mid));
}
forAll(used, facei)
// Patch -> Samples
pointField fcs(patch_.faceCentres());
mapPtr_->reverseDistribute(mapIndices_.size(), fcs);
forAll(mapIndices_, samplei)
{
if (used[facei] == 0)
{
FatalErrorInFunction
<< "On patch " << patch_.name()
<< " patchface " << facei
<< " is never assigned to."
<< abort(FatalError);
}
const label celli = mapIndices_[samplei];
if (celli == -1) continue;
const point& cc = sampleMesh().cellCentres()[celli];
const point mid = 0.51*cc + 0.49*fcs[samplei];
obj.write(linePointRef(cc, mid));
}
}
}
@ -592,33 +463,6 @@ void Foam::mappedPatchBase::calcAMI() const
AMIPtr_.clear();
if (debug)
{
const polyPatch& nbr = samplePolyPatch();
pointField nbrPoints(nbr.localPoints());
OFstream os(patch_.name() + "_neighbourPatch-org.obj");
meshTools::writeOBJ(os, samplePolyPatch().localFaces(), nbrPoints);
// transform neighbour patch to local system
primitivePatch nbrPatch0
(
SubList<face>
(
nbr.localFaces(),
nbr.size()
),
nbrPoints
);
OFstream osN(patch_.name() + "_neighbourPatch-trans.obj");
meshTools::writeOBJ(osN, nbrPatch0, nbrPoints);
OFstream osO(patch_.name() + "_ownerPatch.obj");
meshTools::writeOBJ(osO, patch_.localFaces(), patch_.localPoints());
}
// Get the projection surface, if any
const word surfType(surfDict_.lookupOrDefault<word>("type", "none"));
if (!surfPtr_.valid() && surfType != "none")
@ -741,6 +585,7 @@ Foam::mappedPatchBase::mappedPatchBase
:
patch_(pp),
sampleRegion_(patch_.boundaryMesh().mesh().name()),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mode_(NEARESTPATCHFACE),
samplePatch_(""),
coupleGroup_(),
@ -748,8 +593,8 @@ Foam::mappedPatchBase::mappedPatchBase
offset_(Zero),
offsets_(pp.size(), offset_),
distance_(0),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mapPtr_(nullptr),
mapIndices_(),
AMIPtr_(nullptr),
AMIReverse_(false),
surfPtr_(nullptr),
@ -768,6 +613,7 @@ Foam::mappedPatchBase::mappedPatchBase
:
patch_(pp),
sampleRegion_(sampleRegion),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mode_(mode),
samplePatch_(samplePatch),
coupleGroup_(),
@ -775,8 +621,8 @@ Foam::mappedPatchBase::mappedPatchBase
offset_(Zero),
offsets_(offsets),
distance_(0),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mapPtr_(nullptr),
mapIndices_(),
AMIPtr_(nullptr),
AMIReverse_(false),
surfPtr_(nullptr),
@ -795,6 +641,7 @@ Foam::mappedPatchBase::mappedPatchBase
:
patch_(pp),
sampleRegion_(sampleRegion),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mode_(mode),
samplePatch_(samplePatch),
coupleGroup_(),
@ -802,8 +649,8 @@ Foam::mappedPatchBase::mappedPatchBase
offset_(offset),
offsets_(0),
distance_(0),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mapPtr_(nullptr),
mapIndices_(),
AMIPtr_(nullptr),
AMIReverse_(false),
surfPtr_(nullptr),
@ -822,6 +669,7 @@ Foam::mappedPatchBase::mappedPatchBase
:
patch_(pp),
sampleRegion_(sampleRegion),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mode_(mode),
samplePatch_(samplePatch),
coupleGroup_(),
@ -829,8 +677,8 @@ Foam::mappedPatchBase::mappedPatchBase
offset_(Zero),
offsets_(0),
distance_(distance),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mapPtr_(nullptr),
mapIndices_(),
AMIPtr_(nullptr),
AMIReverse_(false),
surfPtr_(nullptr),
@ -846,6 +694,7 @@ Foam::mappedPatchBase::mappedPatchBase
:
patch_(pp),
sampleRegion_(dict.lookupOrDefault<word>("sampleRegion", "")),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mode_(sampleModeNames_.read(dict.lookup("sampleMode"))),
samplePatch_(dict.lookupOrDefault<word>("samplePatch", "")),
coupleGroup_(dict),
@ -853,8 +702,8 @@ Foam::mappedPatchBase::mappedPatchBase
offset_(Zero),
offsets_(0),
distance_(0.0),
sameRegion_(sampleRegion_ == patch_.boundaryMesh().mesh().name()),
mapPtr_(nullptr),
mapIndices_(),
AMIPtr_(nullptr),
AMIReverse_(dict.lookupOrDefault<bool>("flipNormals", false)),
surfPtr_(nullptr),
@ -923,6 +772,7 @@ Foam::mappedPatchBase::mappedPatchBase
:
patch_(pp),
sampleRegion_(mpb.sampleRegion_),
sameRegion_(mpb.sameRegion_),
mode_(mpb.mode_),
samplePatch_(mpb.samplePatch_),
coupleGroup_(mpb.coupleGroup_),
@ -930,8 +780,8 @@ Foam::mappedPatchBase::mappedPatchBase
offset_(mpb.offset_),
offsets_(mpb.offsets_),
distance_(mpb.distance_),
sameRegion_(mpb.sameRegion_),
mapPtr_(nullptr),
mapIndices_(),
AMIPtr_(nullptr),
AMIReverse_(mpb.AMIReverse_),
surfPtr_(nullptr),
@ -948,6 +798,7 @@ Foam::mappedPatchBase::mappedPatchBase
:
patch_(pp),
sampleRegion_(mpb.sampleRegion_),
sameRegion_(mpb.sameRegion_),
mode_(mpb.mode_),
samplePatch_(mpb.samplePatch_),
coupleGroup_(mpb.coupleGroup_),
@ -960,8 +811,8 @@ Foam::mappedPatchBase::mappedPatchBase
: vectorField(0)
),
distance_(mpb.distance_),
sameRegion_(mpb.sameRegion_),
mapPtr_(nullptr),
mapIndices_(),
AMIPtr_(nullptr),
AMIReverse_(mpb.AMIReverse_),
surfPtr_(nullptr),
@ -1091,6 +942,7 @@ Foam::tmp<Foam::pointField> Foam::mappedPatchBase::samplePoints() const
void Foam::mappedPatchBase::clearOut()
{
mapPtr_.clear();
mapIndices_.clear();
AMIPtr_.clear();
surfPtr_.clear();
}

View File

@ -65,10 +65,6 @@ Description
Note: if offsetMode is \c normal it uses outwards pointing normals. So
supply a negative distance if sampling inside the domain.
Note:
Storage is not optimal. It temporary collects all (patch)face centres
on all processors to keep the addressing calculation simple.
SourceFiles
mappedPatchBase.C
@ -78,8 +74,6 @@ SourceFiles
#define mappedPatchBase_H
#include "pointField.H"
#include "Tuple2.H"
#include "pointIndexHit.H"
#include "AMIInterpolation.H"
#include "coupleGroupIdentifier.H"
@ -109,8 +103,7 @@ public:
NEARESTCELL, // nearest cell containing sample
NEARESTPATCHFACE, // nearest face on selected patch
NEARESTPATCHFACEAMI, // nearest patch face + AMI interpolation
NEARESTFACE, // nearest face
NEARESTONLYCELL // nearest cell (even if not containing cell)
NEARESTFACE // nearest face
};
//- How to project face centres
@ -121,37 +114,11 @@ public:
NORMAL // use face normal + distance
};
static const NamedEnum<sampleMode, 5> sampleModeNames_;
static const NamedEnum<sampleMode, 4> sampleModeNames_;
static const NamedEnum<offsetMode, 3> offsetModeNames_;
//- Helper class for finding nearest. Point and local index, squared
// distance, and processor index.
typedef Tuple2<pointIndexHit, Tuple2<scalar, label>> nearInfo;
//- Equality operator for nearest info
class nearestEqOp
{
public:
void operator()(nearInfo& x, const nearInfo& y) const
{
if (y.first().hit())
{
if (!x.first().hit())
{
x = y;
}
else if (y.second().first() < x.second().first())
{
x = y;
}
}
}
};
protected:
// Protected data
@ -162,6 +129,9 @@ protected:
//- Region to sample
mutable word sampleRegion_;
//- Same region
mutable bool sameRegion_;
//- What to sample
const sampleMode mode_;
@ -183,15 +153,15 @@ protected:
//- Offset distance (normal)
scalar distance_;
//- Same region
mutable bool sameRegion_;
// Derived information
//- Distributor
mutable autoPtr<distributionMap> mapPtr_;
//- Map pre-addressing
mutable labelList mapIndices_;
// AMI interpolator (only for NEARESTPATCHFACEAMI)
@ -210,21 +180,12 @@ protected:
// Protected Member Functions
//- Collect single list of samples and originating processor+face.
void collectSamples
(
pointField& samples,
labelList& patchFaceProcs,
labelList& patchFaces
) const;
//- Find cells/faces containing samples
void findSamples
(
const sampleMode mode, // search mode
const pointField& samples,
labelList& sampleProcs, // processor containing sample
labelList& sampleIndices // local index of cell/face
const globalIndex& patchGlobalIndex,
labelList& sampleGlobalPatchFaces,
labelList& sampleIndices
) const;
//- Return size of mapped mesh/patch/boundary
@ -324,6 +285,10 @@ public:
//- Return reference to the parallel distribution map
inline const distributionMap& map() const;
//- Return reference to the indices that have to be supplied to the
// parallel distribution map
inline const labelList& mapIndices() const;
//- Get the region mesh
const polyMesh& sampleMesh() const;
@ -345,19 +310,38 @@ public:
//- Wrapper around map/interpolate data distribution
template<class Type>
void distribute(List<Type>& lst) const;
//- Wrapper around map/interpolate data distribution with operation
template<class Type, class CombineOp>
void distribute(List<Type>& lst, const CombineOp& cop) const;
tmp<Field<Type>> distribute(const Field<Type>& fld) const;
//- Wrapper around map/interpolate data distribution
template<class Type>
void reverseDistribute(List<Type>& lst) const;
tmp<Field<Type>> distribute(const tmp<Field<Type>>& fld) const;
//- Wrapper around map/interpolate data distribution
template<class Type>
tmp<Field<Type>> reverseDistribute(const Field<Type>& fld) const;
//- Wrapper around map/interpolate data distribution
template<class Type>
tmp<Field<Type>> reverseDistribute
(
const tmp<Field<Type>>& fld
) const;
//- Wrapper around map/interpolate data distribution with operation
template<class Type, class CombineOp>
void reverseDistribute(List<Type>& lst, const CombineOp& cop) const;
tmp<Field<Type>> reverseDistribute
(
const Field<Type>& fld,
const CombineOp& cop
) const;
//- Wrapper around map/interpolate data distribution with operation
template<class Type, class CombineOp>
tmp<Field<Type>> reverseDistribute
(
const tmp<Field<Type>>& fld,
const CombineOp& cop
) const;
// I/O

View File

@ -51,6 +51,8 @@ inline const Foam::word& Foam::mappedPatchBase::sampleRegion() const
const label samplePatchID =
coupleGroup_.findOtherPatchID(patch_, sampleRegion_);
sameRegion_ = sampleRegion_ == patch_.boundaryMesh().mesh().name();
samplePatch_ = sampleMesh().boundaryMesh()[samplePatchID].name();
}
@ -62,20 +64,7 @@ inline const Foam::word& Foam::mappedPatchBase::samplePatch() const
{
if (samplePatch_.empty())
{
if (!coupleGroup_.valid())
{
FatalErrorInFunction
<< "Supply either a patchName or a coupleGroup"
<< " for patch " << patch_.name()
<< " in region " << patch_.boundaryMesh().mesh().name()
<< exit(FatalError);
}
// Try and use patchGroup to find samplePatch and sampleRegion
const label samplePatchID =
coupleGroup_.findOtherPatchID(patch_, sampleRegion_);
samplePatch_ = sampleMesh().boundaryMesh()[samplePatchID].name();
sampleRegion();
}
return samplePatch_;
@ -99,4 +88,15 @@ inline const Foam::distributionMap& Foam::mappedPatchBase::map() const
}
inline const Foam::labelList& Foam::mappedPatchBase::mapIndices() const
{
if (mapPtr_.empty())
{
calcMapping();
}
return mapIndices_;
}
// ************************************************************************* //

View File

@ -25,11 +25,11 @@ License
#include "mappedPatchBase.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::mappedPatchBase::distribute(List<Type>& lst) const
Foam::tmp<Foam::Field<Type>>
Foam::mappedPatchBase::distribute(const Field<Type>& fld) const
{
switch (mode_)
{
@ -39,116 +39,77 @@ void Foam::mappedPatchBase::distribute(List<Type>& lst) const
{
calcAMI();
}
lst = AMIPtr_->interpolateToSource(Field<Type>(move(lst)));
break;
return AMIPtr_->interpolateToSource(fld);
}
default:
{
map().distribute(lst);
}
}
}
template<class Type, class CombineOp>
void Foam::mappedPatchBase::distribute
(
List<Type>& lst,
const CombineOp& cop
) const
{
switch (mode_)
{
case NEARESTPATCHFACEAMI:
{
if (AMIPtr_.empty())
if (mapPtr_.empty())
{
calcAMI();
calcMapping();
}
lst = AMIPtr_->interpolateToSource(Field<Type>(move(lst)), cop);
break;
}
default:
{
distributionMapBase::distribute
(
Pstream::defaultCommsType,
map().schedule(),
map().constructSize(),
map().subMap(),
false,
map().constructMap(),
false,
lst,
cop,
flipOp(),
Type(Zero)
);
tmp<Field<Type>> tResult(new Field<Type>(fld, mapIndices_));
mapPtr_->distribute(tResult.ref());
return tResult;
}
}
}
template<class Type>
void Foam::mappedPatchBase::reverseDistribute(List<Type>& lst) const
Foam::tmp<Foam::Field<Type>>
Foam::mappedPatchBase::distribute(const tmp<Field<Type>>& fld) const
{
tmp<Field<Type>> tResult = distribute(fld());
fld.clear();
return tResult;
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::mappedPatchBase::reverseDistribute(const Field<Type>& fld) const
{
switch (mode_)
{
case NEARESTPATCHFACE:
{
if (mapPtr_.empty())
{
calcMapping();
}
Field<Type> sampleFld(fld);
mapPtr_->reverseDistribute(mapIndices_.size(), sampleFld);
tmp<Field<Type>> tResult(new Field<Type>(samplePolyPatch().size()));
UIndirectList<Type>(tResult.ref(), mapIndices_) = fld;
return tResult;
}
case NEARESTPATCHFACEAMI:
{
if (AMIPtr_.empty())
{
calcAMI();
}
lst = AMIPtr_->interpolateToTarget(Field<Type>(move(lst)));
break;
return AMIPtr_->interpolateToTarget(fld);
}
default:
{
map().reverseDistribute(sampleSize(), lst);
break;
FatalErrorInFunction
<< "Reverse distribute can only be used in "
<< sampleModeNames_[NEARESTPATCHFACE] << " or "
<< sampleModeNames_[NEARESTPATCHFACEAMI] << "mode"
<< exit(FatalError);
}
}
}
template<class Type, class CombineOp>
void Foam::mappedPatchBase::reverseDistribute
(
List<Type>& lst,
const CombineOp& cop
) const
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::mappedPatchBase::reverseDistribute(const tmp<Field<Type>>& fld) const
{
switch (mode_)
{
case NEARESTPATCHFACEAMI:
{
if (AMIPtr_.empty())
{
calcAMI();
}
lst = AMIPtr_->interpolateToTarget(Field<Type>(move(lst)), cop);
break;
}
default:
{
distributionMapBase::distribute
(
Pstream::defaultCommsType,
map().schedule(),
sampleSize(),
map().constructMap(),
false,
map().subMap(),
false,
lst,
cop,
flipOp(),
Type(Zero)
);
break;
}
}
tmp<Field<Type>> tResult = reverseDistribute(fld());
fld.clear();
return tResult;
}

View File

@ -28,10 +28,6 @@ Description
Determines a mapping between patch face centres and mesh cell or face
centres and processors they're on.
Note:
Storage is not optimal. It stores all face centres and cells on all
processors to keep the addressing calculation simple.
SourceFiles
mappedPolyPatch.C

View File

@ -28,10 +28,6 @@ Description
Determines a mapping between patch face centres and mesh cell or face
centres and processors they're on.
Note:
Storage is not optimal. It stores all face centres and cells on all
processors to keep the addressing calculation simple.
SourceFiles
mappedWallPolyPatch.C

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2012-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,7 +26,7 @@ License
#include "regionToFace.H"
#include "polyMesh.H"
#include "faceSet.H"
#include "mappedPatchBase.H"
#include "RemoteData.H"
#include "indirectPrimitivePatch.H"
#include "PatchTools.H"
#include "addToRunTimeSelectionTable.H"
@ -108,45 +108,36 @@ void Foam::regionToFace::combine(topoSet& set, const bool add) const
mesh_.points()
);
mappedPatchBase::nearInfo ni
(
pointIndexHit(false, Zero, -1),
Tuple2<scalar, label>
(
sqr(great),
Pstream::myProcNo()
)
);
RemoteData<scalar> ni;
forAll(patch, i)
{
const point& fc = patch.faceCentres()[i];
scalar d2 = magSqr(fc-nearPoint_);
if (!ni.first().hit() || d2 < ni.second().first())
const scalar dSqr = magSqr(fc - nearPoint_);
if (ni.proci == -1 || dSqr < ni.data)
{
ni.second().first() = d2;
ni.first().setHit();
ni.first().setPoint(fc);
ni.first().setIndex(i);
ni.proci = Pstream::myProcNo();
ni.elementi = i;
ni.data = dSqr;
}
}
// Globally reduce
combineReduce(ni, mappedPatchBase::nearestEqOp());
combineReduce(ni, RemoteData<scalar>::smallestEqOp());
Info<< " Found nearest face at " << ni.first().rawPoint()
<< " on processor " << ni.second().second()
<< " face " << ni.first().index()
<< " distance " << Foam::sqrt(ni.second().first()) << endl;
Info<< " Found nearest face on processor " << ni.proci
<< " face " << ni.elementi
<< " distance " << Foam::sqrt(ni.data) << endl;
labelList faceRegion(patch.size(), -1);
markZone
(
patch,
ni.second().second(), // proci
ni.first().index(), // start face
0, // currentZone
ni.proci,
ni.elementi,
0,
faceRegion
);

View File

@ -142,7 +142,7 @@ Foam::radiationCoupledBase::~radiationCoupledBase()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalarField Foam::radiationCoupledBase::emissivity() const
Foam::tmp<Foam::scalarField> Foam::radiationCoupledBase::emissivity() const
{
switch (method_)
{
@ -168,16 +168,14 @@ Foam::scalarField Foam::radiationCoupledBase::emissivity() const
// NOTE: for an opaqueSolid the absorptionEmission model returns the
// emissivity of the surface rather than the emission coefficient
// and the input specification MUST correspond to this.
scalarField emissivity
(
radiation.absorptionEmission().e()().boundaryField()
[
nbrPatch.index()
]
);
mpp.distribute(emissivity);
return emissivity;
return
mpp.distribute
(
radiation.absorptionEmission().e()().boundaryField()
[
nbrPatch.index()
]
);
}
break;

View File

@ -137,7 +137,7 @@ public:
//- Calculate corresponding emissivity field
scalarField emissivity() const;
tmp<scalarField> emissivity() const;
// Mapping functions

View File

@ -225,7 +225,7 @@ public:
void toPrimary
(
const label regionPatchi,
List<Type>& regionField
Field<Type>& regionField
) const;
//- Convert a local region field to the primary region with op
@ -233,7 +233,7 @@ public:
void toPrimary
(
const label regionPatchi,
List<Type>& regionField,
Field<Type>& regionField,
const CombineOp& cop
) const;
@ -242,7 +242,7 @@ public:
void toRegion
(
const label regionPatchi,
List<Type>& primaryFieldField
Field<Type>& primaryFieldField
) const;
//- Return a primary patch field mapped the region internal field

View File

@ -31,7 +31,7 @@ template<class Type>
void Foam::regionModels::regionModel::toPrimary
(
const label regionPatchi,
List<Type>& regionField
Field<Type>& regionField
) const
{
forAll(intCoupledPatchIDs_, i)
@ -43,35 +43,7 @@ void Foam::regionModels::regionModel::toPrimary
(
regionMesh().boundaryMesh()[regionPatchi]
);
mpb.reverseDistribute(regionField);
return;
}
}
FatalErrorInFunction
<< "Region patch ID " << regionPatchi << " not found in region mesh"
<< abort(FatalError);
}
template<class Type, class CombineOp>
void Foam::regionModels::regionModel::toPrimary
(
const label regionPatchi,
List<Type>& regionField,
const CombineOp& cop
) const
{
forAll(intCoupledPatchIDs_, i)
{
if (intCoupledPatchIDs_[i] == regionPatchi)
{
const mappedPatchBase& mpb =
refCast<const mappedPatchBase>
(
regionMesh().boundaryMesh()[regionPatchi]
);
mpb.reverseDistribute(regionField, cop);
regionField = mpb.reverseDistribute(regionField);
return;
}
}
@ -86,7 +58,7 @@ template<class Type>
void Foam::regionModels::regionModel::toRegion
(
const label regionPatchi,
List<Type>& primaryField
Field<Type>& primaryField
) const
{
forAll(intCoupledPatchIDs_, i)
@ -98,7 +70,7 @@ void Foam::regionModels::regionModel::toRegion
(
regionMesh().boundaryMesh()[regionPatchi]
);
mpb.distribute(primaryField);
primaryField = mpb.distribute(primaryField);
return;
}
}
@ -122,8 +94,8 @@ void Foam::regionModels::regionModel::toRegion
mappedPatchFieldBase<Type> mpf(mpb, primaryPatchField);
UIndirectList<Type>(regionField, regionPatch.faceCells())
= mpf.mappedField();
UIndirectList<Type>(regionField, regionPatch.faceCells()) =
mpf.mappedField();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,7 +26,7 @@ License
#include "patchProbes.H"
#include "volFields.H"
#include "IOmanip.H"
#include "mappedPatchBase.H"
#include "RemoteData.H"
#include "treeBoundBox.H"
#include "treeDataFace.H"
#include "addToRunTimeSelectionTable.H"
@ -65,7 +65,7 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
}
// All the info for nearest. Construct to miss
List<mappedPatchBase::nearInfo> nearest(this->size());
List<RemoteData<scalar>> nearest(this->size());
const polyPatch& pp = bm[patchi];
@ -94,7 +94,6 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
3.0 // duplicity
);
forAll(probeLocations(), probei)
{
const point sample = probeLocations()[probei];
@ -133,26 +132,19 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
{
const point& fc = mesh.faceCentres()[facei];
mappedPatchBase::nearInfo sampleInfo;
sampleInfo.first() = pointIndexHit
(
true,
fc,
facei
);
sampleInfo.second().first() = magSqr(fc-sample);
sampleInfo.second().second() = Pstream::myProcNo();
nearest[probei]= sampleInfo;
nearest[probei].proci = Pstream::myProcNo();
nearest[probei].elementi = facei;
nearest[probei].data = magSqr(fc-sample);
}
}
}
// Find nearest.
Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
Pstream::listCombineGather
(
nearest,
RemoteData<scalar>::smallestEqOp()
);
Pstream::listCombineScatter(nearest);
if (debug)
@ -160,26 +152,20 @@ void Foam::patchProbes::findElements(const fvMesh& mesh)
InfoInFunction << endl;
forAll(nearest, sampleI)
{
label proci = nearest[sampleI].second().second();
label localI = nearest[sampleI].first().index();
Info<< " " << sampleI << " coord:"<< operator[](sampleI)
<< " found on processor:" << proci
<< " in local cell/face:" << localI
<< " with fc:" << nearest[sampleI].first().rawPoint() << endl;
Info<< " " << sampleI << " coord:" << operator[](sampleI)
<< " found on processor:" << nearest[sampleI].proci
<< " in local cell/face:" << nearest[sampleI].elementi
<< endl;
}
}
// Extract any local faces to sample
elementList_.setSize(nearest.size(), -1);
forAll(nearest, sampleI)
{
if (nearest[sampleI].second().second() == Pstream::myProcNo())
if (nearest[sampleI].proci == Pstream::myProcNo())
{
// Store the face to sample
elementList_[sampleI] = nearest[sampleI].first().index();
elementList_[sampleI] = nearest[sampleI].elementi;
}
}
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,7 +30,7 @@ License
#include "treeDataFace.H"
#include "Time.H"
#include "meshTools.H"
#include "mappedPatchBase.H"
#include "RemoteData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -105,42 +105,31 @@ void Foam::sampledSets::boundaryPoints::calcSamples
(void)mesh().tetBasePtIs();
// Generate the nearest patch information for each sampling point
List<mappedPatchBase::nearInfo> nearest(points_.size());
List<RemoteData<Tuple2<scalar, point>>> nearest(points_.size());
forAll(points_, sampleI)
{
const point& sample = points_[sampleI];
pointIndexHit& nearHit = nearest[sampleI].first();
scalar& nearDist = nearest[sampleI].second().first();
label& nearProc = nearest[sampleI].second().second();
const pointIndexHit pih =
patchFaces.size()
? patchTree.findNearest(sample, sqr(maxDistance_))
: pointIndexHit();
// Find the nearest
if (patchFaces.size())
if (pih.hit())
{
nearHit = patchTree.findNearest(sample, sqr(maxDistance_));
}
else
{
nearHit.setMiss();
}
// Fill in the information
if (nearHit.hit())
{
nearHit.setIndex(patchFaces[nearHit.index()]);
nearDist = magSqr(nearHit.hitPoint() - sample);
nearProc = Pstream::myProcNo();
}
else
{
nearHit.setIndex(-1);
nearDist = Foam::sqr(great);
nearProc = Pstream::myProcNo();
nearest[sampleI].proci = Pstream::myProcNo();
nearest[sampleI].elementi = patchFaces[pih.index()];
nearest[sampleI].data.first() = magSqr(pih.hitPoint() - sample);
nearest[sampleI].data.second() = pih.hitPoint();
}
}
// Reduce to get the nearest patch locations globally
Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
Pstream::listCombineGather
(
nearest,
RemoteData<Tuple2<scalar, point>>::smallestFirstEqOp()
);
Pstream::listCombineScatter(nearest);
// Dump connecting lines from the sampling points to the hit locations
@ -150,13 +139,13 @@ void Foam::sampledSets::boundaryPoints::calcSamples
label verti = 0;
forAll(nearest, i)
forAll(nearest, sampleI)
{
if (nearest[i].first().hit())
if (nearest[sampleI].proci != -1)
{
meshTools::writeOBJ(str, points_[i]);
meshTools::writeOBJ(str, points_[sampleI]);
verti++;
meshTools::writeOBJ(str, nearest[i].first().hitPoint());
meshTools::writeOBJ(str, nearest[sampleI].data.second());
verti++;
str << "l " << verti - 1 << ' ' << verti << nl;
}
@ -166,16 +155,13 @@ void Foam::sampledSets::boundaryPoints::calcSamples
// Store the sampling locations on the nearest processor
forAll(nearest, sampleI)
{
const pointIndexHit& nearHit = nearest[sampleI].first();
const label& nearProc = nearest[sampleI].second().second();
if (nearHit.hit())
if (nearest[sampleI].proci != -1)
{
if (nearProc == Pstream::myProcNo())
if (nearest[sampleI].proci == Pstream::myProcNo())
{
label facei = nearHit.index();
const label facei = nearest[sampleI].elementi;
samplingPositions.append(nearHit.hitPoint());
samplingPositions.append(nearest[sampleI].data.second());
samplingSegments.append(sampleI);
samplingCells.append(mesh().faceOwner()[facei]);
samplingFaces.append(facei);

View File

@ -26,6 +26,7 @@ License
#include "sampledPatchInternalField.H"
#include "interpolationCellPoint.H"
#include "PrimitivePatchInterpolation.H"
#include "OBJstream.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -64,6 +65,45 @@ Foam::sampledSurfaces::patchInternalField::sampleField
}
namespace Foam
{
template<class Type>
inline labelList hist(const List<List<Type>>& ll)
{
labelList result;
forAll(ll, i)
{
const label s = ll[i].size();
result.resize(max(result.size(), s + 1), 0);
result[s] ++;
}
return result;
}
template<class Type>
inline List<Type> flatten(const List<List<Type>>& ll)
{
label s = 0;
forAll(ll, i)
{
s += ll[i].size();
}
List<Type> result(s);
label resulti = 0;
forAll(ll, i)
{
SubList<Type>(result, ll[i].size(), resulti) = ll[i];
resulti += ll[i].size();
}
return result;
}
}
template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::sampledSurfaces::patchInternalField::interpolateField
@ -82,35 +122,32 @@ Foam::sampledSurfaces::patchInternalField::interpolateField
forAll(patchIDs(), i)
{
// See mappedFixedValueFvPatchField
// Cells on which samples are generated
const labelList& sampleCells = mappers_[i].mapIndices();
// Send the patch points to the cells
const distributionMap& distMap = mappers_[i].map();
pointField samplePoints(mappers_[i].samplePoints());
distMap.reverseDistribute(sampleCells.size(), samplePoints);
// Send back sample points to processor that holds the cell.
// Mark cells with point::max so we know which ones we need
// to interpolate (since expensive).
vectorField samples(mappers_[i].samplePoints());
distMap.reverseDistribute(mesh().nCells(), point::max, samples);
Field<Type> patchVals(mesh().nCells());
forAll(samples, celli)
// Interpolate values
Field<Type> sampleValues(sampleCells.size());
forAll(sampleCells, i)
{
if (samples[celli] != point::max)
if (sampleCells[i] != -1)
{
patchVals[celli] = interpolator.interpolate
(
samples[celli],
celli
);
sampleValues[i] =
interpolator.interpolate(samplePoints[i], sampleCells[i]);
}
}
distMap.distribute(patchVals);
// Send the values back to the patch
Field<Type> patchValues(sampleValues);
distMap.distribute(patchValues);
// Now patchVals holds the interpolated data in patch face order.
// Collect.
SubList<Type>(allPatchVals, patchVals.size(), sz) = patchVals;
sz += patchVals.size();
// Insert into the full value field
SubList<Type>(allPatchVals, patchValues.size(), sz) = patchValues;
sz += patchValues.size();
}
// Interpolate to points. Reconstruct the patch of all faces to aid

View File

@ -116,12 +116,14 @@ Foam::semiPermeableBaffleMassFractionFvPatchScalarField::calcPhiYp() const
// Get the cell-mass fractions
const scalarField Yc(patchInternalField());
scalarField nbrYc
const scalarField nbrYc
(
nbrPatch.lookupPatchField<volScalarField, scalar>(YName)
.patchInternalField()
mpp.distribute
(
nbrPatch.lookupPatchField<volScalarField, scalar>(YName)
.patchInternalField()
)
);
mpp.distribute(nbrYc);
// Get the patch delta coefficients multiplied by the diffusivity
const thermophysicalTransportModel& ttm =
@ -133,11 +135,13 @@ Foam::semiPermeableBaffleMassFractionFvPatchScalarField::calcPhiYp() const
(
ttm.alphaEff(patch().index())*patch().deltaCoeffs()
);
scalarField nbrAlphaEffDeltap
const scalarField nbrAlphaEffDeltap
(
ttm.alphaEff(nbrPatch.index())*nbrPatch.deltaCoeffs()
mpp.distribute
(
ttm.alphaEff(nbrPatch.index())*nbrPatch.deltaCoeffs()
)
);
mpp.distribute(nbrAlphaEffDeltap);
// Get the specie molecular weight, if needed
scalar Wi = NaN;
@ -152,8 +156,7 @@ Foam::semiPermeableBaffleMassFractionFvPatchScalarField::calcPhiYp() const
if (property_ == moleFraction || property_ == partialPressure)
{
tW = thermo.W(patch().index());
tNbrW = thermo.W(nbrPatch.index());
mpp.distribute(tNbrW.ref());
tNbrW = mpp.distribute(thermo.W(nbrPatch.index()));
}
// Construct coefficients that convert mass fraction to the property that
@ -172,8 +175,8 @@ Foam::semiPermeableBaffleMassFractionFvPatchScalarField::calcPhiYp() const
case molarConcentration:
{
k *= thermo.rho(patch().index())/Wi;
scalarField nbrRhop(thermo.rho(nbrPatch.index()));
mpp.distribute(nbrRhop);
tmp<scalarField> nbrRhop =
mpp.distribute(thermo.rho(nbrPatch.index()));
nbrK *= nbrRhop/Wi;
}
break;
@ -181,8 +184,11 @@ Foam::semiPermeableBaffleMassFractionFvPatchScalarField::calcPhiYp() const
case partialPressure:
{
k *= thermo.p().boundaryField()[patch().index()]*tW/Wi;
scalarField nbrPp(thermo.p().boundaryField()[nbrPatch.index()]);
mpp.distribute(nbrPp);
tmp<scalarField> nbrPp =
mpp.distribute
(
thermo.p().boundaryField()[nbrPatch.index()]
);
nbrK *= nbrPp*tNbrW/Wi;
}
break;