Changed handling of coupled interfaces for ldu

This commit is contained in:
Hrvoje Jasak 2013-05-31 11:10:15 +01:00
parent c6a7e7b70b
commit 542c19da79
72 changed files with 1191 additions and 625 deletions

View file

@ -182,7 +182,8 @@ public:
const lduMatrix&,
const scalarField&,
const direction,
const bool
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{}
@ -194,7 +195,8 @@ public:
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const = 0;
};

View file

@ -234,7 +234,8 @@ public:
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{}
};

View file

@ -48,7 +48,8 @@ template
class Type
>
void ProcessorPointPatchField
<PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::resizeBuf
<PatchField, Mesh, PointPatch, ProcessorPointPatch, MatrixType, Type>::
resizeBuf
(
List<char>& buf,
const label size
@ -961,7 +962,8 @@ initInterfaceMatrixUpdate
const lduMatrix& m,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
tmp<scalarField> tlocalMult(new scalarField(this->size(), 0));
@ -985,6 +987,100 @@ initInterfaceMatrixUpdate
// use the counter.
label coeffI = 0;
if (switchToLhs)
{
// Owner side
// ~~~~~~~~~~
{
const labelList& cutOwn = procPatch_.cutEdgeOwnerIndices();
const labelList& cutOwnStart = procPatch_.cutEdgeOwnerStart();
forAll (mp, pointI)
{
label ownIndex = cutOwnStart[pointI];
label endOwn = cutOwnStart[pointI + 1];
for (; ownIndex < endOwn; ownIndex++)
{
localMult[pointI] +=
coeffs[coeffI]*psiInternal[U[cutOwn[ownIndex]]];
// Multiply the internal side as well using the cut mask
result[U[cutOwn[ownIndex]]] -=
cutMask[coeffI]*coeffs[coeffI]*psiInternal[mp[pointI]];
coeffI++;
}
}
}
// Neighbour side
// ~~~~~~~~~~~~~~
{
const labelList& cutNei = procPatch_.cutEdgeNeighbourIndices();
const labelList& cutNeiStart = procPatch_.cutEdgeNeighbourStart();
forAll (mp, pointI)
{
label neiIndex = cutNeiStart[pointI];
label endNei = cutNeiStart[pointI + 1];
for (; neiIndex < endNei; neiIndex++)
{
localMult[pointI] +=
coeffs[coeffI]*psiInternal[L[cutNei[neiIndex]]];
// Multiply the internal side as well using the cut mask
result[L[cutNei[neiIndex]]] -=
cutMask[coeffI]*coeffs[coeffI]*psiInternal[mp[pointI]];
coeffI++;
}
}
}
// Doubly cut coefficients
// ~~~~~~~~~~~~~~~~~~~~~~~
// There exists a possibility of having an internal edge for a
// point on the processor patch which is in fact connected to
// another point of the same patch. This particular nastiness
// introduces a deformation in the solution because the edge is
// either multiplied twice or not at all. For this purpose, the
// offending edges need to be separated out and multiplied
// appropriately. This will only happen for cell tetrahedral
// decomposition and is generally nasty.
// No need for cut mask here
{
const labelList& doubleCut = procPatch_.doubleCutEdgeIndices();
const labelList& doubleCutOwner = procPatch_.doubleCutOwner();
const labelList& doubleCutNeighbour =
procPatch_.doubleCutNeighbour();
forAll (doubleCut, edgeI)
{
// Owner side
localMult[doubleCutOwner[edgeI]] +=
coeffs[coeffI]*psiInternal[U[doubleCut[edgeI]]];
coeffI++;
// Neighbour side
localMult[doubleCutNeighbour[edgeI]] +=
coeffs[coeffI]*psiInternal[L[doubleCut[edgeI]]];
coeffI++;
}
}
// Add the local multiplication to this side as well
forAll (mp, pointI)
{
result[mp[pointI]] -= localMult[pointI];
}
}
else
{
// Owner side
// ~~~~~~~~~~
{
@ -1051,7 +1147,8 @@ initInterfaceMatrixUpdate
const labelList& doubleCut = procPatch_.doubleCutEdgeIndices();
const labelList& doubleCutOwner = procPatch_.doubleCutOwner();
const labelList& doubleCutNeighbour = procPatch_.doubleCutNeighbour();
const labelList& doubleCutNeighbour =
procPatch_.doubleCutNeighbour();
forAll (doubleCut, edgeI)
{
@ -1073,6 +1170,7 @@ initInterfaceMatrixUpdate
{
result[mp[pointI]] += localMult[pointI];
}
}
// Send the localMult
sendField(tlocalMult, commsType);
@ -1099,9 +1197,13 @@ updateInterfaceMatrix
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Switch to lhs handled in init
// HJ, 22/May/2013
// Get the neighbour side multiplication
tmp<scalarField> tneiMult = receivePointField<scalar>(commsType);
this->addToInternalField(result, tneiMult());

View file

@ -325,7 +325,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -336,7 +337,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
};

View file

@ -1024,7 +1024,8 @@ void GlobalPointPatchField
const lduMatrix& m,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
tmp<scalarField> tlocalMult(new scalarField(this->size(), 0));
@ -1047,11 +1048,142 @@ void GlobalPointPatchField
label coeffI = 0;
scalarField sumOffDiag(this->size(), 0);
if (switchToLhs)
{
// Owner side
// ~~~~~~~~~~
{
const labelList& cutOwn = globalPointPatch_.cutEdgeOwnerIndices();
const labelList& cutOwnStart = globalPointPatch_.cutEdgeOwnerStart();
const labelList& cutOwnStart =
globalPointPatch_.cutEdgeOwnerStart();
forAll (mp, pointI)
{
label ownIndex = cutOwnStart[pointI];
label endOwn = cutOwnStart[pointI + 1];
for (; ownIndex < endOwn; ownIndex++)
{
localMult[pointI] +=
cutMask[coeffI]*coeffs[coeffI]
*psiInternal[U[cutOwn[ownIndex]]];
sumOffDiag[pointI] += cutMask[coeffI]*coeffs[coeffI];
// Multiply the internal side as well
result[U[cutOwn[ownIndex]]] -=
coeffs[coeffI]*psiInternal[mp[pointI]];
coeffI++;
}
}
}
// Neighbour side
// ~~~~~~~~~~~~~~
{
const labelList& cutNei =
globalPointPatch_.cutEdgeNeighbourIndices();
const labelList& cutNeiStart =
globalPointPatch_.cutEdgeNeighbourStart();
forAll (mp, pointI)
{
label neiIndex = cutNeiStart[pointI];
label endNei = cutNeiStart[pointI + 1];
for (; neiIndex < endNei; neiIndex++)
{
localMult[pointI] +=
cutMask[coeffI]*coeffs[coeffI]
*psiInternal[L[cutNei[neiIndex]]];
sumOffDiag[pointI] += cutMask[coeffI]*coeffs[coeffI];
// Multiply the internal side as well
result[L[cutNei[neiIndex]]] -=
coeffs[coeffI]*psiInternal[mp[pointI]];
coeffI++;
}
}
}
// Doubly cut coefficients
// ~~~~~~~~~~~~~~~~~~~~~~~
// There exists a possibility of having an internal edge for a
// point on the processor patch which is in fact connected to
// another point of the same patch. This particular nastiness
// introduces a deformation in the solution because the edge is
// either multiplied twice or not at all. For this purpose, the
// offending edges need to be separated out and multiplied
// appropriately.
{
const labelList& doubleCut =
globalPointPatch_.doubleCutEdgeIndices();
const labelList& doubleCutOwner =
globalPointPatch_.doubleCutOwner();
const labelList& doubleCutNeighbour =
globalPointPatch_.doubleCutNeighbour();
forAll (doubleCut, edgeI)
{
// Owner side
localMult[doubleCutOwner[edgeI]] +=
cutMask[coeffI]*coeffs[coeffI]*
psiInternal[U[doubleCut[edgeI]]];
sumOffDiag[doubleCutOwner[edgeI]] +=
cutMask[coeffI]*coeffs[coeffI];
coeffI++;
// Neighbour side
localMult[doubleCutNeighbour[edgeI]] +=
cutMask[coeffI]*coeffs[coeffI]*
psiInternal[L[doubleCut[edgeI]]];
sumOffDiag[doubleCutNeighbour[edgeI]] +=
cutMask[coeffI]*coeffs[coeffI];
coeffI++;
}
}
// Reduce/extract the result and enforce over all processors
// Requires global sync points to flush buffers before gather-scatter
// communications. Reconsider. HJ, 29/Mar/2011
if (Pstream::defaultCommsType == Pstream::nonBlocking)
{
IPstream::waitRequests();
OPstream::waitRequests();
}
tmp<Field<scalar> > trpf =
reduceExtractPoint<scalar>(localMult);
Field<scalar>& rpf = trpf();
// Get addressing
const labelList& addr = globalPointPatch_.meshPoints();
forAll (addr, i)
{
result[addr[i]] -= rpf[i];
}
}
else
{
// Owner side
// ~~~~~~~~~~
{
const labelList& cutOwn = globalPointPatch_.cutEdgeOwnerIndices();
const labelList& cutOwnStart =
globalPointPatch_.cutEdgeOwnerStart();
forAll (mp, pointI)
{
@ -1078,7 +1210,8 @@ void GlobalPointPatchField
// Neighbour side
// ~~~~~~~~~~~~~~
{
const labelList& cutNei = globalPointPatch_.cutEdgeNeighbourIndices();
const labelList& cutNei =
globalPointPatch_.cutEdgeNeighbourIndices();
const labelList& cutNeiStart =
globalPointPatch_.cutEdgeNeighbourStart();
@ -1115,9 +1248,12 @@ void GlobalPointPatchField
// offending edges need to be separated out and multiplied
// appropriately.
{
const labelList& doubleCut = globalPointPatch_.doubleCutEdgeIndices();
const labelList& doubleCut =
globalPointPatch_.doubleCutEdgeIndices();
const labelList& doubleCutOwner =
globalPointPatch_.doubleCutOwner();
const labelList& doubleCutOwner = globalPointPatch_.doubleCutOwner();
const labelList& doubleCutNeighbour =
globalPointPatch_.doubleCutNeighbour();
@ -1167,6 +1303,7 @@ void GlobalPointPatchField
{
result[addr[i]] += rpf[i];
}
}
}

View file

@ -226,7 +226,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
};

View file

@ -129,7 +129,8 @@ public:
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{}
@ -141,7 +142,8 @@ public:
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{}
@ -155,7 +157,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const = 0;
//- Update result field based on interface functionality
@ -165,7 +168,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const = 0;
};
@ -192,6 +196,7 @@ makeBlockGAMGInterfaceField(blockScalarGAMGInterfaceField, blockGAMGInterfaceFie
makeBlockGAMGInterfaceField(blockVectorGAMGInterfaceField, blockGAMGInterfaceFieldType##Vector); \
makeBlockGAMGInterfaceField(blockTensorGAMGInterfaceField, blockGAMGInterfaceFieldType##Tensor);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

View file

@ -105,7 +105,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const = 0;
//- Update result field based on interface functionality
@ -115,7 +116,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const = 0;
};

View file

@ -49,14 +49,17 @@ Foam::processorBlockGAMGInterfaceField<Type>::processorBlockGAMGInterfaceField
if (isA<processorBlockLduInterfaceField<Type> >(fineInterfaceField))
{
const processorBlockLduInterfaceField<Type>& p =
refCast<const processorBlockLduInterfaceField<Type> >(fineInterfaceField);
refCast<const processorBlockLduInterfaceField<Type> >
(
fineInterfaceField
);
doTransform_ = p.doTransform();
rank_ = p.rank();
}
else
{
FatalErrorIn("Foam::processorBlockGAMGInterfaceField<Type> Constructor")
FatalErrorIn("processorBlockGAMGInterfaceField<Type> Constructor")
<< "fineInterface must be of processor type and either" << endl
<< " processorBlockLduInterfaceField<Type> or " << endl
<< " processorFvPatchField<Type> " << endl
@ -68,7 +71,8 @@ Foam::processorBlockGAMGInterfaceField<Type>::processorBlockGAMGInterfaceField
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
Foam::processorBlockGAMGInterfaceField<Type>::~processorBlockGAMGInterfaceField()
Foam::processorBlockGAMGInterfaceField<Type>::
~processorBlockGAMGInterfaceField()
{}
@ -81,7 +85,8 @@ void Foam::processorBlockGAMGInterfaceField<Type>::initInterfaceMatrixUpdate
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procInterface_.compressedSend
@ -98,7 +103,8 @@ void Foam::processorBlockGAMGInterfaceField<Type>::updateInterfaceMatrix
Field<Type>& result,
const BlockLduMatrix<Type>&,
const CoeffField<Type>& coeffs,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
Field<Type> pnf
@ -108,12 +114,14 @@ void Foam::processorBlockGAMGInterfaceField<Type>::updateInterfaceMatrix
if (coeffs.activeType() == blockCoeffBase::SCALAR)
{
pnf = coeffs.asScalar() *
pnf = coeffs.asScalar()*
procInterface_.compressedReceive<Type>(commsType, this->size())();
}
else if (coeffs.activeType() == blockCoeffBase::LINEAR)
{
pnf = cmptMultiply(coeffs.asLinear(),
pnf = cmptMultiply
(
coeffs.asLinear(),
procInterface_.compressedReceive<Type>(commsType, this->size())()
);
}
@ -125,10 +133,20 @@ void Foam::processorBlockGAMGInterfaceField<Type>::updateInterfaceMatrix
const unallocLabelList& faceCells = procInterface_.faceCells();
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= pnf[elemI];
}
}
}

View file

@ -120,7 +120,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -130,7 +131,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -47,7 +47,8 @@ void processorBlockGAMGInterfaceField<scalar>::initInterfaceMatrixUpdate
Field<scalar>&,
const BlockLduMatrix<scalar>&,
const CoeffField<scalar>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procInterface_.compressedSend
@ -57,6 +58,7 @@ void processorBlockGAMGInterfaceField<scalar>::initInterfaceMatrixUpdate
);
}
template<>
void processorBlockGAMGInterfaceField<scalar>::updateInterfaceMatrix
(
@ -64,7 +66,8 @@ void processorBlockGAMGInterfaceField<scalar>::updateInterfaceMatrix
Field<scalar>& result,
const BlockLduMatrix<scalar>&,
const CoeffField<scalar>& coeffs,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
Field<scalar> pnf
@ -77,12 +80,25 @@ void processorBlockGAMGInterfaceField<scalar>::updateInterfaceMatrix
const unallocLabelList& faceCells = procInterface_.faceCells();
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View file

@ -47,11 +47,15 @@ void processorBlockGAMGInterfaceField<tensor>::initInterfaceMatrixUpdate
Field<tensor>&,
const BlockLduMatrix<tensor>&,
const CoeffField<tensor>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
FatalErrorIn("Foam::tensor Foam::processorBlockGAMGInterfaceField<tensor>::initInterafaceMatrix(...)")
<< "Not implemented" << abort(FatalError);
FatalErrorIn
(
"tensor processorBlockGAMGInterfaceField<tensor>::"
"initInterafaceMatrix(...)"
) << "Not implemented" << abort(FatalError);
}
template<>
@ -61,11 +65,15 @@ void processorBlockGAMGInterfaceField<tensor>::updateInterfaceMatrix
Field<tensor>& result,
const BlockLduMatrix<tensor>&,
const CoeffField<tensor>& coeffs,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
FatalErrorIn("Foam::tensor Foam::processorBlockGAMGInterfaceField<tensor>::initInterafaceMatrix(...)")
<< "Not implemented" << abort(FatalError);
FatalErrorIn
(
"tensor processorBlockGAMGInterfaceField<tensor>::"
"updateInterfaceMatrix(...)"
) << "Not implemented" << abort(FatalError);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -46,7 +46,8 @@ void processorBlockGAMGInterfaceField<vector>::initInterfaceMatrixUpdate
Field<vector>&,
const BlockLduMatrix<vector>&,
const CoeffField<vector>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procInterface_.compressedSend
@ -63,7 +64,8 @@ void processorBlockGAMGInterfaceField<vector>::updateInterfaceMatrix
Field<vector>& result,
const BlockLduMatrix<vector>&,
const CoeffField<vector>& coeffs,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
Field<vector> pnf
@ -76,12 +78,24 @@ void processorBlockGAMGInterfaceField<vector>::updateInterfaceMatrix
const unallocLabelList& faceCells = procInterface_.faceCells();
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= pnf[elemI];
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View file

@ -421,7 +421,8 @@ public:
(
const FieldField<CoeffField, Type>& interfaceCoeffs,
TypeField& result,
const TypeField& psi
const TypeField& psi,
const bool switchToLhs = false
) const;
//- Update coupled interfaces
@ -429,7 +430,8 @@ public:
(
const FieldField<CoeffField, Type>& interfaceCoeffs,
TypeField& result,
const TypeField& psi
const TypeField& psi,
const bool switchToLhs = false
) const;

View file

@ -36,7 +36,8 @@ void Foam::BlockLduMatrix<Type>::initInterfaces
(
const FieldField<CoeffField, Type>& interfaceCoeffs,
TypeField& result,
const TypeField& psi
const TypeField& psi,
const bool switchToLhs
) const
{
if
@ -55,7 +56,8 @@ void Foam::BlockLduMatrix<Type>::initInterfaces
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::defaultCommsType
Pstream::defaultCommsType,
switchToLhs
);
}
}
@ -81,7 +83,8 @@ void Foam::BlockLduMatrix<Type>::initInterfaces
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::blocking
Pstream::blocking,
switchToLhs
);
}
}
@ -101,7 +104,8 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
(
const FieldField<CoeffField, Type>& interfaceCoeffs,
TypeField& result,
const TypeField& psi
const TypeField& psi,
const bool switchToLhs
) const
{
if
@ -127,7 +131,8 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::defaultCommsType
Pstream::defaultCommsType,
switchToLhs
);
}
}
@ -151,7 +156,8 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::scheduled
Pstream::scheduled,
switchToLhs
);
}
else
@ -162,7 +168,8 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::scheduled
Pstream::scheduled,
switchToLhs
);
}
}
@ -185,14 +192,15 @@ void Foam::BlockLduMatrix<Type>::updateInterfaces
result,
*this,
interfaceCoeffs[interfaceI],
Pstream::blocking
Pstream::blocking,
switchToLhs
);
}
}
}
else
{
FatalErrorIn("BlockLduMatrix<Type>::updateMatrixInterfaces")
FatalErrorIn("BlockLduMatrix<Type>::updateInterfaces")
<< "Unsuported communications type "
<< Pstream::commsTypeNames[Pstream::defaultCommsType]
<< exit(FatalError);

View file

@ -32,14 +32,16 @@ Description
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// TODO - This code is currently not called so we have specialized
// initInterfaceMatrixUpdate in processorFvPatchScalarfield. This needs to be fixed
// initInterfaceMatrixUpdate in processorFvPatchScalarfield
// This needs to be fixed
template<>
void Foam::BlockLduMatrix<scalar>::initInterfaces
(
const FieldField<CoeffField, scalar>& coupleCoeffs,
scalarField& result,
const scalarField& psi
const scalarField& psi,
const bool switchToLhs
) const
{
if
@ -58,7 +60,8 @@ void Foam::BlockLduMatrix<scalar>::initInterfaces
result,
*this,
coupleCoeffs[interfaceI].asScalar(),
Pstream::defaultCommsType
Pstream::defaultCommsType,
switchToLhs
);
}
}
@ -84,7 +87,8 @@ void Foam::BlockLduMatrix<scalar>::initInterfaces
result,
*this,
coupleCoeffs[interfaceI].asScalar(),
Pstream::blocking
Pstream::blocking,
switchToLhs
);
}
}
@ -104,7 +108,8 @@ void Foam::BlockLduMatrix<scalar>::updateInterfaces
(
const FieldField<CoeffField, scalar>& coupleCoeffs,
scalarField& result,
const scalarField& psi
const scalarField& psi,
const bool switchToLhs
) const
{
if
@ -130,7 +135,8 @@ void Foam::BlockLduMatrix<scalar>::updateInterfaces
result,
*this,
coupleCoeffs[interfaceI].asScalar(),
Pstream::defaultCommsType
Pstream::defaultCommsType,
switchToLhs
);
}
}
@ -154,7 +160,8 @@ void Foam::BlockLduMatrix<scalar>::updateInterfaces
result,
*this,
coupleCoeffs[interfaceI].asScalar(),
Pstream::scheduled
Pstream::scheduled,
switchToLhs
);
}
else
@ -165,7 +172,8 @@ void Foam::BlockLduMatrix<scalar>::updateInterfaces
result,
*this,
coupleCoeffs[interfaceI].asScalar(),
Pstream::scheduled
Pstream::scheduled,
switchToLhs
);
}
}
@ -188,7 +196,8 @@ void Foam::BlockLduMatrix<scalar>::updateInterfaces
result,
*this,
coupleCoeffs[interfaceI].asScalar(),
Pstream::blocking
Pstream::blocking,
switchToLhs
);
}
}

