mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-12-28 03:37:59 +00:00
2972 lines
71 KiB
C
2972 lines
71 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2011-2017 OpenFOAM Foundation
|
|
Copyright (C) 2016-2023 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "volFields.H"
|
|
#include "surfaceFields.H"
|
|
#include "calculatedFvPatchFields.H"
|
|
#include "extrapolatedCalculatedFvPatchFields.H"
|
|
#include "coupledFvPatchFields.H"
|
|
#include "IndirectList.H"
|
|
#include "UniformList.H"
|
|
#include "demandDrivenData.H"
|
|
|
|
#include "cyclicFvPatchField.H"
|
|
#include "cyclicAMIFvPatchField.H"
|
|
#include "cyclicACMIFvPatchField.H"
|
|
|
|
#include "processorLduInterfaceField.H"
|
|
|
|
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
|
|
|
|
template<class Type>
|
|
template<class Type2>
|
|
void Foam::fvMatrix<Type>::addToInternalField
|
|
(
|
|
const labelUList& addr,
|
|
const Field<Type2>& pf,
|
|
Field<Type2>& intf
|
|
) const
|
|
{
|
|
if (addr.size() != pf.size())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "addressing (" << addr.size()
|
|
<< ") and field (" << pf.size() << ") are different sizes" << endl
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
forAll(addr, facei)
|
|
{
|
|
intf[addr[facei]] += pf[facei];
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
template<class Type2>
|
|
void Foam::fvMatrix<Type>::addToInternalField
|
|
(
|
|
const labelUList& addr,
|
|
const tmp<Field<Type2>>& tpf,
|
|
Field<Type2>& intf
|
|
) const
|
|
{
|
|
addToInternalField(addr, tpf(), intf);
|
|
tpf.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
template<class Type2>
|
|
void Foam::fvMatrix<Type>::subtractFromInternalField
|
|
(
|
|
const labelUList& addr,
|
|
const Field<Type2>& pf,
|
|
Field<Type2>& intf
|
|
) const
|
|
{
|
|
if (addr.size() != pf.size())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "addressing (" << addr.size()
|
|
<< ") and field (" << pf.size() << ") are different sizes" << endl
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
forAll(addr, facei)
|
|
{
|
|
intf[addr[facei]] -= pf[facei];
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
template<class Type2>
|
|
void Foam::fvMatrix<Type>::subtractFromInternalField
|
|
(
|
|
const labelUList& addr,
|
|
const tmp<Field<Type2>>& tpf,
|
|
Field<Type2>& intf
|
|
) const
|
|
{
|
|
subtractFromInternalField(addr, tpf(), intf);
|
|
tpf.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::addBoundaryDiag
|
|
(
|
|
scalarField& diag,
|
|
const direction solveCmpt
|
|
) const
|
|
{
|
|
for (label fieldi = 0; fieldi < nMatrices(); ++fieldi)
|
|
{
|
|
const auto& bpsi = this->psi(fieldi).boundaryField();
|
|
|
|
forAll(bpsi, ptfi)
|
|
{
|
|
const label patchi = globalPatchID(fieldi, ptfi);
|
|
|
|
if (patchi != -1)
|
|
{
|
|
addToInternalField
|
|
(
|
|
lduAddr().patchAddr(patchi),
|
|
internalCoeffs_[patchi].component(solveCmpt),
|
|
diag
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::addCmptAvBoundaryDiag(scalarField& diag) const
|
|
{
|
|
for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
|
|
{
|
|
const auto& bpsi = this->psi(fieldi).boundaryField();
|
|
|
|
forAll(bpsi, ptfi)
|
|
{
|
|
const label patchi = globalPatchID(fieldi, ptfi);
|
|
if (patchi != -1)
|
|
{
|
|
addToInternalField
|
|
(
|
|
lduAddr().patchAddr(patchi),
|
|
cmptAv(internalCoeffs_[patchi]),
|
|
diag
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::addBoundarySource
|
|
(
|
|
Field<Type>& source,
|
|
const bool couples
|
|
) const
|
|
{
|
|
for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
|
|
{
|
|
const auto& bpsi = this->psi(fieldi).boundaryField();
|
|
|
|
forAll(bpsi, ptfi)
|
|
{
|
|
const fvPatchField<Type>& ptf = bpsi[ptfi];
|
|
|
|
const label patchi = globalPatchID(fieldi, ptfi);
|
|
|
|
if (patchi != -1)
|
|
{
|
|
const Field<Type>& pbc = boundaryCoeffs_[patchi];
|
|
|
|
if (!ptf.coupled())
|
|
{
|
|
addToInternalField
|
|
(
|
|
lduAddr().patchAddr(patchi),
|
|
pbc,
|
|
source
|
|
);
|
|
}
|
|
else if (couples)
|
|
{
|
|
const tmp<Field<Type>> tpnf = ptf.patchNeighbourField();
|
|
const Field<Type>& pnf = tpnf();
|
|
|
|
const labelUList& addr = lduAddr().patchAddr(patchi);
|
|
|
|
forAll(addr, facei)
|
|
{
|
|
source[addr[facei]] +=
|
|
cmptMultiply(pbc[facei], pnf[facei]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
template<template<class> class ListType>
|
|
void Foam::fvMatrix<Type>::setValuesFromList
|
|
(
|
|
const labelUList& cellLabels,
|
|
const ListType<Type>& values
|
|
)
|
|
{
|
|
const fvMesh& mesh = psi_.mesh();
|
|
|
|
const cellList& cells = mesh.cells();
|
|
const labelUList& own = mesh.owner();
|
|
const labelUList& nei = mesh.neighbour();
|
|
|
|
scalarField& Diag = diag();
|
|
Field<Type>& psi =
|
|
const_cast
|
|
<
|
|
GeometricField<Type, fvPatchField, volMesh>&
|
|
>(psi_).primitiveFieldRef();
|
|
|
|
|
|
// Following actions:
|
|
// - adjust local field psi
|
|
// - set local matrix to be diagonal (so adjust source)
|
|
// - cut connections to neighbours
|
|
// - make (on non-adjusted cells) contribution explicit
|
|
|
|
if (symmetric() || asymmetric())
|
|
{
|
|
forAll(cellLabels, i)
|
|
{
|
|
const label celli = cellLabels[i];
|
|
const Type& value = values[i];
|
|
|
|
for (const label facei : cells[celli])
|
|
{
|
|
if (mesh.isInternalFace(facei))
|
|
{
|
|
if (symmetric())
|
|
{
|
|
if (celli == own[facei])
|
|
{
|
|
source_[nei[facei]] -= upper()[facei]*value;
|
|
}
|
|
else
|
|
{
|
|
source_[own[facei]] -= upper()[facei]*value;
|
|
}
|
|
|
|
upper()[facei] = 0.0;
|
|
}
|
|
else
|
|
{
|
|
if (celli == own[facei])
|
|
{
|
|
source_[nei[facei]] -= lower()[facei]*value;
|
|
}
|
|
else
|
|
{
|
|
source_[own[facei]] -= upper()[facei]*value;
|
|
}
|
|
|
|
upper()[facei] = 0.0;
|
|
lower()[facei] = 0.0;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
const label patchi = mesh.boundaryMesh().whichPatch(facei);
|
|
|
|
if (internalCoeffs_[patchi].size())
|
|
{
|
|
const label patchFacei =
|
|
mesh.boundaryMesh()[patchi].whichFace(facei);
|
|
|
|
internalCoeffs_[patchi][patchFacei] = Zero;
|
|
boundaryCoeffs_[patchi][patchFacei] = Zero;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Note: above loop might have affected source terms on adjusted cells
|
|
// so make sure to adjust them afterwards
|
|
forAll(cellLabels, i)
|
|
{
|
|
const label celli = cellLabels[i];
|
|
const Type& value = values[i];
|
|
|
|
psi[celli] = value;
|
|
source_[celli] = value*Diag[celli];
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
bool Foam::fvMatrix<Type>::checkImplicit(const label fieldi)
|
|
{
|
|
const auto& bpsi = this->psi(fieldi).boundaryField();
|
|
|
|
word idName;
|
|
forAll(bpsi, patchi)
|
|
{
|
|
if (bpsi[patchi].useImplicit())
|
|
{
|
|
if (debug)
|
|
{
|
|
Pout<< "fvMatrix<Type>::checkImplicit "
|
|
<< " field:" << this->psi(fieldi).name()
|
|
<< " on mesh:"
|
|
<< this->psi(fieldi).mesh().name()
|
|
<< " patch:" << bpsi[patchi].patch().name()
|
|
<< endl;
|
|
}
|
|
|
|
idName += Foam::name(patchi);
|
|
useImplicit_ = true;
|
|
}
|
|
}
|
|
|
|
if (useImplicit_)
|
|
{
|
|
lduAssemblyName_ = word("lduAssembly") + idName;
|
|
}
|
|
|
|
return !idName.empty();
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
template<class Type>
|
|
Foam::fvMatrix<Type>::fvMatrix
|
|
(
|
|
const GeometricField<Type, fvPatchField, volMesh>& psi,
|
|
const dimensionSet& ds
|
|
)
|
|
:
|
|
lduMatrix(psi.mesh()),
|
|
psi_(psi),
|
|
useImplicit_(false),
|
|
lduAssemblyName_(),
|
|
nMatrix_(0),
|
|
dimensions_(ds),
|
|
source_(psi.size(), Zero),
|
|
internalCoeffs_(psi.mesh().boundary().size()),
|
|
boundaryCoeffs_(psi.mesh().boundary().size()),
|
|
faceFluxCorrectionPtr_(nullptr)
|
|
{
|
|
DebugInFunction
|
|
<< "Constructing fvMatrix<Type> for field " << psi_.name() << endl;
|
|
|
|
checkImplicit();
|
|
|
|
forAll(psi.mesh().boundary(), patchi)
|
|
{
|
|
internalCoeffs_.set
|
|
(
|
|
patchi,
|
|
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
|
|
);
|
|
|
|
boundaryCoeffs_.set
|
|
(
|
|
patchi,
|
|
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
|
|
);
|
|
}
|
|
|
|
auto& psiRef = this->psi(0);
|
|
const label currentStatePsi = psiRef.eventNo();
|
|
psiRef.boundaryFieldRef().updateCoeffs();
|
|
psiRef.eventNo() = currentStatePsi;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::fvMatrix<Type>::fvMatrix(const fvMatrix<Type>& fvm)
|
|
:
|
|
lduMatrix(fvm),
|
|
psi_(fvm.psi_),
|
|
useImplicit_(fvm.useImplicit_),
|
|
lduAssemblyName_(fvm.lduAssemblyName_),
|
|
nMatrix_(fvm.nMatrix_),
|
|
dimensions_(fvm.dimensions_),
|
|
source_(fvm.source_),
|
|
internalCoeffs_(fvm.internalCoeffs_),
|
|
boundaryCoeffs_(fvm.boundaryCoeffs_),
|
|
faceFluxCorrectionPtr_(nullptr)
|
|
{
|
|
DebugInFunction
|
|
<< "Copying fvMatrix<Type> for field " << psi_.name() << endl;
|
|
|
|
if (fvm.faceFluxCorrectionPtr_)
|
|
{
|
|
faceFluxCorrectionPtr_ =
|
|
new GeometricField<Type, fvsPatchField, surfaceMesh>
|
|
(
|
|
*(fvm.faceFluxCorrectionPtr_)
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::fvMatrix<Type>::fvMatrix(const tmp<fvMatrix<Type>>& tmat)
|
|
:
|
|
lduMatrix(tmat.constCast(), tmat.movable()),
|
|
psi_(tmat().psi_),
|
|
useImplicit_(tmat().useImplicit_),
|
|
lduAssemblyName_(tmat().lduAssemblyName_),
|
|
nMatrix_(tmat().nMatrix_),
|
|
dimensions_(tmat().dimensions_),
|
|
source_(tmat.constCast().source_, tmat.movable()),
|
|
internalCoeffs_(tmat.constCast().internalCoeffs_, tmat.movable()),
|
|
boundaryCoeffs_(tmat.constCast().boundaryCoeffs_, tmat.movable()),
|
|
faceFluxCorrectionPtr_(nullptr)
|
|
{
|
|
DebugInFunction
|
|
<< "Copy/move fvMatrix<Type> for field " << psi_.name() << endl;
|
|
|
|
if (tmat().faceFluxCorrectionPtr_)
|
|
{
|
|
if (tmat.movable())
|
|
{
|
|
faceFluxCorrectionPtr_ = tmat().faceFluxCorrectionPtr_;
|
|
tmat().faceFluxCorrectionPtr_ = nullptr;
|
|
}
|
|
else
|
|
{
|
|
faceFluxCorrectionPtr_ =
|
|
new GeometricField<Type, fvsPatchField, surfaceMesh>
|
|
(
|
|
*(tmat().faceFluxCorrectionPtr_)
|
|
);
|
|
}
|
|
}
|
|
|
|
tmat.clear();
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * * //
|
|
|
|
template<class Type>
|
|
Foam::fvMatrix<Type>::~fvMatrix()
|
|
{
|
|
DebugInFunction
|
|
<< "Destroying fvMatrix<Type> for field " << psi_.name() << endl;
|
|
|
|
deleteDemandDrivenData(faceFluxCorrectionPtr_);
|
|
subMatrices_.clear();
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::setInterfaces
|
|
(
|
|
lduInterfaceFieldPtrsList& interfaces,
|
|
PtrDynList<lduInterfaceField>& newInterfaces
|
|
)
|
|
{
|
|
interfaces.setSize(internalCoeffs_.size());
|
|
for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
|
|
{
|
|
const auto& bpsi = this->psi(fieldi).boundaryField();
|
|
lduInterfaceFieldPtrsList fieldInterfaces(bpsi.scalarInterfaces());
|
|
|
|
forAll (fieldInterfaces, patchi)
|
|
{
|
|
label globalPatchID = lduMeshPtr()->patchMap()[fieldi][patchi];
|
|
|
|
if (globalPatchID != -1)
|
|
{
|
|
if (fieldInterfaces.set(patchi))
|
|
{
|
|
if (isA<cyclicLduInterfaceField>(bpsi[patchi]))
|
|
{
|
|
newInterfaces.append
|
|
(
|
|
new cyclicFvPatchField<Type>
|
|
(
|
|
refCast<const fvPatch>
|
|
(
|
|
lduMeshPtr()->interfaces()[globalPatchID]
|
|
),
|
|
bpsi[patchi].internalField()
|
|
)
|
|
);
|
|
interfaces.set(globalPatchID, &newInterfaces.last());
|
|
|
|
}
|
|
else if (isA<cyclicAMILduInterfaceField>(bpsi[patchi]))
|
|
{
|
|
newInterfaces.append
|
|
(
|
|
new cyclicAMIFvPatchField<Type>
|
|
(
|
|
refCast<const fvPatch>
|
|
(
|
|
lduMeshPtr()->interfaces()[globalPatchID]
|
|
),
|
|
bpsi[patchi].internalField()
|
|
)
|
|
);
|
|
interfaces.set(globalPatchID, &newInterfaces.last());
|
|
}
|
|
else if (isA<cyclicACMILduInterfaceField>(bpsi[patchi]))
|
|
{
|
|
newInterfaces.append
|
|
(
|
|
new cyclicACMIFvPatchField<Type>
|
|
(
|
|
refCast<const fvPatch>
|
|
(
|
|
lduMeshPtr()->interfaces()[globalPatchID]
|
|
),
|
|
bpsi[patchi].internalField()
|
|
)
|
|
);
|
|
interfaces.set(globalPatchID, &newInterfaces.last());
|
|
}
|
|
else
|
|
{
|
|
interfaces.set(globalPatchID, &fieldInterfaces[patchi]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::mapContributions
|
|
(
|
|
label fieldi,
|
|
const FieldField<Field, Type>& fluxContrib,
|
|
FieldField<Field, Type>& contrib,
|
|
bool internal
|
|
) const
|
|
{
|
|
const lduPrimitiveMeshAssembly* ptr = lduMeshPtr();
|
|
|
|
const labelList& patchMap = ptr->patchMap()[fieldi];
|
|
|
|
forAll(contrib, patchi)
|
|
{
|
|
const label globalPtchId = patchMap[patchi];
|
|
|
|
if (globalPtchId != -1)
|
|
{
|
|
// Cache contrib before overwriting
|
|
const Field<Type> saveContrib(fluxContrib[globalPtchId]);
|
|
contrib[patchi].setSize(psi_.boundaryField()[patchi].size()),
|
|
contrib[patchi] = pTraits<Type>::zero;
|
|
|
|
if (internal)
|
|
{
|
|
contrib[patchi] =
|
|
cmptMultiply
|
|
(
|
|
saveContrib,
|
|
psi_.boundaryField()[patchi].patchInternalField()
|
|
);
|
|
}
|
|
else
|
|
{
|
|
if (this->psi(fieldi).boundaryField()[patchi].coupled())
|
|
{
|
|
contrib[patchi] =
|
|
cmptMultiply
|
|
(
|
|
saveContrib,
|
|
psi_.boundaryField()[patchi].patchNeighbourField()
|
|
);
|
|
}
|
|
}
|
|
}
|
|
else if (globalPtchId == -1)
|
|
{
|
|
const polyPatch& pp =
|
|
this->psi(fieldi).mesh().boundaryMesh()[patchi];
|
|
|
|
if (pp.masterImplicit())
|
|
{
|
|
label virtualPatch =
|
|
ptr->patchLocalToGlobalMap()[fieldi][patchi];
|
|
|
|
const label nbrPatchId = pp.neighbPolyPatchID();
|
|
|
|
// Copy contrib before overwriting
|
|
const Field<Type> saveContrib(fluxContrib[virtualPatch]);
|
|
|
|
Field<Type>& coeffs = contrib[patchi];
|
|
Field<Type>& nbrCoeffs = contrib[nbrPatchId];
|
|
|
|
coeffs.setSize(psi_.boundaryField()[patchi].size());
|
|
nbrCoeffs.setSize(psi_.boundaryField()[nbrPatchId].size());
|
|
|
|
coeffs = pTraits<Type>::zero;
|
|
nbrCoeffs = pTraits<Type>::zero;
|
|
|
|
// nrb cells
|
|
const labelList& nbrCellIds =
|
|
ptr->cellBoundMap()[fieldi][patchi];
|
|
|
|
const labelList& cellIds =
|
|
ptr->cellBoundMap()[fieldi][nbrPatchId];
|
|
|
|
const GeometricField<Type, fvPatchField, volMesh>& psi =
|
|
this->psi(fieldi);
|
|
|
|
forAll(saveContrib, subFaceI)
|
|
{
|
|
const label faceId =
|
|
ptr->facePatchFaceMap()[fieldi][patchi][subFaceI];
|
|
const label nbrFaceId =
|
|
ptr->facePatchFaceMap()[fieldi][nbrPatchId][subFaceI];
|
|
|
|
const label nbrCellId = nbrCellIds[subFaceI];
|
|
const label cellId = cellIds[subFaceI];
|
|
|
|
if (internal)
|
|
{
|
|
coeffs[faceId] +=
|
|
cmptMultiply(saveContrib[subFaceI], psi[cellId]);
|
|
|
|
nbrCoeffs[nbrFaceId] +=
|
|
cmptMultiply(saveContrib[subFaceI], psi[nbrCellId]);
|
|
}
|
|
else //boundary
|
|
{
|
|
coeffs[faceId] +=
|
|
cmptMultiply(saveContrib[subFaceI], psi[nbrCellId]);
|
|
|
|
nbrCoeffs[nbrFaceId] +=
|
|
cmptMultiply(saveContrib[subFaceI], psi[cellId]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::setBounAndInterCoeffs()
|
|
{
|
|
// If it is a multi-fvMatrix needs correct internalCoeffs and
|
|
// boundaryCoeffs size
|
|
if (nMatrix_ > 0)
|
|
{
|
|
label interfaceI(0);
|
|
for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
|
|
{
|
|
const auto& psi = this->psi(fieldi);
|
|
|
|
forAll(psi.mesh().boundary(), patchi)
|
|
{
|
|
interfaceI++;
|
|
}
|
|
}
|
|
internalCoeffs_.setSize(interfaceI);
|
|
boundaryCoeffs_.setSize(interfaceI);
|
|
|
|
interfaceI = 0;
|
|
for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
|
|
{
|
|
const auto& psi = this->psi(fieldi);
|
|
|
|
forAll(psi.mesh().boundary(), patchi)
|
|
{
|
|
internalCoeffs_.set
|
|
(
|
|
interfaceI,
|
|
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
|
|
);
|
|
|
|
boundaryCoeffs_.set
|
|
(
|
|
interfaceI,
|
|
new Field<Type>(psi.mesh().boundary()[patchi].size(), Zero)
|
|
);
|
|
interfaceI++;
|
|
}
|
|
}
|
|
}
|
|
|
|
for (label i=0; i < nMatrices(); ++i)
|
|
{
|
|
const auto& bpsi = this->psi(i).boundaryField();
|
|
|
|
// Cache to-be implicit internal/boundary
|
|
FieldField<Field, Type> boundary(bpsi.size());
|
|
FieldField<Field, Type> internal(bpsi.size());
|
|
|
|
label implicit = 0;
|
|
forAll(bpsi, patchI)
|
|
{
|
|
label globalPatchId = lduMeshPtr()->patchMap()[i][patchI];
|
|
if (globalPatchId == -1)
|
|
{
|
|
boundary.set
|
|
(
|
|
implicit,
|
|
matrix(i).boundaryCoeffs()[patchI].clone()
|
|
);
|
|
internal.set
|
|
(
|
|
implicit,
|
|
matrix(i).internalCoeffs()[patchI].clone()
|
|
);
|
|
implicit++;
|
|
}
|
|
}
|
|
|
|
// Update non-implicit patches (re-order)
|
|
forAll(bpsi, patchI)
|
|
{
|
|
label globalPatchId = lduMeshPtr()->patchMap()[i][patchI];
|
|
if (globalPatchId != -1)
|
|
{
|
|
if (matrix(i).internalCoeffs().set(patchI))
|
|
{
|
|
internalCoeffs_.set
|
|
(
|
|
globalPatchId,
|
|
matrix(i).internalCoeffs()[patchI].clone()
|
|
);
|
|
}
|
|
|
|
if (matrix(i).boundaryCoeffs().set(patchI))
|
|
{
|
|
boundaryCoeffs_.set
|
|
(
|
|
globalPatchId,
|
|
matrix(i).boundaryCoeffs()[patchI].clone()
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Store implicit patches at the end of the list
|
|
implicit = 0;
|
|
forAll(bpsi, patchI)
|
|
{
|
|
label globalPatchId = lduMeshPtr()->patchMap()[i][patchI];
|
|
if (globalPatchId == -1)
|
|
{
|
|
const label implicitPatchId =
|
|
lduMeshPtr()->patchLocalToGlobalMap()[i][patchI];
|
|
|
|
internalCoeffs_.set
|
|
(
|
|
implicitPatchId, internal[implicit].clone()
|
|
);
|
|
boundaryCoeffs_.set
|
|
(
|
|
implicitPatchId, boundary[implicit].clone()
|
|
);
|
|
|
|
implicit++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// forAll(internalCoeffs_, patchI)
|
|
// {
|
|
// DebugVar(patchI)
|
|
// DebugVar(internalCoeffs_[patchI])
|
|
// DebugVar(boundaryCoeffs_[patchI])
|
|
// }
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::manipulateMatrix(direction cmp)
|
|
{
|
|
for (label i=0; i < nMatrices(); ++i)
|
|
{
|
|
forAll(psi(i).boundaryField(), patchI)
|
|
{
|
|
label globalPatchId = lduMeshPtr()->patchMap()[i][patchI];
|
|
|
|
if (globalPatchId == -1)
|
|
{
|
|
psi(i).boundaryFieldRef()[patchI].manipulateMatrix
|
|
(
|
|
*this,
|
|
i,
|
|
cmp
|
|
);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::transferFvMatrixCoeffs()
|
|
{
|
|
const labelListList& faceMap = lduMeshPtr()->faceMap();
|
|
const labelList& cellMap = lduMeshPtr()->cellOffsets();
|
|
|
|
label newFaces = lduMeshPtr()->lduAddr().upperAddr().size();
|
|
label newCells = lduMeshPtr()->lduAddr().size();
|
|
|
|
scalarField lowerAssemb(newFaces, Zero);
|
|
scalarField upperAssemb(newFaces, Zero);
|
|
scalarField diagAssemb(newCells, Zero);
|
|
Field<Type> sourceAssemb(newCells, Zero);
|
|
|
|
bool asymmetricAssemby = false;
|
|
for (label i=0; i < nMatrices(); ++i)
|
|
{
|
|
if (matrix(i).asymmetric())
|
|
{
|
|
asymmetricAssemby = true;
|
|
}
|
|
}
|
|
// Move append contents into intermediate list
|
|
for (label i=0; i < nMatrices(); ++i)
|
|
{
|
|
if (asymmetricAssemby)
|
|
{
|
|
const scalarField lowerSub(matrix(i).lower());
|
|
forAll(lowerSub, facei)
|
|
{
|
|
lowerAssemb[faceMap[i][facei]] = lowerSub[facei];
|
|
}
|
|
}
|
|
|
|
const scalarField upperSub(matrix(i).upper());
|
|
const scalarField diagSub(matrix(i).diag());
|
|
const Field<Type> sourceSub(matrix(i).source());
|
|
|
|
forAll(upperSub, facei)
|
|
{
|
|
upperAssemb[faceMap[i][facei]] = upperSub[facei];
|
|
}
|
|
|
|
forAll(diagSub, celli)
|
|
{
|
|
const label globalCelli = cellMap[i] + celli;
|
|
diagAssemb[globalCelli] = diagSub[celli];
|
|
sourceAssemb[globalCelli] = sourceSub[celli];
|
|
}
|
|
}
|
|
|
|
if (asymmetricAssemby)
|
|
{
|
|
lower().setSize(newFaces, Zero);
|
|
lower() = lowerAssemb;
|
|
}
|
|
upper().setSize(newFaces, Zero);
|
|
upper() = upperAssemb;
|
|
|
|
diag().setSize(newCells, Zero);
|
|
diag() = diagAssemb;
|
|
|
|
source().setSize(newCells, Zero);
|
|
source() = sourceAssemb;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::lduPrimitiveMeshAssembly* Foam::fvMatrix<Type>::lduMeshPtr()
|
|
{
|
|
const lduPrimitiveMeshAssembly* lduAssemMeshPtr =
|
|
psi_.mesh().thisDb().objectRegistry::template findObject
|
|
<
|
|
lduPrimitiveMeshAssembly
|
|
> (lduAssemblyName_);
|
|
|
|
return const_cast<lduPrimitiveMeshAssembly*>(lduAssemMeshPtr);
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
const Foam::lduPrimitiveMeshAssembly* Foam::fvMatrix<Type>::lduMeshPtr() const
|
|
{
|
|
return
|
|
(
|
|
psi_.mesh().thisDb().objectRegistry::template cfindObject
|
|
<
|
|
lduPrimitiveMeshAssembly
|
|
> (lduAssemblyName_)
|
|
);
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::createOrUpdateLduPrimitiveAssembly()
|
|
{
|
|
lduPrimitiveMeshAssembly* ptr = lduMeshPtr();
|
|
|
|
IOobject io
|
|
(
|
|
lduAssemblyName_,
|
|
psi_.mesh().time().timeName(),
|
|
psi_.mesh().thisDb(),
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
IOobject::REGISTER
|
|
);
|
|
|
|
UPtrList<lduMesh> uMeshPtr(nMatrices());
|
|
|
|
UPtrList<GeometricField<Type, fvPatchField, volMesh>>
|
|
uFieldPtr(nMatrices());
|
|
|
|
for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
|
|
{
|
|
const fvMesh& meshi = this->psi(fieldi).mesh();
|
|
uMeshPtr.set
|
|
(
|
|
fieldi,
|
|
&const_cast<fvMesh&>(meshi)
|
|
);
|
|
uFieldPtr.set(fieldi, &this->psi(fieldi));
|
|
}
|
|
|
|
if (!ptr)
|
|
{
|
|
lduPrimitiveMeshAssembly* lduAssemMeshPtr =
|
|
new lduPrimitiveMeshAssembly(io, uMeshPtr);
|
|
|
|
lduAssemMeshPtr->store();
|
|
lduAssemMeshPtr->update(uFieldPtr);
|
|
|
|
Info
|
|
<< "Creating lduPrimitiveAssembly: " << lduAssemblyName_ << endl;
|
|
}
|
|
else if
|
|
(
|
|
psi_.mesh().changing() && !psi_.mesh().upToDatePoints(*ptr)
|
|
)
|
|
{
|
|
// Clear losortPtr_, ownerStartPtr_, losortStartPtr_
|
|
ptr->lduAddr().clearOut();
|
|
ptr->update(uFieldPtr);
|
|
psi_.mesh().setUpToDatePoints(*ptr);
|
|
|
|
Info
|
|
<< "Updating lduPrimitiveAssembly: " << lduAssemblyName_ << endl;
|
|
}
|
|
else
|
|
{
|
|
Info
|
|
<< "Using lduPrimitiveAssembly: " << lduAssemblyName_ << endl;
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::setValues
|
|
(
|
|
const labelUList& cellLabels,
|
|
const Type& value
|
|
)
|
|
{
|
|
this->setValuesFromList(cellLabels, UniformList<Type>(value));
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::setValues
|
|
(
|
|
const labelUList& cellLabels,
|
|
const UList<Type>& values
|
|
)
|
|
{
|
|
this->setValuesFromList(cellLabels, values);
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::setValues
|
|
(
|
|
const labelUList& cellLabels,
|
|
const UIndirectList<Type>& values
|
|
)
|
|
{
|
|
this->setValuesFromList(cellLabels, values);
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::setReference
|
|
(
|
|
const label celli,
|
|
const Type& value,
|
|
const bool forceReference
|
|
)
|
|
{
|
|
if ((forceReference || psi_.needReference()) && celli >= 0)
|
|
{
|
|
source()[celli] += diag()[celli]*value;
|
|
diag()[celli] += diag()[celli];
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::setReferences
|
|
(
|
|
const labelUList& cellLabels,
|
|
const Type& value,
|
|
const bool forceReference
|
|
)
|
|
{
|
|
if (forceReference || psi_.needReference())
|
|
{
|
|
forAll(cellLabels, celli)
|
|
{
|
|
const label cellId = cellLabels[celli];
|
|
if (cellId >= 0)
|
|
{
|
|
source()[cellId] += diag()[cellId]*value;
|
|
diag()[cellId] += diag()[cellId];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::setReferences
|
|
(
|
|
const labelUList& cellLabels,
|
|
const UList<Type>& values,
|
|
const bool forceReference
|
|
)
|
|
{
|
|
if (forceReference || psi_.needReference())
|
|
{
|
|
forAll(cellLabels, celli)
|
|
{
|
|
const label cellId = cellLabels[celli];
|
|
if (cellId >= 0)
|
|
{
|
|
source()[cellId] += diag()[cellId]*values[celli];
|
|
diag()[cellId] += diag()[cellId];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::addFvMatrix(fvMatrix& matrix)
|
|
{
|
|
subMatrices_.append(matrix.clone());
|
|
++nMatrix_;
|
|
|
|
if (dimensions_ != matrix.dimensions())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "incompatible dimensions for matrix addition "
|
|
<< endl << " "
|
|
<< "[" << dimensions_ << " ] "
|
|
<< " [" << matrix.dimensions() << " ]"
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
for (label fieldi = 0; fieldi < nMatrices(); fieldi++)
|
|
{
|
|
if (checkImplicit(fieldi))
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
internalCoeffs_.clear();
|
|
boundaryCoeffs_.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::relax(const scalar alpha)
|
|
{
|
|
if (alpha <= 0)
|
|
{
|
|
return;
|
|
}
|
|
|
|
DebugInFunction
|
|
<< "Relaxing " << psi_.name() << " by " << alpha << endl;
|
|
|
|
Field<Type>& S = source();
|
|
scalarField& D = diag();
|
|
|
|
// Store the current unrelaxed diagonal for use in updating the source
|
|
scalarField D0(D);
|
|
|
|
// Calculate the sum-mag off-diagonal from the interior faces
|
|
scalarField sumOff(D.size(), Zero);
|
|
sumMagOffDiag(sumOff);
|
|
|
|
// Handle the boundary contributions to the diagonal
|
|
forAll(psi_.boundaryField(), patchi)
|
|
{
|
|
const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
|
|
|
|
if (ptf.size())
|
|
{
|
|
const labelUList& pa = lduAddr().patchAddr(patchi);
|
|
Field<Type>& iCoeffs = internalCoeffs_[patchi];
|
|
|
|
if (ptf.coupled())
|
|
{
|
|
const Field<Type>& pCoeffs = boundaryCoeffs_[patchi];
|
|
|
|
// For coupled boundaries add the diagonal and
|
|
// off-diagonal contributions
|
|
forAll(pa, face)
|
|
{
|
|
D[pa[face]] += component(iCoeffs[face], 0);
|
|
sumOff[pa[face]] += mag(component(pCoeffs[face], 0));
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// For non-coupled boundaries add the maximum magnitude diagonal
|
|
// contribution to ensure stability
|
|
forAll(pa, face)
|
|
{
|
|
D[pa[face]] += cmptMax(cmptMag(iCoeffs[face]));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
if (debug)
|
|
{
|
|
// Calculate amount of non-dominance.
|
|
label nNon = 0;
|
|
scalar maxNon = 0.0;
|
|
scalar sumNon = 0.0;
|
|
forAll(D, celli)
|
|
{
|
|
scalar d = (sumOff[celli] - D[celli])/mag(D[celli]);
|
|
|
|
if (d > 0)
|
|
{
|
|
nNon++;
|
|
maxNon = max(maxNon, d);
|
|
sumNon += d;
|
|
}
|
|
}
|
|
|
|
reduce(nNon, sumOp<label>(), UPstream::msgType(), psi_.mesh().comm());
|
|
reduce
|
|
(
|
|
maxNon,
|
|
maxOp<scalar>(),
|
|
UPstream::msgType(),
|
|
psi_.mesh().comm()
|
|
);
|
|
reduce
|
|
(
|
|
sumNon,
|
|
sumOp<scalar>(),
|
|
UPstream::msgType(),
|
|
psi_.mesh().comm()
|
|
);
|
|
sumNon /= returnReduce
|
|
(
|
|
D.size(),
|
|
sumOp<label>(),
|
|
UPstream::msgType(),
|
|
psi_.mesh().comm()
|
|
);
|
|
|
|
InfoInFunction
|
|
<< "Matrix dominance test for " << psi_.name() << nl
|
|
<< " number of non-dominant cells : " << nNon << nl
|
|
<< " maximum relative non-dominance : " << maxNon << nl
|
|
<< " average relative non-dominance : " << sumNon << nl
|
|
<< endl;
|
|
}
|
|
|
|
|
|
// Ensure the matrix is diagonally dominant...
|
|
// Assumes that the central coefficient is positive and ensures it is
|
|
forAll(D, celli)
|
|
{
|
|
D[celli] = max(mag(D[celli]), sumOff[celli]);
|
|
}
|
|
|
|
// ... then relax
|
|
D /= alpha;
|
|
|
|
// Now remove the diagonal contribution from coupled boundaries
|
|
forAll(psi_.boundaryField(), patchi)
|
|
{
|
|
const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
|
|
|
|
if (ptf.size())
|
|
{
|
|
const labelUList& pa = lduAddr().patchAddr(patchi);
|
|
Field<Type>& iCoeffs = internalCoeffs_[patchi];
|
|
|
|
if (ptf.coupled())
|
|
{
|
|
forAll(pa, face)
|
|
{
|
|
D[pa[face]] -= component(iCoeffs[face], 0);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
forAll(pa, face)
|
|
{
|
|
D[pa[face]] -= cmptMin(iCoeffs[face]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Finally add the relaxation contribution to the source.
|
|
S += (D - D0)*psi_.primitiveField();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::relax()
|
|
{
|
|
word name = psi_.select
|
|
(
|
|
psi_.mesh().data::template getOrDefault<bool>
|
|
("finalIteration", false)
|
|
);
|
|
|
|
if (psi_.mesh().relaxEquation(name))
|
|
{
|
|
relax(psi_.mesh().equationRelaxationFactor(name));
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::boundaryManipulate
|
|
(
|
|
typename GeometricField<Type, fvPatchField, volMesh>::
|
|
Boundary& bFields
|
|
)
|
|
{
|
|
forAll(bFields, patchi)
|
|
{
|
|
bFields[patchi].manipulateMatrix(*this);
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::scalarField> Foam::fvMatrix<Type>::D() const
|
|
{
|
|
auto tdiag = tmp<scalarField>::New(diag());
|
|
addCmptAvBoundaryDiag(tdiag.ref());
|
|
return tdiag;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::Field<Type>> Foam::fvMatrix<Type>::DD() const
|
|
{
|
|
tmp<Field<Type>> tdiag(pTraits<Type>::one*diag());
|
|
|
|
forAll(psi_.boundaryField(), patchi)
|
|
{
|
|
const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
|
|
|
|
if (!ptf.coupled() && ptf.size())
|
|
{
|
|
addToInternalField
|
|
(
|
|
lduAddr().patchAddr(patchi),
|
|
internalCoeffs_[patchi],
|
|
tdiag.ref()
|
|
);
|
|
}
|
|
}
|
|
|
|
return tdiag;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const
|
|
{
|
|
auto tAphi = volScalarField::New
|
|
(
|
|
"A(" + psi_.name() + ')',
|
|
psi_.mesh(),
|
|
dimensions_/psi_.dimensions()/dimVol,
|
|
fvPatchFieldBase::extrapolatedCalculatedType()
|
|
);
|
|
|
|
tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V();
|
|
tAphi.ref().correctBoundaryConditions();
|
|
|
|
return tAphi;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
|
|
Foam::fvMatrix<Type>::H() const
|
|
{
|
|
auto tHphi = GeometricField<Type, fvPatchField, volMesh>::New
|
|
(
|
|
"H(" + psi_.name() + ')',
|
|
psi_.mesh(),
|
|
dimensions_/dimVol,
|
|
fvPatchFieldBase::extrapolatedCalculatedType()
|
|
);
|
|
auto& Hphi = tHphi.ref();
|
|
|
|
// Loop over field components
|
|
for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
|
|
{
|
|
scalarField psiCmpt(psi_.primitiveField().component(cmpt));
|
|
|
|
scalarField boundaryDiagCmpt(psi_.size(), Zero);
|
|
addBoundaryDiag(boundaryDiagCmpt, cmpt);
|
|
boundaryDiagCmpt.negate();
|
|
addCmptAvBoundaryDiag(boundaryDiagCmpt);
|
|
|
|
Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
|
|
}
|
|
|
|
Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
|
|
addBoundarySource(Hphi.primitiveFieldRef());
|
|
|
|
Hphi.primitiveFieldRef() /= psi_.mesh().V();
|
|
Hphi.correctBoundaryConditions();
|
|
|
|
typename Type::labelType validComponents
|
|
(
|
|
psi_.mesh().template validComponents<Type>()
|
|
);
|
|
|
|
for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
|
|
{
|
|
if (validComponents[cmpt] == -1)
|
|
{
|
|
Hphi.replace
|
|
(
|
|
cmpt,
|
|
dimensionedScalar(Hphi.dimensions(), Zero)
|
|
);
|
|
}
|
|
}
|
|
|
|
return tHphi;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::H1() const
|
|
{
|
|
auto tH1 = volScalarField::New
|
|
(
|
|
"H(1)",
|
|
psi_.mesh(),
|
|
dimensions_/(dimVol*psi_.dimensions()),
|
|
fvPatchFieldBase::extrapolatedCalculatedType()
|
|
);
|
|
auto& H1_ = tH1.ref();
|
|
|
|
H1_.primitiveFieldRef() = lduMatrix::H1();
|
|
|
|
forAll(psi_.boundaryField(), patchi)
|
|
{
|
|
const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
|
|
|
|
if (ptf.coupled() && ptf.size())
|
|
{
|
|
addToInternalField
|
|
(
|
|
lduAddr().patchAddr(patchi),
|
|
boundaryCoeffs_[patchi].component(0),
|
|
H1_
|
|
);
|
|
}
|
|
}
|
|
|
|
H1_.primitiveFieldRef() /= psi_.mesh().V();
|
|
H1_.correctBoundaryConditions();
|
|
|
|
return tH1;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::GeometricField<Type, Foam::fvsPatchField, Foam::surfaceMesh>>
|
|
Foam::fvMatrix<Type>::
|
|
flux() const
|
|
{
|
|
if (!psi_.mesh().fluxRequired(psi_.name()))
|
|
{
|
|
FatalErrorInFunction
|
|
<< "flux requested but " << psi_.name()
|
|
<< " not specified in the fluxRequired sub-dictionary"
|
|
" of fvSchemes."
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
if (nMatrices() > 1)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Flux requested but " << psi_.name()
|
|
<< " can't handle multiple fvMatrix."
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
auto tfieldFlux = GeometricField<Type, fvsPatchField, surfaceMesh>::New
|
|
(
|
|
"flux(" + psi_.name() + ')',
|
|
psi_.mesh(),
|
|
dimensions()
|
|
);
|
|
auto& fieldFlux = tfieldFlux.ref();
|
|
|
|
fieldFlux.setOriented();
|
|
|
|
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
|
|
{
|
|
fieldFlux.primitiveFieldRef().replace
|
|
(
|
|
cmpt,
|
|
lduMatrix::faceH(psi_.primitiveField().component(cmpt))
|
|
);
|
|
}
|
|
|
|
FieldField<Field, Type> InternalContrib = internalCoeffs_;
|
|
|
|
label fieldi = 0;
|
|
if (!useImplicit_)
|
|
{
|
|
forAll(InternalContrib, patchi)
|
|
{
|
|
InternalContrib[patchi] =
|
|
cmptMultiply
|
|
(
|
|
InternalContrib[patchi],
|
|
psi_.boundaryField()[patchi].patchInternalField()
|
|
);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
FieldField<Field, Type> fluxInternalContrib(internalCoeffs_);
|
|
|
|
mapContributions(fieldi, fluxInternalContrib, InternalContrib, true);
|
|
}
|
|
|
|
FieldField<Field, Type> NeighbourContrib = boundaryCoeffs_;
|
|
|
|
if (!useImplicit_)
|
|
{
|
|
forAll(NeighbourContrib, patchi)
|
|
{
|
|
if (psi_.boundaryField()[patchi].coupled())
|
|
{
|
|
NeighbourContrib[patchi] =
|
|
cmptMultiply
|
|
(
|
|
NeighbourContrib[patchi],
|
|
psi_.boundaryField()[patchi].patchNeighbourField()
|
|
);
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
FieldField<Field, Type> fluxBoundaryContrib(boundaryCoeffs_);
|
|
|
|
mapContributions(fieldi, fluxBoundaryContrib, NeighbourContrib, false);
|
|
}
|
|
|
|
typename GeometricField<Type, fvsPatchField, surfaceMesh>::
|
|
Boundary& ffbf = fieldFlux.boundaryFieldRef();
|
|
|
|
forAll(ffbf, patchi)
|
|
{
|
|
ffbf[patchi] = InternalContrib[patchi] - NeighbourContrib[patchi];
|
|
//DebugVar(gSum(ffbf[patchi]))
|
|
}
|
|
|
|
if (faceFluxCorrectionPtr_)
|
|
{
|
|
fieldFlux += *faceFluxCorrectionPtr_;
|
|
}
|
|
|
|
return tfieldFlux;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
const Foam::dictionary& Foam::fvMatrix<Type>::solverDict() const
|
|
{
|
|
return psi_.mesh().solverDict
|
|
(
|
|
psi_.select
|
|
(
|
|
psi_.mesh().data::template getOrDefault<bool>
|
|
("finalIteration", false)
|
|
)
|
|
);
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator=(const fvMatrix<Type>& fvmv)
|
|
{
|
|
if (this == &fvmv)
|
|
{
|
|
return; // Self-assignment is a no-op
|
|
}
|
|
|
|
if (&psi_ != &(fvmv.psi_))
|
|
{
|
|
FatalErrorInFunction
|
|
<< "different fields"
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
dimensions_ = fvmv.dimensions_;
|
|
lduMatrix::operator=(fvmv);
|
|
source_ = fvmv.source_;
|
|
internalCoeffs_ = fvmv.internalCoeffs_;
|
|
boundaryCoeffs_ = fvmv.boundaryCoeffs_;
|
|
|
|
if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
|
|
{
|
|
*faceFluxCorrectionPtr_ = *fvmv.faceFluxCorrectionPtr_;
|
|
}
|
|
else if (fvmv.faceFluxCorrectionPtr_)
|
|
{
|
|
faceFluxCorrectionPtr_ =
|
|
new GeometricField<Type, fvsPatchField, surfaceMesh>
|
|
(*fvmv.faceFluxCorrectionPtr_);
|
|
}
|
|
|
|
useImplicit_ = fvmv.useImplicit_;
|
|
lduAssemblyName_ = fvmv.lduAssemblyName_;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator=(const tmp<fvMatrix<Type>>& tfvmv)
|
|
{
|
|
operator=(tfvmv());
|
|
tfvmv.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::negate()
|
|
{
|
|
lduMatrix::negate();
|
|
source_.negate();
|
|
internalCoeffs_.negate();
|
|
boundaryCoeffs_.negate();
|
|
|
|
if (faceFluxCorrectionPtr_)
|
|
{
|
|
faceFluxCorrectionPtr_->negate();
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator+=(const fvMatrix<Type>& fvmv)
|
|
{
|
|
checkMethod(*this, fvmv, "+=");
|
|
|
|
dimensions_ += fvmv.dimensions_;
|
|
lduMatrix::operator+=(fvmv);
|
|
source_ += fvmv.source_;
|
|
internalCoeffs_ += fvmv.internalCoeffs_;
|
|
boundaryCoeffs_ += fvmv.boundaryCoeffs_;
|
|
|
|
useImplicit_ = fvmv.useImplicit_;
|
|
lduAssemblyName_ = fvmv.lduAssemblyName_;
|
|
nMatrix_ = fvmv.nMatrix_;
|
|
|
|
if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
|
|
{
|
|
*faceFluxCorrectionPtr_ += *fvmv.faceFluxCorrectionPtr_;
|
|
}
|
|
else if (fvmv.faceFluxCorrectionPtr_)
|
|
{
|
|
faceFluxCorrectionPtr_ = new
|
|
GeometricField<Type, fvsPatchField, surfaceMesh>
|
|
(
|
|
*fvmv.faceFluxCorrectionPtr_
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator+=(const tmp<fvMatrix<Type>>& tfvmv)
|
|
{
|
|
operator+=(tfvmv());
|
|
tfvmv.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator-=(const fvMatrix<Type>& fvmv)
|
|
{
|
|
checkMethod(*this, fvmv, "-=");
|
|
|
|
dimensions_ -= fvmv.dimensions_;
|
|
lduMatrix::operator-=(fvmv);
|
|
source_ -= fvmv.source_;
|
|
internalCoeffs_ -= fvmv.internalCoeffs_;
|
|
boundaryCoeffs_ -= fvmv.boundaryCoeffs_;
|
|
|
|
useImplicit_ = fvmv.useImplicit_;
|
|
lduAssemblyName_ = fvmv.lduAssemblyName_;
|
|
nMatrix_ = fvmv.nMatrix_;
|
|
|
|
if (faceFluxCorrectionPtr_ && fvmv.faceFluxCorrectionPtr_)
|
|
{
|
|
*faceFluxCorrectionPtr_ -= *fvmv.faceFluxCorrectionPtr_;
|
|
}
|
|
else if (fvmv.faceFluxCorrectionPtr_)
|
|
{
|
|
faceFluxCorrectionPtr_ =
|
|
new GeometricField<Type, fvsPatchField, surfaceMesh>
|
|
(-*fvmv.faceFluxCorrectionPtr_);
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator-=(const tmp<fvMatrix<Type>>& tfvmv)
|
|
{
|
|
operator-=(tfvmv());
|
|
tfvmv.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator+=
|
|
(
|
|
const DimensionedField<Type, volMesh>& su
|
|
)
|
|
{
|
|
checkMethod(*this, su, "+=");
|
|
source() -= su.mesh().V()*su.field();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator+=
|
|
(
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu
|
|
)
|
|
{
|
|
operator+=(tsu());
|
|
tsu.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator+=
|
|
(
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
|
|
)
|
|
{
|
|
operator+=(tsu());
|
|
tsu.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator-=
|
|
(
|
|
const DimensionedField<Type, volMesh>& su
|
|
)
|
|
{
|
|
checkMethod(*this, su, "-=");
|
|
source() += su.mesh().V()*su.field();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator-=
|
|
(
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu
|
|
)
|
|
{
|
|
operator-=(tsu());
|
|
tsu.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator-=
|
|
(
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
|
|
)
|
|
{
|
|
operator-=(tsu());
|
|
tsu.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator+=
|
|
(
|
|
const dimensioned<Type>& su
|
|
)
|
|
{
|
|
source() -= psi().mesh().V()*su;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator-=
|
|
(
|
|
const dimensioned<Type>& su
|
|
)
|
|
{
|
|
source() += psi().mesh().V()*su;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator*=
|
|
(
|
|
const volScalarField::Internal& dsf
|
|
)
|
|
{
|
|
dimensions_ *= dsf.dimensions();
|
|
lduMatrix::operator*=(dsf.field());
|
|
source_ *= dsf.field();
|
|
|
|
forAll(boundaryCoeffs_, patchi)
|
|
{
|
|
scalarField pisf
|
|
(
|
|
dsf.mesh().boundary()[patchi].patchInternalField(dsf.field())
|
|
);
|
|
|
|
internalCoeffs_[patchi] *= pisf;
|
|
boundaryCoeffs_[patchi] *= pisf;
|
|
}
|
|
|
|
if (faceFluxCorrectionPtr_)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "cannot scale a matrix containing a faceFluxCorrection"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator*=
|
|
(
|
|
const tmp<volScalarField::Internal>& tfld
|
|
)
|
|
{
|
|
operator*=(tfld());
|
|
tfld.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator*=
|
|
(
|
|
const tmp<volScalarField>& tfld
|
|
)
|
|
{
|
|
operator*=(tfld());
|
|
tfld.clear();
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::fvMatrix<Type>::operator*=
|
|
(
|
|
const dimensioned<scalar>& ds
|
|
)
|
|
{
|
|
dimensions_ *= ds.dimensions();
|
|
lduMatrix::operator*=(ds.value());
|
|
source_ *= ds.value();
|
|
internalCoeffs_ *= ds.value();
|
|
boundaryCoeffs_ *= ds.value();
|
|
|
|
if (faceFluxCorrectionPtr_)
|
|
{
|
|
*faceFluxCorrectionPtr_ *= ds.value();
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
|
|
|
|
template<class Type>
|
|
void Foam::checkMethod
|
|
(
|
|
const fvMatrix<Type>& mat1,
|
|
const fvMatrix<Type>& mat2,
|
|
const char* op
|
|
)
|
|
{
|
|
if (&mat1.psi() != &mat2.psi())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Incompatible fields for operation\n "
|
|
<< "[" << mat1.psi().name() << "] "
|
|
<< op
|
|
<< " [" << mat2.psi().name() << "]"
|
|
<< abort(FatalError);
|
|
}
|
|
|
|
if
|
|
(
|
|
dimensionSet::checking()
|
|
&& mat1.dimensions() != mat2.dimensions()
|
|
)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Incompatible dimensions for operation\n "
|
|
<< "[" << mat1.psi().name() << mat1.dimensions()/dimVolume << " ] "
|
|
<< op
|
|
<< " [" << mat2.psi().name() << mat2.dimensions()/dimVolume << " ]"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::checkMethod
|
|
(
|
|
const fvMatrix<Type>& mat,
|
|
const DimensionedField<Type, volMesh>& fld,
|
|
const char* op
|
|
)
|
|
{
|
|
if
|
|
(
|
|
dimensionSet::checking()
|
|
&& mat.dimensions()/dimVolume != fld.dimensions()
|
|
)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Incompatible dimensions for operation\n "
|
|
<< "[" << mat.psi().name() << mat.dimensions()/dimVolume << " ] "
|
|
<< op
|
|
<< " [" << fld.name() << fld.dimensions() << " ]"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
void Foam::checkMethod
|
|
(
|
|
const fvMatrix<Type>& mat,
|
|
const dimensioned<Type>& dt,
|
|
const char* op
|
|
)
|
|
{
|
|
if
|
|
(
|
|
dimensionSet::checking()
|
|
&& mat.dimensions()/dimVolume != dt.dimensions()
|
|
)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Incompatible dimensions for operation\n "
|
|
<< "[" << mat.psi().name() << mat.dimensions()/dimVolume << " ] "
|
|
<< op
|
|
<< " [" << dt.name() << dt.dimensions() << " ]"
|
|
<< abort(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::SolverPerformance<Type> Foam::solve
|
|
(
|
|
fvMatrix<Type>& fvm,
|
|
const dictionary& solverControls
|
|
)
|
|
{
|
|
return fvm.solve(solverControls);
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::SolverPerformance<Type> Foam::solve
|
|
(
|
|
const tmp<fvMatrix<Type>>& tmat,
|
|
const dictionary& solverControls
|
|
)
|
|
{
|
|
SolverPerformance<Type> solverPerf(tmat.constCast().solve(solverControls));
|
|
|
|
tmat.clear();
|
|
|
|
return solverPerf;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::SolverPerformance<Type> Foam::solve(fvMatrix<Type>& fvm)
|
|
{
|
|
return fvm.solve();
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::SolverPerformance<Type> Foam::solve(const tmp<fvMatrix<Type>>& tmat)
|
|
{
|
|
SolverPerformance<Type> solverPerf(tmat.constCast().solve());
|
|
|
|
tmat.clear();
|
|
|
|
return solverPerf;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::correction
|
|
(
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
tmp<Foam::fvMatrix<Type>> tAcorr = A - (A & A.psi());
|
|
|
|
// Delete the faceFluxCorrection from the correction matrix
|
|
// as it does not have a clear meaning or purpose
|
|
deleteDemandDrivenData(tAcorr.ref().faceFluxCorrectionPtr());
|
|
|
|
return tAcorr;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::correction
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
tmp<Foam::fvMatrix<Type>> tAcorr = tA - (tA() & tA().psi());
|
|
|
|
// Delete the faceFluxCorrection from the correction matrix
|
|
// as it does not have a clear meaning or purpose
|
|
deleteDemandDrivenData(tAcorr.ref().faceFluxCorrectionPtr());
|
|
|
|
return tAcorr;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const fvMatrix<Type>& B
|
|
)
|
|
{
|
|
checkMethod(A, B, "==");
|
|
return (A - B);
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const fvMatrix<Type>& B
|
|
)
|
|
{
|
|
checkMethod(tA(), B, "==");
|
|
return (tA - B);
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<fvMatrix<Type>>& tB
|
|
)
|
|
{
|
|
checkMethod(A, tB(), "==");
|
|
return (A - tB);
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<fvMatrix<Type>>& tB
|
|
)
|
|
{
|
|
checkMethod(tA(), tB(), "==");
|
|
return (tA - tB);
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const DimensionedField<Type, volMesh>& su
|
|
)
|
|
{
|
|
checkMethod(A, su, "==");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() += su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "==");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() += tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "==");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const DimensionedField<Type, volMesh>& su
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "==");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() += su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "==");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() += tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "==");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const dimensioned<Type>& su
|
|
)
|
|
{
|
|
checkMethod(A, su, "==");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() += A.psi().mesh().V()*su.value();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const dimensioned<Type>& su
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "==");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() += tC().psi().mesh().V()*su.value();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const Foam::zero
|
|
)
|
|
{
|
|
return A;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const Foam::zero
|
|
)
|
|
{
|
|
return tA;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().negate();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().negate();
|
|
return tC;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const fvMatrix<Type>& B
|
|
)
|
|
{
|
|
checkMethod(A, B, "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref() += B;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const fvMatrix<Type>& B
|
|
)
|
|
{
|
|
checkMethod(tA(), B, "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref() += B;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<fvMatrix<Type>>& tB
|
|
)
|
|
{
|
|
checkMethod(A, tB(), "+");
|
|
tmp<fvMatrix<Type>> tC(tB.ptr());
|
|
tC.ref() += A;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<fvMatrix<Type>>& tB
|
|
)
|
|
{
|
|
checkMethod(tA(), tB(), "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref() += tB();
|
|
tB.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const DimensionedField<Type, volMesh>& su
|
|
)
|
|
{
|
|
checkMethod(A, su, "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() -= su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const DimensionedField<Type, volMesh>& su
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() -= su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const DimensionedField<Type, volMesh>& su,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
checkMethod(A, su, "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() -= su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const DimensionedField<Type, volMesh>& su,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() -= su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const fvMatrix<Type>& B
|
|
)
|
|
{
|
|
checkMethod(A, B, "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref() -= B;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const fvMatrix<Type>& B
|
|
)
|
|
{
|
|
checkMethod(tA(), B, "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref() -= B;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<fvMatrix<Type>>& tB
|
|
)
|
|
{
|
|
checkMethod(A, tB(), "-");
|
|
tmp<fvMatrix<Type>> tC(tB.ptr());
|
|
tC.ref() -= A;
|
|
tC.ref().negate();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<fvMatrix<Type>>& tB
|
|
)
|
|
{
|
|
checkMethod(tA(), tB(), "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref() -= tB();
|
|
tB.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const DimensionedField<Type, volMesh>& su
|
|
)
|
|
{
|
|
checkMethod(A, su, "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() += su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() += tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const DimensionedField<Type, volMesh>& su
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() += su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() += tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const DimensionedField<Type, volMesh>& su,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
checkMethod(A, su, "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().negate();
|
|
tC.ref().source() -= su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().negate();
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
checkMethod(A, tsu(), "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().negate();
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const DimensionedField<Type, volMesh>& su,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().negate();
|
|
tC.ref().source() -= su.mesh().V()*su.field();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<DimensionedField<Type, volMesh>>& tsu,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().negate();
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().field();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
checkMethod(tA(), tsu(), "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().negate();
|
|
tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField();
|
|
tsu.clear();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const dimensioned<Type>& su
|
|
)
|
|
{
|
|
checkMethod(A, su, "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() -= su.value()*A.psi().mesh().V();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const dimensioned<Type>& su
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() -= su.value()*tC().psi().mesh().V();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const dimensioned<Type>& su,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
checkMethod(A, su, "+");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() -= su.value()*A.psi().mesh().V();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+
|
|
(
|
|
const dimensioned<Type>& su,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "+");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() -= su.value()*tC().psi().mesh().V();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const fvMatrix<Type>& A,
|
|
const dimensioned<Type>& su
|
|
)
|
|
{
|
|
checkMethod(A, su, "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().source() += su.value()*tC().psi().mesh().V();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const tmp<fvMatrix<Type>>& tA,
|
|
const dimensioned<Type>& su
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().source() += su.value()*tC().psi().mesh().V();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const dimensioned<Type>& su,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
checkMethod(A, su, "-");
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref().negate();
|
|
tC.ref().source() -= su.value()*A.psi().mesh().V();
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator-
|
|
(
|
|
const dimensioned<Type>& su,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
checkMethod(tA(), su, "-");
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref().negate();
|
|
tC.ref().source() -= su.value()*tC().psi().mesh().V();
|
|
return tC;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
|
|
(
|
|
const volScalarField::Internal& dsf,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref() *= dsf;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
|
|
(
|
|
const tmp<volScalarField::Internal>& tdsf,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref() *= tdsf;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
|
|
(
|
|
const tmp<volScalarField>& tvsf,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref() *= tvsf;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
|
|
(
|
|
const volScalarField::Internal& dsf,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref() *= dsf;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
|
|
(
|
|
const tmp<volScalarField::Internal>& tdsf,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref() *= tdsf;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
|
|
(
|
|
const tmp<volScalarField>& tvsf,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref() *= tvsf;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
|
|
(
|
|
const dimensioned<scalar>& ds,
|
|
const fvMatrix<Type>& A
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
|
|
tC.ref() *= ds;
|
|
return tC;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator*
|
|
(
|
|
const dimensioned<scalar>& ds,
|
|
const tmp<fvMatrix<Type>>& tA
|
|
)
|
|
{
|
|
tmp<fvMatrix<Type>> tC(tA.ptr());
|
|
tC.ref() *= ds;
|
|
return tC;
|
|
}
|
|
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
|
|
Foam::operator&
|
|
(
|
|
const fvMatrix<Type>& M,
|
|
const DimensionedField<Type, volMesh>& psi
|
|
)
|
|
{
|
|
auto tMphi = GeometricField<Type, fvPatchField, volMesh>::New
|
|
(
|
|
"M&" + psi.name(),
|
|
psi.mesh(),
|
|
M.dimensions()/dimVol,
|
|
fvPatchFieldBase::extrapolatedCalculatedType()
|
|
);
|
|
auto& Mphi = tMphi.ref();
|
|
|
|
// Loop over field components
|
|
if (M.hasDiag())
|
|
{
|
|
for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
|
|
{
|
|
scalarField psiCmpt(psi.field().component(cmpt));
|
|
scalarField boundaryDiagCmpt(M.diag());
|
|
M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
|
|
Mphi.primitiveFieldRef().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
Mphi.primitiveFieldRef() = Zero;
|
|
}
|
|
|
|
Mphi.primitiveFieldRef() += M.lduMatrix::H(psi.field()) + M.source();
|
|
M.addBoundarySource(Mphi.primitiveFieldRef());
|
|
|
|
Mphi.primitiveFieldRef() /= -psi.mesh().V();
|
|
Mphi.correctBoundaryConditions();
|
|
|
|
return tMphi;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
|
|
Foam::operator&
|
|
(
|
|
const fvMatrix<Type>& M,
|
|
const tmp<DimensionedField<Type, volMesh>>& tpsi
|
|
)
|
|
{
|
|
tmp<GeometricField<Type, fvPatchField, volMesh>> tMpsi = M & tpsi();
|
|
tpsi.clear();
|
|
return tMpsi;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
|
|
Foam::operator&
|
|
(
|
|
const fvMatrix<Type>& M,
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tpsi
|
|
)
|
|
{
|
|
tmp<GeometricField<Type, fvPatchField, volMesh>> tMpsi = M & tpsi();
|
|
tpsi.clear();
|
|
return tMpsi;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
|
|
Foam::operator&
|
|
(
|
|
const tmp<fvMatrix<Type>>& tM,
|
|
const DimensionedField<Type, volMesh>& psi
|
|
)
|
|
{
|
|
tmp<GeometricField<Type, fvPatchField, volMesh>> tMpsi = tM() & psi;
|
|
tM.clear();
|
|
return tMpsi;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
|
|
Foam::operator&
|
|
(
|
|
const tmp<fvMatrix<Type>>& tM,
|
|
const tmp<DimensionedField<Type, volMesh>>& tpsi
|
|
)
|
|
{
|
|
tmp<GeometricField<Type, fvPatchField, volMesh>> tMpsi = tM() & tpsi();
|
|
tM.clear();
|
|
tpsi.clear();
|
|
return tMpsi;
|
|
}
|
|
|
|
template<class Type>
|
|
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
|
|
Foam::operator&
|
|
(
|
|
const tmp<fvMatrix<Type>>& tM,
|
|
const tmp<GeometricField<Type, fvPatchField, volMesh>>& tpsi
|
|
)
|
|
{
|
|
tmp<GeometricField<Type, fvPatchField, volMesh>> tMpsi = tM() & tpsi();
|
|
tM.clear();
|
|
tpsi.clear();
|
|
return tMpsi;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * //
|
|
|
|
template<class Type>
|
|
Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
|
|
{
|
|
os << static_cast<const lduMatrix&>(fvm) << nl
|
|
<< fvm.dimensions_ << nl
|
|
<< fvm.source_ << nl
|
|
<< fvm.internalCoeffs_ << nl
|
|
<< fvm.boundaryCoeffs_ << endl;
|
|
|
|
os.check(FUNCTION_NAME);
|
|
|
|
return os;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
|
|
|
|
#include "fvMatrixSolve.C"
|
|
|
|
// ************************************************************************* //
|