zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-et' into 'develop'

Browse Source 
ENH: Expressions: basic expression template support

See merge request Development/openfoam!763
This commit is contained in:
Andrew Heather
2025-10-22 15:45:08 +00:00
parent 78cf2ce5db b0066fc550
commit f186b0098b
25 changed files with 7801 additions and 89 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
applications/test/expressionTemplates/Make/files Normal file
View File

@ -0,0 +1,5 @@
Test-expressionTemplates-volFields.C
/* Test-expressionTemplates-pointFields.C */
/* Test-expressionTemplates-Fields.C */
EXE = $(FOAM_USER_APPBIN)/Test-expressionTemplates

11
applications/test/expressionTemplates/Make/options Normal file
View File

@ -0,0 +1,11 @@
EXE_INC = \
-O0 \
-I$(FOAM_SRC)/meshTools/lnInclude \
-I$(FOAM_SRC)/finiteVolume/lnInclude \
EXE_LIBS = \
/* -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms */ \
/* -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib */ \
/* -rpath /Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib/$(FOAM_MPI) */\
-L/Volumes/case_sensitive/OpenFOAM/work/feature-et/platforms/darwin64GccDPInt32Opt/lib \
-lfiniteVolume

126
applications/test/expressionTemplates/Test-expressionTemplates-pointFields.C Normal file
View File

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2025 M. Janssens
-------------------------------------------------------------------------------
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/>.
Application
Test-expressionTemplates
\*---------------------------------------------------------------------------*/
#include "Time.H"
#include "argList.H"
#include "polyMesh.H"
#include "pointMesh.H"
#include "pointFields.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createPolyMesh.H"
{
scalarField vals0({1.0, 2.0, 3.0});
const Expression::ListConstRefWrap<scalar> wvals0(vals0);
scalarField vals1;
Expression::ListRefWrap<scalar> wvals1(vals1.size(), vals1);
wvals1 = wvals0;
return 0;
}
const pointMesh& pMesh = pointMesh::New(mesh);
// Field, dimensionedScalar as List expression
{
scalarField vals({1.0, 2.0, 3.0});
dimensionedScalar d(dimless, 4.0);
scalarField result(d.expr(vals.size())*vals.expr());
Pout<< "result:" << result << endl;
}
Info<< "Reading field p\n" << endl;
pointScalarField p
(
IOobject
(
"pointDisplacement",
runTime.timeName(),
pMesh.thisDb(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
pMesh
);
pointScalarField result
(
IOobject
(
"result",
runTime.timeName(),
pMesh.thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
pMesh,
dimensionedScalar("zero", Foam::sqr(p().dimensions()), 0)
);
// Construct value + dimensions with correct sizing
const auto oneField(dimensionedScalar(dimless, 1.0).expr(p));
// Make expression
auto expression = oneField * p.expr() + p.expr();
// Combine expressions
auto newExpression = sqr(expression);
// Assign values
result = newExpression;
DebugVar(result);
// Construct from expression
pointScalarField result2("result2", pMesh, sqrt(p.expr()));
DebugVar(result2);
return 0;
}
// ************************************************************************* //

451
applications/test/expressionTemplates/Test-expressionTemplates-volFields.C Normal file
View File

@ -0,0 +1,451 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
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/>.
Application
volPointInterpolationTest
\*---------------------------------------------------------------------------*/
#include "Time.H"
#include "argList.H"
#include "fvMesh.H"
#include "ListExpression.H"
#include "GeometricFieldExpression.H"
#include "fvCFD.H"
#include "fvMatrixExpression.H"
#include <ratio>
#include <chrono>
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
tmp<volScalarField> someFunction(const volScalarField& fld)
{
return fld*1.0001;
}
template<class Type>
void fusedGaussFvmLaplacian
(
fvMatrix<Type>& fvm,
const surfaceInterpolationScheme<scalar>& interpGammaScheme,
const fv::snGradScheme<Type>& snGradScheme,
const GeometricField<scalar, fvPatchField, volMesh>& gamma,
const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
// Replacement for gaussLaplacianScheme::fvmLaplacian with scalar gamma
typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfaceType;
const auto& mesh = vf.mesh();
// Expression for weights
const auto weights = interpGammaScheme.weights(gamma).expr();
// Expression for gamma_face * magSf
const auto gammaMagSf =
Expression::interpolate(gamma.expr(), weights, mesh)
* mesh.magSf().expr();
// Expression for deltaCoeffs
const auto deltaCoeffs = snGradScheme.deltaCoeffs(vf).expr();
// Construct matrix
Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs);
if (snGradScheme.corrected())
{
// Wrap correction
const auto corr(snGradScheme.correction(vf).expr());
const auto V = mesh.V().expr();
if (mesh.fluxRequired(vf.name()))
{
fvm.faceFluxCorrectionPtr() = std::make_unique<surfaceType>
(
IOobject
(
"faceFluxCorr",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh,
gamma.dimensions()
*mesh.magSf().dimensions()
*corr.data().dimensions()
);
auto& faceFluxCorr = *fvm.faceFluxCorrectionPtr();
faceFluxCorr = gammaMagSf*corr;
fvm.source() =
fvm.source().expr()
- (
V * fvc::div
(
faceFluxCorr
)().primitiveField().expr()
);
}
else
{
// Temporary field
surfaceType faceFluxCorr
(
IOobject
(
"faceFluxCorr",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh,
gamma.dimensions()
*mesh.magSf().dimensions()
*corr.data().dimensions()
);
faceFluxCorr = gammaMagSf*corr;
fvm.source() =
fvm.source().expr()
- (
V * fvc::div
(
faceFluxCorr
)().primitiveField().expr()
);
}
}
}
using namespace std::chrono;
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh.thisDb(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
{
volScalarField p2("p2", p);
DebugVar(p2.boundaryFieldRef());
for (auto& pf : p.boundaryFieldRef())
{
scalarField newVals(pf.size());
forAll(pf, i)
{
newVals[i] = scalar(i);
}
pf == newVals;
}
//std::vector<const scalarField*> ptrs;
// List<const scalarField*> ptrs;
// for (const auto& pp : p.boundaryField())
// {
// ptrs.push_back(&pp);
// }
DebugVar(sizeof(long));
DebugVar(p.boundaryField());
Expression::ListsConstRefWrap<scalarField> expr(p.boundaryField());
const auto twoA = expr + expr;
Expression::ListsRefWrap<scalarField> expr2(p2.boundaryFieldRef());
Pout<< "**before assignment twoA:" << twoA.size()
<< " expr2:" << expr2.size() << endl;
expr2 = twoA;
Pout<< "**after assignment twoA:" << twoA.size()
<< " expr2:" << expr2.size() << endl;
DebugVar(p2.boundaryField());
// forAll(expr, i)
// {
// Pout<< "i:" << i
// //<< " expr:" << expr[i]
// //<< " twoA:" << twoA[i]
// << " expr2:" << expr2[i]
// << endl;
// }
return 0;
}
{
//DebugVar(linearInterpolate(p));
//auto tweights = linear<scalar>(mesh).weights(p);
//DebugVar(tweights);
surfaceScalarField result
(
IOobject
(
"result",
runTime.timeName(),
mesh.thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh,
dimensionedScalar(p.dimensions(), 0)
);
//result = Expression::interpolate
//(
// p.expr(),
// tweights().expr(),
// mesh
//);
result = Expression::linearInterpolate(p.expr(), mesh);
DebugVar(result);
return 0;
}
volScalarField p2
(
IOobject
(
"p",
runTime.timeName(),
mesh.thisDb(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE,
IOobject::NO_REGISTER
),
mesh
);
// Expresions of volFields
{
volScalarField result
(
IOobject
(
"result",
runTime.timeName(),
mesh.thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh,
dimensionedScalar(p.dimensions(), 0)
);
{
Pout<< "No expression templates:" << endl;
const high_resolution_clock::time_point t1 =
high_resolution_clock::now();
result = p + p2;
const high_resolution_clock::time_point t2 =
high_resolution_clock::now();
const duration<double> time_span = t2 - t1;
Pout<< "Operation time:" << time_span.count() << endl;
}
{
Pout<< "With expression templates:" << endl;
const high_resolution_clock::time_point t1 =
high_resolution_clock::now();
result = p.expr() + p2.expr();
const high_resolution_clock::time_point t2 =
high_resolution_clock::now();
const duration<double> time_span = t2 - t1;
Pout<< "Operation time:" << time_span.count() << endl;
}
// const auto oldDimensions = p.dimensions();
// p.dimensions().reset(dimless);
// p2.dimensions().reset(dimless);
// result.dimensions().reset(dimless);
// {
//
// Pout<< "Complex expression : No expression templates:" << endl;
// const high_resolution_clock::time_point t1 =
// high_resolution_clock::now();
// result = cos(p + 0.5*sqrt(p2-sin(p)));
// const high_resolution_clock::time_point t2 =
// high_resolution_clock::now();
// const duration<double> time_span = t2 - t1;
// Pout<< "Operation time:" << time_span.count() << endl;
// }
// {
// Pout<< "Complex expression : With expression templates:" << endl;
// const high_resolution_clock::time_point t1 =
// high_resolution_clock::now();
// const auto zeroDotFive
// (
// dimensionedScalar(dimless, 0.5).expr(p)
// );
// result = cos(p.expr() + zeroDotFive*sqrt(p2.expr()-sin(p.expr())));
// const high_resolution_clock::time_point t2 =
// high_resolution_clock::now();
// const duration<double> time_span = t2 - t1;
// Pout<< "Operation time:" << time_span.count() << endl;
// }
// p.dimensions().reset(oldDimensions);
// p2.dimensions().reset(oldDimensions);
// result.dimensions().reset(oldDimensions);
return 0;
// auto expression = someFunction(p).expr() + someFunction(p).expr();
// result = expression;
// DebugVar(result);
}
{
volScalarField result
(
IOobject
(
"result",
runTime.timeName(),
mesh.thisDb(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh,
dimensionedScalar(p.dimensions(), 0)
);
auto expression = someFunction(p).expr() + someFunction(p).expr();
result = expression;
DebugVar(result);
}
// Expresions of volFields
{
volScalarField result
(
"result",
mesh,
sqr(p.expr() + p.expr())
);
DebugVar(result);
}
{
// Fill p with some values
forAll(p, celli)
{
p[celli] = celli;
}
p.correctBoundaryConditions();
// Interpolate to surface field
surfaceScalarField result
(
"result",
mesh,
Expression::interpolate //<volExpr, surfaceExpr>
(
p.expr(),
mesh.surfaceInterpolation::weights().expr(),
mesh
)
);
DebugVar(result);
}
{
// For testing as a replacement of laplacian weights
const volScalarField gamma
(
IOobject
(
"gamma",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh,
dimensionedScalar(dimless, 1.0)
);
fvMatrix<scalar> fvm
(
p,
gamma.dimensions()*mesh.magSf().dimensions()*p.dimensions()
);
const linear<scalar> interpGammaScheme(mesh);
const fv::correctedSnGrad<scalar> snGradScheme(mesh);
fusedGaussFvmLaplacian
(
fvm,
interpGammaScheme,
snGradScheme,
gamma,
p
);
DebugVar(fvm.source());
}
// Expressions of fvMatrix
{
tmp<fvMatrix<scalar>> tm0(fvm::laplacian(p));
const fvMatrix<scalar>& m0 = tm0();
DebugVar(m0.dimensions());
tmp<fvMatrix<scalar>> tm1(fvm::laplacian(p));
const fvMatrix<scalar>& m1 = tm1();
DebugVar(m1.dimensions());
fvMatrix<scalar> m2(p, m0.expr() + m1.expr());
DebugVar(m2.dimensions());
}
return 0;
}
// ************************************************************************* //

4
etc/config.sh/settings
View File

@ -81,8 +81,8 @@ Darwin)
case "$WM_COMPILE_CONTROL" in
(*version=*) ;;
(*)
WM_COMPILER=Clang
echo "openfoam (darwin): using clang instead of gcc" 1>&2
#WM_COMPILER=Clang
#echo "openfoam (darwin): using clang instead of gcc" 1>&2
;;
esac
fi

2
src/Allwmake
View File

@ -55,8 +55,6 @@ wmake $targetType fileFormats
wmake $targetType surfMesh
wmake $targetType meshTools
#------------------------------------------------------------------------------
wmake $targetType finiteArea
wmake $targetType finiteVolume
wmake $targetType fused/finiteVolume

33
src/OpenFOAM/containers/Lists/List/List.H
View File

@ -45,6 +45,7 @@ SourceFiles
#include "autoPtr.H"
#include "UList.H"
#include "SLListFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -65,7 +66,6 @@ typedef List<bool> boolList; //!< A List of bools
typedef List<char> charList; //!< A List of chars
typedef List<label> labelList; //!< A List of labels
/*---------------------------------------------------------------------------*\
Class List Declaration
\*---------------------------------------------------------------------------*/
@ -395,6 +395,37 @@ public:
//- Alias for resize()
void setSize(label n, const T& val) { this->resize(n, val); }
// Expression templates
//- Construct from value expression
template<typename E>
explicit List
(
const Expression::ListExpression
<
E,
typename E::value_type
>& expr
)
{
expr.evaluate(*this);
}
//- Assign values from expression
template<typename E>
void operator=
(
const Expression::ListExpression
<
E,
typename E::value_type
>& expr
)
{
expr.evaluate(*this);
}
};

38
src/OpenFOAM/containers/Lists/List/UList.H
View File

@ -80,6 +80,11 @@ typedef UList<bool> boolUList; //!< A UList of bools
typedef UList<char> charUList; //!< A UList of chars
typedef UList<label> labelUList; //!< A UList of labels
namespace Expression
{
template<class E, typename ValType> class ListExpression;
template<class T> class ListConstRefWrap;
}
/*---------------------------------------------------------------------------*\
Class UList Declaration
@ -697,6 +702,39 @@ public:
{
return this->contains(val, pos);
}
// Expression templates
//- Have unique tag
using is_List = void;
//- Wrap value as expression
auto expr() const
{
return Expression::ListConstRefWrap<T>
(
reinterpret_cast<const List<T>&>(*this)
);
}
//- Assign values from expression
template<typename E>
void operator=
(
const Expression::ListExpression
<
E,
typename E::value_type
>& expr
)
{
// bypass expr.evaluate to avoid resize ...
for (label i = 0; i < expr.size(); ++i)
{
operator[](i) = expr[i];
}
}
};

28
src/OpenFOAM/dimensionedTypes/dimensionedType/dimensionedType.H
View File

@ -60,6 +60,12 @@ template<class Type> class dimensioned;
template<class Type>
Istream& operator>>(Istream& is, dimensioned<Type>& dt);
namespace Expression
{
template<class T> class UniformListWrap;
template<class GeoField> class UniformGeometricFieldWrap;
}
/*---------------------------------------------------------------------------*\
Class dimensioned Declaration
\*---------------------------------------------------------------------------*/
@ -449,6 +455,28 @@ public:
{
return getOrAddToDict(name, dict, deflt);
}
// Expression templates
//- Wrap value as constant-value List expression
auto expr(const label size) const
{
return Expression::UniformListWrap<Type>(size, value_);
}
//- Wrap value as constant-value GeometricField expression. Supplied
// fields is used for sizing only and needs to be valid at time of
// evaluation.
template<class GeoField>
auto expr(const GeoField& fld) const
{
return Expression::UniformGeometricFieldWrap<GeoField>
(
fld,
*this
);
}
};

252
src/OpenFOAM/expressionTemplates/GenericExpression.H Normal file
View File

@ -0,0 +1,252 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2025 M. Janssens
-------------------------------------------------------------------------------
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/>.
Namespace
Foam::GenericExpression
Description
Expression templates.
SourceFiles
GenericExpression.H
\*---------------------------------------------------------------------------*/
#ifndef Foam_GenericExpression_H
#define Foam_GenericExpression_H
#include <cassert>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace Expression
{
/*---------------------------------------------------------------------------*\
Class GenericExpression Declarations
\*---------------------------------------------------------------------------*/
template<typename E>
class GenericExpression
{
public:
static constexpr bool is_leaf = false;
auto evaluate() const
{
// Delegation to the actual expression type. This avoids dynamic
// polymorphism (a.k.a. virtual functions in C++)
return static_cast<E const&>(*this).evaluate();
}
};
// Macros
// ~~~~~~
#define EXPRESSION_FUNCTION1(BaseClass, Func, BaseFunc, OpFunc) \
template<typename E1> \
class OpFunc \
: \
public BaseClass<OpFunc<E1>> \
{ \
typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_; \
\
public: \
static constexpr bool is_leaf = false; \
\
OpFunc(const E1& u) \
: \
u_(u) \
{} \
\
auto evaluate() const \
{ \
return BaseFunc(u_.evaluate()); \
} \
}; \
template<typename E1> \
OpFunc<E1> Func \
( \
const BaseClass<E1>& u \
) \
{ \
return OpFunc<E1>(static_cast<const E1&>(u)); \
}
#define EXPRESSION_FUNCTION2(BaseClass, Func, BaseFunc, OpFunc) \
template<typename E1, typename E2> \
class OpFunc \
: \
public BaseClass<OpFunc<E1, E2>> \
{ \
typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_; \
typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_; \
\
public: \
static constexpr bool is_leaf = false; \
\
OpFunc(const E1& u, const E2& v) \
: \
u_(u), v_(v) \
{} \
auto evaluate() const \
{ \
return BaseFunc(u_.evaluate(), v_.evaluate()); \
} \
}; \
template<typename E1, typename E2> \
OpFunc<E1, E2> Func \
( \
const BaseClass<E1>& u, \
const BaseClass<E2>& v \
) \
{ \
return OpFunc<E1, E2> \
( \
static_cast<const E1&>(u), \
static_cast<const E2&>(v) \
); \
}
#define EXPRESSION_OPERATOR(BaseClass, Op, BaseOp, OpFunc) \
template<typename E1, typename E2> \
class OpFunc \
: \
public BaseClass \
< \
OpFunc<E1, E2> \
> \
{ \
/* cref if leaf, copy otherwise */ \
typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_; \
typename std::conditional<E2::is_leaf, const E2&, const E2>::type v_; \
\
public: \
static constexpr bool is_leaf = false; \
\
OpFunc(const E1& u, const E2& v) \
: \
u_(u), v_(v) \
{} \
auto evaluate() const \
{ \
return u_.evaluate() BaseOp v_.evaluate(); \
} \
}; \
template<typename E1, typename E2> \
OpFunc<E1, E2> \
operator Op \
( \
BaseClass<E1> const& u, \
BaseClass<E2> const& v \
) \
{ \
return OpFunc<E1, E2> \
( \
static_cast<const E1&>(u), \
static_cast<const E2&>(v) \
); \
}
// Do '-' separately until we work out macro expansion...
template<typename E1>
class g_negate
:
public GenericExpression
<
g_negate<E1>
>
{
typename std::conditional<E1::is_leaf, const E1&, const E1>::type u_;
public:
static constexpr bool is_leaf = false;
g_negate(const E1& u)
:
u_(u)
{}
auto evaluate() const
{
return -u_.evaluate();
}
};
template<typename E1>
g_negate<E1> operator-
(
const GenericExpression<E1>& u
)
{
return g_negate<E1>(static_cast<const E1&>(u));
}
EXPRESSION_FUNCTION1(GenericExpression, sqr, Foam::sqr, g_sqr)
EXPRESSION_FUNCTION1(GenericExpression, sqrt, Foam::sqrt, g_sqrt)
EXPRESSION_FUNCTION1(GenericExpression, magSqr, Foam::magSqr, g_magSqr)
EXPRESSION_FUNCTION1(GenericExpression, symm, Foam::symm, g_symm)
EXPRESSION_FUNCTION1(GenericExpression, pow2, Foam::pow2, g_pow2)
EXPRESSION_FUNCTION1(GenericExpression, pow3, Foam::pow3, g_pow3)
EXPRESSION_FUNCTION1(GenericExpression, pow4, Foam::pow4, g_pow4)
EXPRESSION_FUNCTION1(GenericExpression, operator~, Foam::operator~, g_compl)
//TBD. Parse problem
//EXPRESSION_FUNCTION1(GenericExpression, operator!, Foam::operator!, g_not)
#undef EXPRESSION_FUNCTION1
EXPRESSION_FUNCTION2(GenericExpression, operator+, Foam::operator+, g_add);
EXPRESSION_FUNCTION2(GenericExpression, operator-, Foam::operator-, g_subtract);
EXPRESSION_FUNCTION2(GenericExpression, operator*, Foam::operator*, g_multiply);
EXPRESSION_FUNCTION2(GenericExpression, operator/, Foam::operator/, g_divide);
#undef EXPRESSION_FUNCTION2
EXPRESSION_OPERATOR(GenericExpression, ||, ||, g_or)
EXPRESSION_OPERATOR(GenericExpression, &&, &&, g_and)
EXPRESSION_OPERATOR(GenericExpression, &, &, g_bitand)
EXPRESSION_OPERATOR(GenericExpression, |, |, g_bitor)
EXPRESSION_OPERATOR(GenericExpression, ^, ^, g_xor)
#undef EXPRESSION_OPERATOR
} // End namespace Expression
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

1879
src/OpenFOAM/expressionTemplates/GeometricFieldExpression.H Normal file
View File

File diff suppressed because it is too large Load Diff

2411
src/OpenFOAM/expressionTemplates/ListExpression.H Normal file
View File

File diff suppressed because it is too large Load Diff

145
src/OpenFOAM/expressionTemplates/boolExpression.H Normal file
View File

@ -0,0 +1,145 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2025 M. Janssens
-------------------------------------------------------------------------------
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/>.
Namespace
Foam::boolExpression
Description
Expression templates for bool
SourceFiles
boolExpression.H
\*---------------------------------------------------------------------------*/
#ifndef Foam_boolExpression_H
#define Foam_boolExpression_H
#include "GenericExpression.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace Expression
{
/*---------------------------------------------------------------------------*\
Namespace Expression Declarations
\*---------------------------------------------------------------------------*/
// Wrap of non-const bool
// ~~~~~~~~~~~~~~~~~~~~~~
class boolWrap
:
public GenericExpression<boolWrap>
{
bool elems_;
public:
static constexpr bool is_leaf = false; //true;
//- Construct from components
boolWrap(const bool elems)
:
elems_(elems)
{}
// Construct from GenericExpression, forcing its evaluation.
template<typename E>
boolWrap
(
const GenericExpression<E>& expr
)
:
elems_(expr.evaluate())
{}
auto evaluate() const
{
return elems_;
}
template<typename E>
auto evaluate
(
const GenericExpression<E>& expr
)
{
elems_ = expr.evaluate();
return elems_;
}
//- Assignment
template<typename E>
void operator=
(
const GenericExpression<E>& expr
)
{
elems_ = expr.evaluate();
}
};
// Wrap of const bool
// ~~~~~~~~~~~~~~~~~~
class boolConstWrap
:
public GenericExpression<boolConstWrap>
{
const bool elems_;
public:
// ! Store as copy (since holds reference)
static constexpr bool is_leaf = false;
// construct from components
boolConstWrap(const bool elems)
:
elems_(elems)
{}
auto evaluate() const
{
return elems_;
}
};
} // End namespace Expression
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

148
src/OpenFOAM/expressionTemplates/dimensionSetExpression.H Normal file
View File

@ -0,0 +1,148 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2025 M. Janssens
-------------------------------------------------------------------------------
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/>.
Namespace
Foam::dimensionSetExpression
Description
Expression templates for dimensionSet
SourceFiles
dimensionSetExpression.H
\*---------------------------------------------------------------------------*/
#ifndef Foam_dimensionSetExpression_H
#define Foam_dimensionSetExpression_H
#include "GenericExpression.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace Expression
{
/*---------------------------------------------------------------------------*\
Namespace Expression Declarations
\*---------------------------------------------------------------------------*/
// Wrap of non-const reference to dimensionSet
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class dimensionSetRefWrap
:
public GenericExpression<dimensionSetRefWrap>
{
dimensionSet& elems_;
public:
static constexpr bool is_leaf = false; //true;
//- Construct from components
dimensionSetRefWrap(dimensionSet& elems)
:
elems_(elems)
{}
// Construct from GenericExpression, forcing its evaluation.
template<typename E>
dimensionSetRefWrap
(
dimensionSet& elems,
const GenericExpression<E>& expr
)
:
elems_(elems)
{
elems_ = expr.evaluate();
}
auto evaluate() const
{
return elems_;
}
template<typename E>
auto& evaluate
(
const GenericExpression<E>& expr
)
{
elems_ = expr.evaluate();
return elems_;
}
//- Assignment
template<typename E>
void operator=
(
const GenericExpression<E>& expr
)
{
elems_ = expr.evaluate();
}
};
// Wrap of const reference to dimensionSet
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class dimensionSetConstRefWrap
:
public GenericExpression<dimensionSetConstRefWrap>
{
const dimensionSet& elems_;
public:
// ! Store as copy (since holds reference)
static constexpr bool is_leaf = false;
// construct from components
dimensionSetConstRefWrap(const dimensionSet& elems)
:
elems_(elems)
{}
auto evaluate() const
{
return elems_;
}
};
} // End namespace Expression
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

31
src/OpenFOAM/fields/Fields/Field/Field.H
View File

@ -582,6 +582,37 @@ public:
friend Ostream& operator<< <Type>
(Ostream&, const tmp<Field<Type>>&);
// Expression templates
//- Construct from value expression
template<typename E>
explicit Field
(
const Expression::ListExpression
<
E,
typename E::value_type
>& expr
)
{
expr.evaluate(*this);
}
//- Assign values from expression
template<typename E>
void operator=
(
const Expression::ListExpression
<
E,
typename E::value_type
>& expr
)
{
expr.evaluate(*this);
}
};

118
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C
View File

@ -1675,6 +1675,124 @@ Foam::Ostream& Foam::operator<<
}
// * * * * * * * * * * * * * Expression Templates * * * * * * * * * * * * * //
template<class Type, template<class> class PatchField, class GeoMesh>
template<typename E>
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
(
const IOobject& io,
const Mesh& mesh,
const Expression::GeometricFieldExpression
<
E,
typename E::IntExpr,
typename E::UncoupledPatchExpr,
typename E::CoupledPatchExpr,
typename E::value_type
>& expr
)
:
Internal(io, mesh, expr.dimensions(), false),
timeIndex_(this->time().timeIndex()),
boundaryField_(mesh.boundary(), *this, PatchField<Type>::calculatedType())
{
DebugInFunction
<< "Creating from expression " << nl << this->info() << endl;
bool hasRead = readIfPresent();
if (!hasRead)
{
expr.evaluate(*this);
}
}
template<class Type, template<class> class PatchField, class GeoMesh>
template<typename E>
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
(
const word& name,
const Mesh& mesh,
const Expression::GeometricFieldExpression
<
E,
typename E::IntExpr,
typename E::UncoupledPatchExpr,
typename E::CoupledPatchExpr,
typename E::value_type
>& expr
)
:
Internal
(
IOobject
(
name,
mesh.time().timeName(),
mesh.thisDb()
),
mesh,
expr.dimensions(),
false
),
timeIndex_(this->time().timeIndex()),
boundaryField_(mesh.boundary(), *this, PatchField<Type>::calculatedType())
{
DebugInFunction
<< "Creating from expression " << nl << this->info() << endl;
expr.evaluate(*this);
}
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::Expression::GeometricFieldConstRefWrap
<
Foam::GeometricField<Type, PatchField, GeoMesh>
>
Foam::GeometricField<Type, PatchField, GeoMesh>::expr() const
{
return Expression::GeometricFieldConstRefWrap<this_type>(*this);
}
template<class Type, template<class> class PatchField, class GeoMesh>
template<typename E>
void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
(
const Expression::GeometricFieldExpression
<
E,
typename E::IntExpr,
typename E::UncoupledPatchExpr,
typename E::CoupledPatchExpr,
typename E::value_type
>& expr
)
{
Expression::GeometricFieldRefWrap<this_type>(*this, expr);
}
template<class Type, template<class> class PatchField, class GeoMesh>
template<typename E>
void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
(
const Expression::GeometricFieldExpression
<
E,
typename E::IntExpr,
typename E::UncoupledPatchExpr,
typename E::CoupledPatchExpr,
typename E::value_type
>& expr
)
{
expr.evaluate(*this, true);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#undef checkField

74
src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H
View File

@ -42,6 +42,7 @@ SourceFiles
#define Foam_GeometricField_H
#include "GeometricBoundaryField.H"
#include "GeometricFieldExpression.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -1142,6 +1143,79 @@ public:
{
return this->clamp_range(lo.value(), hi.value());
}
// Expression templates
//- Have unique tag
using is_GeometricField = void;
//- Construct from IOobject and value expression. Supports
// read_if_present
template<typename E>
GeometricField
(
const IOobject& io,
const Mesh& mesh,
const Expression::GeometricFieldExpression
<
E,
typename E::IntExpr,
typename E::UncoupledPatchExpr,
typename E::CoupledPatchExpr,
typename E::value_type
>& expr
);
//- Construct from name, mesh and value expression.
template<typename E>
GeometricField
(
const word& name,
const Mesh& mesh,
const Expression::GeometricFieldExpression
<
E,
typename E::IntExpr,
typename E::UncoupledPatchExpr,
typename E::CoupledPatchExpr,
typename E::value_type
>& expr
);
//- Wrap value as expression
typename Expression::GeometricFieldConstRefWrap
<
this_type
> expr() const;
//- Assign values from expression
template<typename E>
void operator=
(
const Expression::GeometricFieldExpression
<
E,
typename E::IntExpr,
typename E::UncoupledPatchExpr,
typename E::CoupledPatchExpr,
typename E::value_type
>& expr
);
//- Assign values from expression. Override asssignable.
template<typename E>
void operator==
(
const Expression::GeometricFieldExpression
<
E,
typename E::IntExpr,
typename E::UncoupledPatchExpr,
typename E::CoupledPatchExpr,
typename E::value_type
>& expr
);
};

36
src/OpenFOAM/memory/tmp/tmp.H
View File

@ -59,8 +59,17 @@ See also
namespace Foam
{
namespace Expression
{
template<class GeoField> class GeometricFieldConstTmpWrap;
template<class T> class ListConstTmpWrap;
}
// Forward Declarations
template<class T> class refPtr;
template<class T> class tmp;
template<class Type, template<class> class PatchField, class GeoMesh>
class GeometricField;
/*---------------------------------------------------------------------------*\
Class tmp Declaration
@ -352,6 +361,33 @@ public:
// \deprecated(2020-07) - use bool operator
FOAM_DEPRECATED_FOR(2020-07, "bool operator")
bool empty() const noexcept { return !ptr_; }
// Expression templates
//- Wrap value as expression
auto expr() const
{
// Add 'using is_GeometricField = void' to GeometricField
// Add 'using is_List = void' to List
if constexpr
(
std::is_void_v<typename T::is_GeometricField>
)
{
return Expression::GeometricFieldConstTmpWrap<T>(*this);
}
else //if constexpr (std::is_void_v<T::is_List>)
{
return Expression::ListConstTmpWrap<T>(*this);
}
//else
//{
// //static_assert("XXX",YYY);
// return;
//}
}
};

4
src/OpenFOAM/memory/tmp/tmpI.H
View File

@ -131,7 +131,7 @@ inline Foam::tmp<T>::tmp(const tmp<T>& rhs)
if (ptr_)
{
ptr_->refCount::operator++();
this->checkUseCount();
//this->checkUseCount();
}
else
{
@ -162,7 +162,7 @@ inline Foam::tmp<T>::tmp(const tmp<T>& rhs, bool reuse)
else
{
ptr_->refCount::operator++();
this->checkUseCount();
//this->checkUseCount();
}
}
else

1747
src/finiteVolume/expressionTemplates/fvMatrixExpression.H Normal file
View File

File diff suppressed because it is too large Load Diff

71
src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
View File

@ -2984,6 +2984,77 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const fvMatrix<Type>& fvm)
}
// * * * * * * * * * * * * * Expression Templates * * * * * * * * * * * * * //
template<class Type>
template<typename E>
Foam::fvMatrix<Type>::fvMatrix
(
const GeometricField<Type, fvPatchField, volMesh>& psi,
const Expression::fvMatrixExpression
<
E,
typename E::DiagExpr,
typename E::UpperExpr,
typename E::LowerExpr,
typename E::FaceFluxExpr,
typename E::SourceExpr
>& expr
)
:
lduMatrix(psi.mesh()),
psi_(psi),
useImplicit_(false),
lduAssemblyName_(),
nMatrix_(0),
dimensions_(expr.dimensions()),
source_(psi.size(), Zero),
internalCoeffs_(psi.mesh().boundary().size()),
boundaryCoeffs_(psi.mesh().boundary().size())
{
DebugInFunction
<< "Constructing fvMatrix<Type> from expression for field "
<< psi_.name() << endl;
checkImplicit();
// Fill in diag,upper,lower etc.
expr.evaluate(*this);
auto& psiRef = this->psi(0);
const label currentStatePsi = psiRef.eventNo();
psiRef.boundaryFieldRef().updateCoeffs();
psiRef.eventNo() = currentStatePsi;
}
template<class Type>
Foam::Expression::fvMatrixConstRefWrap<Foam::fvMatrix<Type>>
Foam::fvMatrix<Type>::expr() const
{
return Expression::fvMatrixConstRefWrap<Foam::fvMatrix<Type>>(*this);
}
template<class Type>
template<typename E>
void Foam::fvMatrix<Type>::operator=
(
const Expression::fvMatrixExpression
<
E,
typename E::DiagExpr,
typename E::UpperExpr,
typename E::LowerExpr,
typename E::FaceFluxExpr,
typename E::SourceExpr
>& expr
)
{
expr.evaluate(*this);
}
// * * * * * * * * * * * * * * * * Solvers * * * * * * * * * * * * * * * * * //
#include "fvMatrixSolve.C"

43
src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H
View File

@ -53,6 +53,8 @@ SourceFiles
#include "lduPrimitiveMeshAssembly.H"
#include "lduMesh.H"
#include "fvMatrixExpression.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
@ -741,6 +743,47 @@ public:
Ostream&,
const fvMatrix<Type>&
);
// Expression templates
//- Construct given a field to solve for and expressions for
//- all the components
template<typename E>
fvMatrix
(
const GeometricField<Type, fvPatchField, volMesh>& psi,
const Expression::fvMatrixExpression
<
E,
typename E::DiagExpr,
typename E::UpperExpr,
typename E::LowerExpr,
typename E::FaceFluxExpr,
typename E::SourceExpr
>& expr
);
//- Wrap value as expression
typename Expression::fvMatrixConstRefWrap
<
fvMatrix<Type>
> expr() const;
//- Assign values from expression
template<typename E>
void operator=
(
const Expression::fvMatrixExpression
<
E,
typename E::DiagExpr,
typename E::UpperExpr,
typename E::LowerExpr,
typename E::FaceFluxExpr,
typename E::SourceExpr
>& expr
);
};

81
src/fused/finiteVolume/fusedGaussLaplacianScheme.C
View File

@ -45,6 +45,87 @@ namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type, class GType>
template<class E1, class E2>
void fusedGaussLaplacianScheme<Type, GType>::fvmCorrection
(
fvMatrix<Type>& fvm,
const dimensionSet& gammaDim,
const E1& gammaMagSf,
const E2& corr
)
{
typedef GeometricField<Type, fvsPatchField, surfaceMesh> surfaceType;
const auto& vf = fvm.psi();
const fvMesh& mesh = vf.mesh();
const auto V = mesh.V().expr();
if (mesh.fluxRequired(vf.name()))
{
fvm.faceFluxCorrectionPtr() = std::make_unique<surfaceType>
(
IOobject
(
"faceFluxCorr",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh,
gammaDim
*mesh.magSf().dimensions()
*corr.data().dimensions()
);
auto& faceFluxCorr = *fvm.faceFluxCorrectionPtr();
faceFluxCorr = corr*gammaMagSf;
fvm.source() =
fvm.source().expr()
- (
fvc::div
(
faceFluxCorr
)().primitiveField().expr()
*V
);
}
else
{
surfaceType faceFluxCorr
(
IOobject
(
"faceFluxCorr",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
mesh,
gammaDim
*mesh.magSf().dimensions()
*corr.data().dimensions()
);
faceFluxCorr = corr*gammaMagSf;
fvm.source() =
fvm.source().expr()
- (
fvc::div
(
faceFluxCorr
)().primitiveField().expr()
*V
);
}
}
template<class Type, class GType>
tmp<fvMatrix<Type>>
fusedGaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected

27
src/fused/finiteVolume/fusedGaussLaplacianScheme.H
View File

@ -65,6 +65,17 @@ class fusedGaussLaplacianScheme
{
// Private Member Functions
//- Helper for calculate source correction. Partially expression
// templated. Move to fvMatrixExpression.H?
template<class E1, class E2>
static void fvmCorrection
(
fvMatrix<Type>& fvm,
const dimensionSet& gammaDim,
const E1& gammaMagSf,
const E2& corr
);
//- See gaussLaplacianScheme
tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> gammaSnGradCorr
(
@ -198,6 +209,8 @@ defineFvmLaplacianScalarGamma(sphericalTensor);
defineFvmLaplacianScalarGamma(symmTensor);
defineFvmLaplacianScalarGamma(tensor);
#undef defineFvmLaplacianScalarGamma
// Unweighted laplacian
template<>
tmp<GeometricField<scalar, fvPatchField, volMesh>>
@ -227,20 +240,6 @@ fusedGaussLaplacianScheme<vector, scalar>::fvcLaplacian
const GeometricField<scalar, fvPatchField, volMesh>&,
const GeometricField<vector, fvPatchField, volMesh>&
);
template<>
tmp<fvMatrix<scalar>>
fusedGaussLaplacianScheme<scalar, scalar>::fvmLaplacian
(
const GeometricField<scalar, fvPatchField, volMesh>&,
const GeometricField<scalar, fvPatchField, volMesh>&
);
template<>
tmp<fvMatrix<vector>>
fusedGaussLaplacianScheme<vector, scalar>::fvmLaplacian
(
const GeometricField<scalar, fvPatchField, volMesh>&,
const GeometricField<vector, fvPatchField, volMesh>&
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

125
src/fused/finiteVolume/fusedGaussLaplacianSchemes.C
View File

@ -28,6 +28,7 @@ License
#include "fusedGaussLaplacianScheme.H"
#include "fvMesh.H"
#include "fvMatrixExpression.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -40,61 +41,81 @@ Foam::tmp<Foam::fvMatrix<Foam::Type>> \
Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>:: \
fvmLaplacian \
( \
const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma, \
const GeometricField<scalar, fvPatchField, volMesh>& gamma, \
const GeometricField<Type, fvPatchField, volMesh>& vf \
) \
{ \
DebugPout<< "fusedGaussLaplacianScheme::fvmLaplacian on " << vf.name() \
<< " with scalar gamma " << gamma.name() << endl; \
<< " with scalar gamma " << gamma.name() << " and ET" << endl; \
\
const fvMesh& mesh = this->mesh(); \
\
GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf \
( \
gamma*mesh.magSf() \
); \
const auto weights = this->tinterpGammaScheme_().weights(gamma).expr(); \
const auto gammaMagSf = \
Expression::interpolate(gamma.expr(), weights, mesh) \
* mesh.magSf().expr(); \
/* For compatibility with linearInterpolate : avoid orientation. TBD. */ \
const_cast<orientedType&>(gammaMagSf.oriented()).setOriented(false); \
const auto deltaCoeffs = this->tsnGradScheme_().deltaCoeffs(vf).expr(); \
\
tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected \
tmp<fvMatrix<Type>> tfvm \
( \
gammaMagSf, \
this->tsnGradScheme_().deltaCoeffs(vf), \
vf \
new fvMatrix<Type> \
( \
vf, \
gamma.dimensions()*mesh.magSf().dimensions()*vf.dimensions() \
) \
); \
fvMatrix<Type>& fvm = tfvm.ref(); \
\
Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs); \
\
if (this->tsnGradScheme_().corrected()) \
{ \
if (mesh.fluxRequired(vf.name())) \
{ \
fvm.faceFluxCorrectionPtr() = std::make_unique \
< \
GeometricField<Type, fvsPatchField, surfaceMesh> \
> \
( \
gammaMagSf*this->tsnGradScheme_().correction(vf) \
); \
\
fvm.source() -= \
mesh.V()* \
fvc::div \
( \
*fvm.faceFluxCorrectionPtr() \
)().primitiveField(); \
} \
else \
{ \
fvm.source() -= \
mesh.V()* \
fvc::div \
( \
gammaMagSf*this->tsnGradScheme_().correction(vf) \
)().primitiveField(); \
} \
const auto corr(this->tsnGradScheme_().correction(vf).expr()); \
fvmCorrection(fvm, gamma.dimensions(), gammaMagSf, corr); \
} \
\
return tfvm; \
} \
\
template<> \
Foam::tmp<Foam::fvMatrix<Foam::Type>> \
Foam::fv::fusedGaussLaplacianScheme<Foam::Type, Foam::scalar>:: \
fvmLaplacian \
( \
const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma, \
const GeometricField<Type, fvPatchField, volMesh>& vf \
) \
{ \
DebugPout<< "fusedGaussLaplacianScheme::fvmLaplacian on " << vf.name() \
<< " with interpolated gamma " << gamma.name() << " and ET" << endl; \
\
const fvMesh& mesh = this->mesh(); \
\
const auto gammaMagSf = gamma.expr()* mesh.magSf().expr(); \
const auto deltaCoeffs = this->tsnGradScheme_().deltaCoeffs(vf).expr(); \
\
tmp<fvMatrix<Type>> tfvm \
( \
new fvMatrix<Type> \
( \
vf, \
gamma.dimensions()*mesh.magSf().dimensions()*vf.dimensions() \
) \
); \
fvMatrix<Type>& fvm = tfvm.ref(); \
\
Expression::fvmLaplacianUncorrected(fvm, gammaMagSf, deltaCoeffs); \
\
if (this->tsnGradScheme_().corrected()) \
{ \
const auto corr(this->tsnGradScheme_().correction(vf).expr()); \
fvmCorrection(fvm, gamma.dimensions(), gammaMagSf, corr); \
} \
\
return tfvm; \
} \
\
template<> \
Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh>> \
@ -500,38 +521,6 @@ Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvcLaplacian
}
template<>
Foam::tmp<Foam::fvMatrix<Foam::scalar>>
Foam::fv::fusedGaussLaplacianScheme<Foam::scalar, Foam::scalar>::fvmLaplacian
(
const GeometricField<scalar, fvPatchField, volMesh>& gamma,
const GeometricField<scalar, fvPatchField, volMesh>& vf
)
{
// TBD
DebugPout
<< "fusedGaussLaplacianScheme<scalar, scalar>::fvmLaplacian"
<< " on " << vf.name() << " with gamma " << gamma.name() << endl;
return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
}
template<>
Foam::tmp<Foam::fvMatrix<Foam::vector>>
Foam::fv::fusedGaussLaplacianScheme<Foam::vector, Foam::scalar>::fvmLaplacian
(
const GeometricField<scalar, fvPatchField, volMesh>& gamma,
const GeometricField<vector, fvPatchField, volMesh>& vf
)
{
// TBD
DebugPout
<< "fusedGaussLaplacianScheme<vector, scalar>::fvmLaplacian"
<< " on " << vf.name() << " with gamma " << gamma.name() << endl;
return fvmLaplacian(this->tinterpGammaScheme_().interpolate(gamma)(), vf);
}
template<>
Foam::tmp
<