View file

@ -64,50 +64,24 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
typedef typename TypeCoeffField::linearTypeField linearTypeField;
typedef typename TypeCoeffField::squareTypeField squareTypeField;
FieldField<CoeffField, Type> mBouCoeffs(this->matrix_.coupleUpper().size());
forAll(mBouCoeffs, patchi)
{
const FieldField<CoeffField,Type>& coupleUpperTemp(this->matrix_.coupleUpper());
if (const_cast<BlockLduMatrix<Type>& >(this->matrix_).interfaces().set(patchi))
{
mBouCoeffs.set(patchi, coupleUpperTemp[patchi]);
if (mBouCoeffs[patchi].activeType() == blockCoeffBase::SQUARE)
{
squareTypeField& activeMBouCoeffs = mBouCoeffs[patchi].asSquare();
forAll (activeMBouCoeffs, intI)
{
activeMBouCoeffs[intI] = -activeMBouCoeffs[intI];
}
}
else if (mBouCoeffs[patchi].activeType() == blockCoeffBase::LINEAR)
{
linearTypeField& activeMBouCoeffs = mBouCoeffs[patchi].asLinear();
forAll (activeMBouCoeffs, intI)
{
activeMBouCoeffs[intI] = -activeMBouCoeffs[intI];
}
}
}
}
for (label sweep = 0; sweep < nSweeps_; sweep++)
{
bPrime_ = b;
this->matrix_.initInterfaces
(
mBouCoeffs,
this->matrix_.coupleUpper(),
bPrime_,
x
x,
true // switch to lhs of system
);
this->matrix_.updateInterfaces
(
mBouCoeffs,
this->matrix_.coupleUpper(),
bPrime_,
x
x,
true // switch to lhs of system
);
register label fStart, fEnd, curCoeff;
@ -205,49 +179,24 @@ void Foam::BlockGaussSeidelPrecon<Type>::BlockSweep
typedef typename TypeCoeffField::linearTypeField linearTypeField;
typedef typename TypeCoeffField::squareTypeField squareTypeField;
FieldField<CoeffField, Type> mBouCoeffs(this->matrix_.coupleUpper().size());
forAll(mBouCoeffs, patchi)
{
const FieldField<CoeffField,Type>& coupleUpperTemp(this->matrix_.coupleUpper());
if (const_cast<BlockLduMatrix<Type>& >(this->matrix_).interfaces().set(patchi))
{
mBouCoeffs.set(patchi, coupleUpperTemp[patchi]);
if (mBouCoeffs[patchi].activeType() == blockCoeffBase::SQUARE)
{
squareTypeField& activeMBouCoeffs = mBouCoeffs[patchi].asSquare();
forAll (activeMBouCoeffs, intI)
{
activeMBouCoeffs[intI] = -activeMBouCoeffs[intI];
}
}
else if (mBouCoeffs[patchi].activeType() == blockCoeffBase::LINEAR)
{
linearTypeField& activeMBouCoeffs = mBouCoeffs[patchi].asLinear();
forAll (activeMBouCoeffs, intI)
{
activeMBouCoeffs[intI] = -activeMBouCoeffs[intI];
}
}
}
}
for (label sweep = 0; sweep < nSweeps_; sweep++)
{
bPrime_ = b;
this->matrix_.initInterfaces
(
mBouCoeffs,
this->matrix_.coupleUpper(),
bPrime_,
x
x,
true // switch to lhs of system
);
this->matrix_.updateInterfaces
(
mBouCoeffs,
this->matrix_.coupleUpper(),
bPrime_,
x
x,
true // switch to lhs of system
);
register label fStart, fEnd, curCoeff;

View file

@ -117,7 +117,8 @@ public:
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{}
@ -129,7 +130,8 @@ public:
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const = 0;
};

View file

@ -937,7 +937,8 @@ public:
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& xif,
scalarField& result,
const direction cmpt
const direction cmpt,
const bool switchToLhs = false
) const;
//- Update coupled interfaces for matrix operations
@ -947,7 +948,8 @@ public:
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& xif,
scalarField& result,
const direction cmpt
const direction cmpt,
const bool switchToLhs = false
) const;
//- Initialise the update of coupled interfaces
@ -958,7 +960,8 @@ public:
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& xif,
scalarField& result,
const direction cmpt
const direction cmpt,
const bool switchToLhs = false
) const;
//- Update coupled interfaces for matrix operations
@ -969,7 +972,8 @@ public:
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& xif,
scalarField& result,
const direction cmpt
const direction cmpt,
const bool switchToLhs = false
) const;

