33 template<
class Type,
class DType,
class LUType>
51 template<
class Type,
class DType,
class LUType>
54 const Field<LUType>& Lower =
const_cast<const LduMatrix&
>(*this).lower();
55 const Field<LUType>& Upper =
const_cast<const LduMatrix&
>(*this).upper();
56 Field<DType>& Diag =
diag();
61 for (label face=0; face<l.size(); face++)
63 Diag[l[face]] -= Lower[face];
69 template<
class Type,
class DType,
class LUType>
75 const Field<LUType>& Lower =
const_cast<const LduMatrix&
>(*this).lower();
76 const Field<LUType>& Upper =
const_cast<const LduMatrix&
>(*this).upper();
81 for (label face = 0; face < l.size(); face++)
83 sumOff[u[face]] +=
cmptMag(Lower[face]);
89 template<
class Type,
class DType,
class LUType>
93 tmp<Field<Type>> tHpsi
95 new Field<Type>(lduAddr().size(),
Zero)
98 if (lowerPtr_ || upperPtr_)
100 Field<Type> & Hpsi = tHpsi();
102 Type* __restrict__ HpsiPtr = Hpsi.begin();
104 const Type* __restrict__ psiPtr =
psi.begin();
106 const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
107 const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
109 const LUType* __restrict__ lowerPtr =
lower().begin();
110 const LUType* __restrict__ upperPtr =
upper().begin();
112 const label nFaces =
upper().size();
114 for (label face=0; face<nFaces; face++)
116 HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
117 HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
124 template<
class Type,
class DType,
class LUType>
128 tmp<Field<Type>> tHpsi(
H(tpsi()));
134 template<
class Type,
class DType,
class LUType>
157 template<
class Type,
class DType,
class LUType>
169 template<
class Type,
class DType,
class LUType>
204 source() =
A.source();
207 interfacesUpper_ =
A.interfacesUpper_;
208 interfacesLower_ =
A.interfacesLower_;
212 template<
class Type,
class DType,
class LUType>
232 sourcePtr_->negate();
240 template<
class Type,
class DType,
class LUType>
250 source() +=
A.source();
253 if (symmetric() &&
A.symmetric())
257 else if (symmetric() &&
A.asymmetric())
271 else if (asymmetric() &&
A.symmetric())
285 else if (asymmetric() &&
A.asymmetric())
302 else if (
A.diagonal())
308 <<
"Unknown matrix type combination" 312 interfacesUpper_ +=
A.interfacesUpper_;
313 interfacesLower_ +=
A.interfacesLower_;
317 template<
class Type,
class DType,
class LUType>
327 source() -=
A.source();
330 if (symmetric() &&
A.symmetric())
334 else if (symmetric() &&
A.asymmetric())
348 else if (asymmetric() &&
A.symmetric())
362 else if (asymmetric() &&
A.asymmetric())
379 else if (
A.diagonal())
385 <<
"Unknown matrix type combination" 389 interfacesUpper_ -=
A.interfacesUpper_;
390 interfacesLower_ -=
A.interfacesLower_;
394 template<
class Type,
class DType,
class LUType>
412 if (symmetric() || asymmetric())
420 for (label face=0; face<
upper.size(); face++)
422 upper[face] *= sf[l[face]];
425 for (label face=0; face<
lower.size(); face++)
427 lower[face] *= sf[u[face]];
432 <<
"Scaling a matrix by scalarField is not currently supported\n" 433 "because scaling interfacesUpper_ and interfacesLower_ " 434 "require special transfers" 442 template<
class Type,
class DType,
class LUType>
465 interfacesUpper_ *=
s;
466 interfacesLower_ *=
s;
void size(const label n)
Older name for setAddressableSize.
A face is a list of labels corresponding to mesh vertices.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void sumMagOffDiag(Field< LUType > &sumOff) const
Field< LUType > & lower()
string upper(const std::string &s)
Return string copy transformed with std::toupper on each character.
UList< label > labelUList
A UList of labels.
volScalarField H(IOobject("H", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, dimensionedScalar(dimLength, Zero))
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
tmp< Field< Type > > H(const Field< Type > &) const
void operator*=(const scalarField &)
Field< LUType > & upper()
errorManip< error > abort(error &err)
void operator-=(const LduMatrix< Type, DType, LUType > &)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
void negate(FieldField< Field, Type > &res, const FieldField< Field, Type > &f)
string lower(const std::string &s)
Return string copy transformed with std::tolower on each character.
LduMatrix is a general matrix class in which the coefficients are stored as three arrays...
void operator+=(const LduMatrix< Type, DType, LUType > &)
tmp< Field< Type > > faceH(const Field< Type > &) const
const volScalarField & psi
A class for managing temporary objects.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
void operator=(const LduMatrix< Type, DType, LUType > &)
static constexpr const zero Zero
Global zero (0)