View file

@ -79,7 +79,10 @@ template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::lduMatrix::faceH(const Field<Type>& psi) const
{
tmp<Field<Type> > tfaceHpsi(new Field<Type> (lduAddr().lowerAddr().size()));
tmp<Field<Type> > tfaceHpsi
(
new Field<Type> (lduAddr().lowerAddr().size())
);
Field<Type>& faceHpsi = tfaceHpsi();
if (lowerPtr_ || upperPtr_)

View file

@ -34,7 +34,8 @@ void Foam::lduMatrix::initMatrixInterfaces
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& psiif,
scalarField& result,
const direction cmpt
const direction cmpt,
const bool switchToLhs
) const
{
if
@ -54,7 +55,8 @@ void Foam::lduMatrix::initMatrixInterfaces
*this,
coupleCoeffs[interfaceI],
cmpt,
Pstream::defaultCommsType
Pstream::defaultCommsType,
switchToLhs
);
}
}
@ -81,7 +83,8 @@ void Foam::lduMatrix::initMatrixInterfaces
*this,
coupleCoeffs[interfaceI],
cmpt,
Pstream::blocking
Pstream::blocking,
switchToLhs
);
}
}
@ -102,7 +105,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
const lduInterfaceFieldPtrsList& interfaces,
const scalarField& psiif,
scalarField& result,
const direction cmpt
const direction cmpt,
const bool switchToLhs
) const
{
if
@ -129,7 +133,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
*this,
coupleCoeffs[interfaceI],
cmpt,
Pstream::defaultCommsType
Pstream::defaultCommsType,
switchToLhs
);
}
}
@ -154,7 +159,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
*this,
coupleCoeffs[interfaceI],
cmpt,
Pstream::scheduled
Pstream::scheduled,
switchToLhs
);
}
else
@ -166,7 +172,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
*this,
coupleCoeffs[interfaceI],
cmpt,
Pstream::scheduled
Pstream::scheduled,
switchToLhs
);
}
}
@ -190,7 +197,8 @@ void Foam::lduMatrix::updateMatrixInterfaces
*this,
coupleCoeffs[interfaceI],
cmpt,
Pstream::blocking
Pstream::blocking,
switchToLhs
);
}
}

View file

@ -93,7 +93,7 @@ void Foam::GaussSeidelSmoother::smooth
matrix_.lduAddr().ownerStartAddr().begin();
// Parallel boundary initialisation. The parallel boundary is treated
// Coupled boundary initialisation. The coupled boundary is treated
// as an effective jacobi interface in the boundary.
// Note: there is a change of sign in the coupled
// interface update. The reason for this is that the
@ -105,43 +105,40 @@ void Foam::GaussSeidelSmoother::smooth
// To compensate for this, it is necessary to turn the
// sign of the contribution.
FieldField<Field, scalar> mBouCoeffs(coupleBouCoeffs_.size());
forAll(mBouCoeffs, patchi)
{
if (interfaces_.set(patchi))
{
mBouCoeffs.set(patchi, -coupleBouCoeffs_[patchi]);
}
}
// Handled by LHS switch on initMatrixInterfaces and updateMatrixInterfaces
// HJ, 22/May/2013
for (label sweep = 0; sweep < nSweeps; sweep++)
{
bPrime = b;
// Update from lhs
matrix_.initMatrixInterfaces
(
mBouCoeffs,
coupleBouCoeffs_,
interfaces_,
x,
bPrime,
cmpt
cmpt,
true // switch to lhs
);
// Update from lhs
matrix_.updateMatrixInterfaces
(
mBouCoeffs,
coupleBouCoeffs_,
interfaces_,
x,
bPrime,
cmpt
cmpt,
true // switch to lhs
);
register scalar curX;
register label fStart;
register label fEnd = ownStartPtr[0];
for (register label cellI=0; cellI<nCells; cellI++)
for (register label cellI = 0; cellI < nCells; cellI++)
{
// Start and end of this row
fStart = fEnd;
@ -151,7 +148,7 @@ void Foam::GaussSeidelSmoother::smooth
curX = bPrimePtr[cellI];
// Accumulate the owner product side
for (register label curFace=fStart; curFace<fEnd; curFace++)
for (register label curFace = fStart; curFace < fEnd; curFace++)
{
curX -= upperPtr[curFace]*xPtr[uPtr[curFace]];
}
@ -160,7 +157,7 @@ void Foam::GaussSeidelSmoother::smooth
curX /= diagPtr[cellI];
// Distribute the neighbour side using current x
for (register label curFace=fStart; curFace<fEnd; curFace++)
for (register label curFace = fStart; curFace < fEnd; curFace++)
{
bPrimePtr[uPtr[curFace]] -= lowerPtr[curFace]*curX;
}

View file

@ -78,7 +78,8 @@ void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{
scalarField pnf(size());
@ -95,10 +96,20 @@ void Foam::cyclicGAMGInterfaceField::updateInterfaceMatrix
transformCoupleField(pnf, cmpt);
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}

View file

@ -117,7 +117,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -81,7 +81,8 @@ void Foam::ggiGAMGInterfaceField::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// This must have a reduce in it. HJ, 15/May/2009
@ -96,7 +97,8 @@ void Foam::ggiGAMGInterfaceField::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Get expanded data to zone size. No global reduce allowed
@ -117,21 +119,20 @@ void Foam::ggiGAMGInterfaceField::updateInterfaceMatrix
<< abort(FatalError);
}
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*zonePnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*zonePnf[elemI];
}
// Old treatment
#if(0)
const labelList& za = ggiInterface_.zoneAddressing();
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*zonePnf[za[elemI]];
}
#endif
}

View file

@ -123,7 +123,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -134,7 +135,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -87,7 +87,8 @@ void Foam::mixingPlaneGAMGInterfaceField::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
mixingPlaneInterface_.initInternalFieldTransfer(commsType, psiInternal);
@ -101,7 +102,8 @@ void Foam::mixingPlaneGAMGInterfaceField::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
scalarField pnf =
@ -110,10 +112,20 @@ void Foam::mixingPlaneGAMGInterfaceField::updateInterfaceMatrix
const unallocLabelList& faceCells = mixingPlaneInterface_.faceCells();
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}

View file

@ -126,7 +126,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -137,7 +138,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -78,7 +78,8 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procInterface_.compressedSend
@ -96,7 +97,8 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
scalarField pnf
@ -107,10 +109,20 @@ void Foam::processorGAMGInterfaceField::updateInterfaceMatrix
const unallocLabelList& faceCells = procInterface_.faceCells();
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}

View file

@ -117,7 +117,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -128,7 +129,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -72,7 +72,8 @@ void Foam::regionCoupleGAMGInterfaceField::initInterfaceMatrixUpdate
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// This must have a reduce in it. HJ, 15/May/2009
@ -85,7 +86,8 @@ void Foam::regionCoupleGAMGInterfaceField::initInterfaceMatrixUpdate
m,
coeffs,
cmpt,
commsType
commsType,
switchToLhs
);
}
}
@ -98,7 +100,8 @@ void Foam::regionCoupleGAMGInterfaceField::updateInterfaceMatrix
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Get expanded data to zone size. No global reduce allowed
@ -112,7 +115,8 @@ void Foam::regionCoupleGAMGInterfaceField::updateInterfaceMatrix
m,
coeffs,
cmpt,
commsType
commsType,
switchToLhs
);
}
}

View file

@ -93,7 +93,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -104,7 +105,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
};

View file

@ -45,7 +45,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class coupledFaPatch Declaration
Class coupledFaPatchField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
@ -183,7 +183,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const = 0;
//- Write

View file

@ -185,7 +185,8 @@ void cyclicFaPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{
scalarField pnf(this->size());
@ -203,10 +204,20 @@ void cyclicFaPatchField<Type>::updateInterfaceMatrix
transformCoupleField(pnf, cmpt);
// Multiply the field by coefficients and add into the result
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}

View file

@ -163,7 +163,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -215,7 +215,8 @@ void processorFaPatchField<Type>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procPatch_.compressedSend
@ -234,7 +235,8 @@ void processorFaPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
scalarField pnf
@ -249,10 +251,20 @@ void processorFaPatchField<Type>::updateInterfaceMatrix
const unallocLabelList& edgeFaces = this->patch().edgeFaces();
if (switchToLhs)
{
forAll(edgeFaces, elemI)
{
result[edgeFaces[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(edgeFaces, elemI)
{
result[edgeFaces[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}

View file

@ -178,7 +178,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -189,7 +190,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Processor coupled interface functions

View file

@ -41,7 +41,8 @@ void processorFaPatchField<scalar>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procPatch_.compressedSend
@ -60,7 +61,8 @@ void processorFaPatchField<scalar>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
scalarField pnf
@ -70,10 +72,20 @@ void processorFaPatchField<scalar>::updateInterfaceMatrix
const unallocLabelList& edgeFaces = patch().edgeFaces();
if (switchToLhs)
{
forAll(edgeFaces, facei)
{
result[edgeFaces[facei]] -= coeffs[facei]*pnf[facei];
}
}
else
{
forAll(edgeFaces, facei)
{
result[edgeFaces[facei]] -= coeffs[facei]*pnf[facei];
}
}
}

View file

@ -44,7 +44,8 @@ void processorFaPatchField<scalar>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhsswitchToLhs
) const;
@ -56,7 +57,8 @@ void processorFaPatchField<scalar>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -24,8 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "basicSymmetryFvPatchFields.H"
#include "fvPatchFields.H"
#include "basicSymmetryFvPatchFields.H"
#include "volMesh.H"
#include "addToRunTimeSelectionTable.H"

View file

@ -157,6 +157,7 @@ tmp<Field<Type> > coupledFvPatchField<Type>::valueInternalCoeffs
return Type(pTraits<Type>::one)*w;
}
template<class Type>
tmp<Field<Type> > coupledFvPatchField<Type>::valueBoundaryCoeffs
(

View file

@ -38,6 +38,8 @@ SourceFiles
#include "BlockLduInterfaceField.H"
#include "CoeffField.H"
#include "lduInterfaceField.H"
#include "fvPatchField.H"
#include "coupledFvPatch.H"
@ -47,7 +49,7 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class coupledFvPatch Declaration
Class coupledFvPatchField Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
@ -186,9 +188,11 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const = 0;
// Block coupled interface functionality
//- Initialise neighbour matrix update
@ -198,7 +202,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{}
@ -209,10 +214,14 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
notImplemented("coupledFvPatchField<Type>::updateInterfaceMatrix for block matrices")
notImplemented
(
"coupledFvPatchField<Type>::updateInterfaceMatrix"
);
}
//- Write

View file

@ -185,7 +185,8 @@ void cyclicFvPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{
scalarField pnf(this->size());
@ -203,10 +204,20 @@ void cyclicFvPatchField<Type>::updateInterfaceMatrix
transformCoupleField(pnf, cmpt);
// Multiply the field by coefficients and add into the result
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}

View file

@ -162,7 +162,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
@ -175,7 +176,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
notImplemented

View file

@ -228,7 +228,8 @@ void cyclicGgiFvPatchField<Type>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Communication is allowed either before or after processor
@ -256,10 +257,20 @@ void cyclicGgiFvPatchField<Type>::initInterfaceMatrixUpdate
// Multiply the field by coefficients and add into the result
const unallocLabelList& fc = cyclicGgiPatch_.faceCells();
if (switchToLhs)
{
forAll(fc, elemI)
{
result[fc[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(fc, elemI)
{
result[fc[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}
@ -272,7 +283,8 @@ void cyclicGgiFvPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{}

View file

@ -152,7 +152,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -163,7 +164,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -246,7 +246,8 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Communication is allowed either before or after processor
@ -267,10 +268,20 @@ void ggiFvPatchField<Type>::initInterfaceMatrixUpdate
// Multiply the field by coefficients and add into the result
const unallocLabelList& fc = ggiPatch_.faceCells();
if (switchToLhs)
{
forAll(fc, elemI)
{
result[fc[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(fc, elemI)
{
result[fc[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}
@ -282,7 +293,8 @@ void ggiFvPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{}

View file

@ -155,7 +155,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -166,7 +167,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -146,7 +146,8 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{
scalarField pnf(this->size());
@ -178,10 +179,20 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
this->transformCoupleField(pnf, cmpt);
// Multiply the field by coefficients and add into the result
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}

View file

@ -126,7 +126,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
};

View file

@ -462,7 +462,8 @@ void mixingPlaneFvPatchField<Type>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Communication is allowed either before or after processor
@ -484,10 +485,20 @@ void mixingPlaneFvPatchField<Type>::initInterfaceMatrixUpdate
scalarField pnf = mixingPlanePatch_.interpolate(sField);
// Multiply the field by coefficients and add into the result
if (switchToLhs)
{
forAll(fc, elemI)
{
result[fc[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(fc, elemI)
{
result[fc[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}
@ -499,7 +510,8 @@ void mixingPlaneFvPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{}

View file

@ -234,7 +234,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -245,7 +246,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -194,7 +194,8 @@ void overlapGgiFvPatchField<Type>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Communication is allowed either before or after processor
@ -215,10 +216,20 @@ void overlapGgiFvPatchField<Type>::initInterfaceMatrixUpdate
// Multiply the field by coefficients and add into the result
const unallocLabelList& fc = overlapGgiPatch_.faceCells();
if (switchToLhs)
{
forAll(fc, elemI)
{
result[fc[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(fc, elemI)
{
result[fc[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}
@ -231,7 +242,8 @@ void overlapGgiFvPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{}

View file

@ -152,7 +152,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -163,7 +164,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -218,7 +218,8 @@ void processorFvPatchField<Type>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procPatch_.compressedSend
@ -237,7 +238,8 @@ void processorFvPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
scalarField pnf
@ -249,13 +251,22 @@ void processorFvPatchField<Type>::updateInterfaceMatrix
transformCoupleField(pnf, cmpt);
// Multiply the field by coefficients and add into the result
const unallocLabelList& faceCells = this->patch().faceCells();
if (switchToLhs)
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(faceCells, elemI)
{
result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}

View file

@ -180,7 +180,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -191,9 +192,11 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
// Block coupled interface functionality
//- Initialise neighbour matrix update
@ -203,7 +206,8 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{}
@ -214,12 +218,18 @@ public:
Field<Type>&,
const BlockLduMatrix<Type>&,
const CoeffField<Type>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
notImplemented("processorFvPatchField<Type>::updateInterfaceMatrix for block matrices")
notImplemented
(
"processorFvPatchField<Type>::updateInterfaceMatrix "
"for block matrices"
);
}
//- Processor coupled interface functions
//- Return processor number

View file

@ -41,7 +41,8 @@ void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procPatch_.compressedSend
@ -60,7 +61,8 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
scalarField pnf
@ -70,10 +72,20 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
const unallocLabelList& faceCells = patch().faceCells();
forAll(faceCells, facei)
if (switchToLhs)
{
forAll (faceCells, facei)
{
result[faceCells[facei]] += coeffs[facei]*pnf[facei];
}
}
else
{
forAll (faceCells, facei)
{
result[faceCells[facei]] -= coeffs[facei]*pnf[facei];
}
}
}
@ -84,7 +96,8 @@ void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
scalarField&,
const BlockLduMatrix<scalar>&,
const CoeffField<scalar>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
procPatch_.compressedSend
@ -102,7 +115,8 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
scalarField& result,
const BlockLduMatrix<scalar>&,
const CoeffField<scalar>& coeffs,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
scalarField pnf
@ -113,12 +127,23 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
const unallocLabelList& faceCells = patch().faceCells();
const scalarField& scalarCoeffs = coeffs.asScalar();
forAll(faceCells, facei)
if (switchToLhs)
{
forAll (faceCells, facei)
{
result[faceCells[facei]] += scalarCoeffs[facei]*pnf[facei];
}
}
else
{
forAll (faceCells, facei)
{
result[faceCells[facei]] -= scalarCoeffs[facei]*pnf[facei];
}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View file

@ -45,7 +45,8 @@ void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
@ -57,7 +58,8 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
@ -68,7 +70,8 @@ void processorFvPatchField<scalar>::initInterfaceMatrixUpdate
scalarField&,
const BlockLduMatrix<scalar>&,
const CoeffField<scalar>&,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
@ -79,7 +82,8 @@ void processorFvPatchField<scalar>::updateInterfaceMatrix
scalarField& result,
const BlockLduMatrix<scalar>&,
const CoeffField<scalar>& coeffs,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -456,7 +456,8 @@ void regionCouplingFvPatchField<Type>::initInterfaceMatrixUpdate
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{
if (regionCouplePatch_.coupled())
@ -494,7 +495,8 @@ void regionCouplingFvPatchField<Type>::updateInterfaceMatrix
const lduMatrix&,
const scalarField& coeffs,
const direction ,
const Pstream::commsTypes
const Pstream::commsTypes,
const bool switchToLhs
) const
{
if (regionCouplePatch_.coupled())
@ -508,11 +510,21 @@ void regionCouplingFvPatchField<Type>::updateInterfaceMatrix
// Multiply the field by coefficients and add into the result
const unallocLabelList& fc = regionCouplePatch_.faceCells();
if (switchToLhs)
{
forAll(fc, elemI)
{
result[fc[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll(fc, elemI)
{
result[fc[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}
else
{
FatalErrorIn

View file

@ -232,7 +232,8 @@ public:
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
@ -243,7 +244,8 @@ public:
const lduMatrix&,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;

View file

@ -26,7 +26,7 @@ License
#include "fixedInternalValueFvPatchField.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
#include "fvMatrices.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //

View file

@ -24,8 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "fixedInternalValueFvPatchFields.H"
#include "fvPatchFields.H"
#include "fixedInternalValueFvPatchFields.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"

View file

@ -24,8 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "volFields.H"
#include "surfaceFields.H"
#include "fvPatchFields.H"
#include "freestreamFvPatchFields.H"
#include "addToRunTimeSelectionTable.H"

View file

@ -24,11 +24,9 @@ License
\*---------------------------------------------------------------------------*/
#include "fvPatchFields.H"
#include "syringePressureFvPatchScalarField.H"
#include "volMesh.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -184,15 +184,7 @@ tmp<scalarField> waveTransmissiveFvPatchField<Type>::supercritical() const
template<class Type>
void waveTransmissiveFvPatchField<Type>::write(Ostream& os) const
{
fvPatchField<Type>::write(os);
if (this->phiName_ != "phi")
{
os.writeKeyword("phi") << this->phiName_ << token::END_STATEMENT << nl;
}
if (this->rhoName_ != "rho")
{
os.writeKeyword("rho") << this->rhoName_ << token::END_STATEMENT << nl;
}
advectiveFvPatchField<Type>::write(os);
if (this->UName_ != "U")
{
os.writeKeyword("U") << this->UName_ << token::END_STATEMENT << nl;
@ -203,14 +195,6 @@ void waveTransmissiveFvPatchField<Type>::write(Ostream& os) const
}
os.writeKeyword("gamma") << gamma_ << token::END_STATEMENT << nl;
if (this->lInf_ > SMALL)
{
os.writeKeyword("fieldInf") << this->fieldInf_
<< token::END_STATEMENT << nl;
os.writeKeyword("lInf") << this->lInf_
<< token::END_STATEMENT << nl;
}
}

View file

@ -164,7 +164,8 @@ const Foam::objectRegistry& Foam::fvPatchField<Type>::db() const
template<class Type>
template<class GeometricField, class Type2>
const typename GeometricField::PatchFieldType& Foam::fvPatchField<Type>::lookupPatchField
const typename GeometricField::PatchFieldType&
Foam::fvPatchField<Type>::lookupPatchField
(
const word& name,
const GeometricField*,
@ -247,6 +248,47 @@ void Foam::fvPatchField<Type>::manipulateMatrix(fvMatrix<Type>& matrix)
}
template<class Type>
void Foam::fvPatchField<Type>::patchFlux
(
GeometricField<Type, fvsPatchField, surfaceMesh>& pFlux,
const fvMatrix<Type>& matrix
) const
{
// Code moved from fvMatrix.C
// HJ, 29/May/2013
const label patchI = this->patch().index();
// Virtual function for patch flux. HJ, 29/May/2013
if (this->coupled())
{
// Coupled patch
pFlux.boundaryField()[patchI] =
cmptMultiply
(
matrix.internalCoeffs()[patchI],
this->patchInternalField()
)
- cmptMultiply
(
matrix.boundaryCoeffs()[patchI],
this->patchNeighbourField()
);
}
else
{
// Uncoupled patch
pFlux.boundaryField()[patchI] =
cmptMultiply
(
matrix.internalCoeffs()[patchI],
this->patchInternalField()
)
- matrix.boundaryCoeffs()[patchI];
}
}
template<class Type>
void Foam::fvPatchField<Type>::write(Ostream& os) const
{

View file

@ -59,13 +59,19 @@ class objectRegistry;
class dictionary;
class fvPatchFieldMapper;
class volMesh;
class surfaceMesh;
// Forward declaration of friend functions and operators
template<class Type>
class fvPatchField;
template<class Type>
class fvsPatchField;
template<class Type, template<class> class PatchField, class GeoMesh>
class GeometricField;
template<class Type>
class fvMatrix;
@ -433,6 +439,13 @@ public:
//- Manipulate matrix
virtual void manipulateMatrix(fvMatrix<Type>& matrix);
//- Calculate patch flux
virtual void patchFlux
(
GeometricField<Type, fvsPatchField, surfaceMesh>& pFlux,
const fvMatrix<Type>& matrix
) const;
// I-O

View file

@ -30,6 +30,10 @@ License
#include "fvPatchField.H"
#include "fvPatchFieldsFwd.H"
// Added for instantiation of patchFlux templates. HJ, 29/May/2013
#include "fvMatrices.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif

View file

@ -41,6 +41,7 @@ SourceFiles
#include "fvPatchField.H"
#include "volFieldsFwd.H"
#include "calculatedFvPatchFields.H"
#include "fvMatrices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

View file

@ -877,39 +877,18 @@ flux() const
);
}
FieldField<Field, Type> InternalContrib = internalCoeffs_;
forAll(InternalContrib, patchI)
// This needs to go into virtual functions for all coupled patches
// in order to simplify handling of overset meshes
// HJ, 29/May/2013
forAll (psi_.boundaryField(), patchI)
{
InternalContrib[patchI] =
cmptMultiply
psi_.boundaryField()[patchI].patchFlux
(
InternalContrib[patchI],
psi_.boundaryField()[patchI].patchInternalField()
fieldFlux,
*this
);
}
FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_;
forAll(NeighbourContrib, patchI)
{
if (psi_.boundaryField()[patchI].coupled())
{
NeighbourContrib[patchI] =
cmptMultiply
(
NeighbourContrib[patchI],
psi_.boundaryField()[patchI].patchNeighbourField()
);
}
}
forAll(fieldFlux.boundaryField(), patchI)
{
fieldFlux.boundaryField()[patchI] =
InternalContrib[patchI] - NeighbourContrib[patchI];
}
if (faceFluxCorrectionPtr_)
{
fieldFlux += *faceFluxCorrectionPtr_;

View file

@ -195,6 +195,10 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Type>::solve()
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::fvMatrix<Type>::residual() const
{
// Complete matrix assembly. HJ, 17/Apr/2012
fvMatrix& m = const_cast<fvMatrix&>(*this);
m.completeAssembly();
// Bug fix: Creating a tmp out of a const reference will change the field
// HJ, 15/Apr/2011
tmp<Field<Type> > tres(new Field<Type>(source_));

View file

@ -167,6 +167,7 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
solverPerf.print();
// Diagonal has been restored, clear complete assembly flag?
diag() = saveDiag;
psi_.correctBoundaryConditions();
@ -178,6 +179,10 @@ Foam::lduMatrix::solverPerformance Foam::fvMatrix<Foam::scalar>::solve
template<>
Foam::tmp<Foam::scalarField> Foam::fvMatrix<Foam::scalar>::residual() const
{
// Complete matrix assembly. HJ, 17/Apr/2012
fvMatrix<scalar>& m = const_cast<fvMatrix<scalar>&>(*this);
m.completeAssembly();
scalarField boundaryDiag(psi_.size(), 0.0);
addBoundaryDiag(boundaryDiag, 0);