zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

interpolation: Added field evaluation

Browse Source 
The field evaluations have been optimised using an additional
fieldInterpolation base class so that a virtual call happens only once
per field. This is the same pattern as that used to optimise Function1.
This commit is contained in:
Will Bainbridge
2020-09-16 16:15:35 +01:00
parent ea777f806b
commit 1d9f39761e
17 changed files with 225 additions and 136 deletions
Download Patch File Download Diff File Expand all files Collapse all files

123
src/finiteVolume/interpolation/interpolation/interpolation/interpolation.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,16 +37,123 @@ Foam::interpolation<Type>::interpolation
)
:
psi_(psi),
pMesh_(psi.mesh()),
pMeshPoints_(pMesh_.points()),
pMeshFaces_(pMesh_.faces()),
pMeshFaceCentres_(pMesh_.faceCentres()),
pMeshFaceAreas_(pMesh_.faceAreas())
mesh_(psi.mesh())
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
template<class Type>
Foam::autoPtr<Foam::interpolation<Type>> Foam::interpolation<Type>::New
(
const word& interpolationType,
const GeometricField<Type, fvPatchField, volMesh>& psi
)
{
typename dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(interpolationType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorInFunction
<< "Unknown interpolation type " << interpolationType
<< " for field " << psi.name() << nl << nl
<< "Valid interpolation types : " << endl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<interpolation<Type>>(cstrIter()(psi));
}
template<class Type>
Foam::autoPtr<Foam::interpolation<Type>> Foam::interpolation<Type>::New
(
const dictionary& interpolationSchemes,
const GeometricField<Type, fvPatchField, volMesh>& psi
)
{
return New(word(interpolationSchemes.lookup(psi.name())), psi);
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
template<class Type>
Type Foam::interpolation<Type>::interpolate
(
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei
) const
{
return
interpolate
(
tetIs.tet(mesh_).barycentricToPoint(coordinates),
tetIs.cell(),
facei
);
}
template<class Type, class InterpolationType>
Foam::tmp<Foam::Field<Type>>
Foam::fieldInterpolation<Type, InterpolationType>::interpolate
(
const vectorField& position,
const labelField& celli,
const labelField& facei
) const
{
tmp<Field<Type>> tField(new Field<Type>(position.size()));
Field<Type>& field = tField.ref();
forAll(field, i)
{
field[i] =
static_cast<const InterpolationType&>(*this).interpolate
(
position[i],
celli[i],
isNull(facei) ? -1 : facei[i]
);
}
return tField;
}
template<class Type, class InterpolationType>
Foam::tmp<Foam::Field<Type>>
Foam::fieldInterpolation<Type, InterpolationType>::interpolate
(
const Field<barycentric>& coordinates,
const labelField& celli,
const labelField& tetFacei,
const labelField& tetPti,
const labelField& facei
) const
{
tmp<Field<Type>> tField(new Field<Type>(coordinates.size()));
Field<Type>& field = tField.ref();
forAll(field, i)
{
field[i] =
static_cast<const InterpolationType&>(*this).interpolate
(
coordinates[i],
tetIndices(celli[i], tetFacei[i], tetPti[i]),
isNull(facei) ? -1 : facei[i]
);
}
return tField;
}
#include "interpolationNew.C"
// ************************************************************************* //

81
src/finiteVolume/interpolation/interpolation/interpolation/interpolation.H
View File

@ -54,18 +54,15 @@ class polyMesh;
template<class Type>
class interpolation
{
protected:
// Protected data
//- The vol field to interpolate
const GeometricField<Type, fvPatchField, volMesh>& psi_;
const polyMesh& pMesh_;
const vectorField& pMeshPoints_;
const faceList& pMeshFaces_;
const vectorField& pMeshFaceCentres_;
const vectorField& pMeshFaceAreas_;
//- Reference to the mesh
const polyMesh& mesh_;
public:
@ -135,25 +132,73 @@ public:
const label facei = -1
) const = 0;
//- As above, but for a field
virtual tmp<Field<Type>> interpolate
(
const vectorField& position,
const labelField& celli,
const labelField& facei = NullObjectRef<labelField>()
) const = 0;
//- Interpolate field to the given coordinates in the tetrahedron
// defined by the given indices. Calls interpolate function
// above here execpt where overridden by derived
// defined by the given indices. Calls interpolate function
// above here except where overridden by derived
// interpolation types.
virtual Type interpolate
(
const barycentric& coordinates,
const tetIndices& tetIs,
const label facei = -1
) const
{
return
interpolate
(
tetIs.tet(pMesh_).barycentricToPoint(coordinates),
tetIs.cell(),
facei
);
}
) const;
//- As above, but for a field
virtual tmp<Field<Type>> interpolate
(
const Field<barycentric>& coordinates,
const labelField& celli,
const labelField& tetFacei,
const labelField& tetPti,
const labelField& facei = NullObjectRef<labelField>()
) const = 0;
};
/*---------------------------------------------------------------------------*\
Class fieldInterpolation Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class InterpolationType>
class fieldInterpolation
:
public interpolation<Type>
{
public:
// Constructors
//- Inherit base class constructors
using interpolation<Type>::interpolation;
// Member Functions
//- Interpolate field to the given points in the given cells
virtual tmp<Field<Type>> interpolate
(
const vectorField& position,
const labelField& celli,
const labelField& facei = NullObjectRef<labelField>()
) const;
//- Interpolate field to the given coordinates in the given tetrahedra
virtual tmp<Field<Type>> interpolate
(
const Field<barycentric>& coordinates,
const labelField& celli,
const labelField& tetFacei,
const labelField& tetPti,
const labelField& facei = NullObjectRef<labelField>()
) const;
};

66
src/finiteVolume/interpolation/interpolation/interpolation/interpolationNew.C
View File

@ -1,66 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "interpolation.H"
#include "volFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::autoPtr<Foam::interpolation<Type>> Foam::interpolation<Type>::New
(
const word& interpolationType,
const GeometricField<Type, fvPatchField, volMesh>& psi
)
{
typename dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(interpolationType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorInFunction
<< "Unknown interpolation type " << interpolationType
<< " for field " << psi.name() << nl << nl
<< "Valid interpolation types : " << endl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<interpolation<Type>>(cstrIter()(psi));
}
template<class Type>
Foam::autoPtr<Foam::interpolation<Type>> Foam::interpolation<Type>::New
(
const dictionary& interpolationSchemes,
const GeometricField<Type, fvPatchField, volMesh>& psi
)
{
return New(word(interpolationSchemes.lookup(psi.name())), psi);
}
// ************************************************************************* //

4
src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,7 +34,7 @@ Foam::interpolationCell<Type>::interpolationCell
const GeometricField<Type, fvPatchField, volMesh>& psi
)
:
interpolation<Type>(psi)
fieldInterpolation<Type, interpolationCell<Type>>(psi)
{}

2
src/finiteVolume/interpolation/interpolation/interpolationCell/interpolationCell.H
View File

@ -48,7 +48,7 @@ class fvMesh;
template<class Type>
class interpolationCell
:
public interpolation<Type>
public fieldInterpolation<Type, interpolationCell<Type>>
{
public:

4
src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/interpolationCellPatchConstrained.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,7 +34,7 @@ Foam::interpolationCellPatchConstrained<Type>::interpolationCellPatchConstrained
const GeometricField<Type, fvPatchField, volMesh>& psi
)
:
interpolation<Type>(psi)
fieldInterpolation<Type, interpolationCellPatchConstrained<Type>>(psi)
{}

2
src/finiteVolume/interpolation/interpolation/interpolationCellPatchConstrained/interpolationCellPatchConstrained.H
View File

@ -50,7 +50,7 @@ class fvMesh;
template<class Type>
class interpolationCellPatchConstrained
:
public interpolation<Type>
public fieldInterpolation<Type, interpolationCellPatchConstrained<Type>>
{
public:

6
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,7 +34,7 @@ Foam::interpolationCellPoint<Type>::interpolationCellPoint
const GeometricField<Type, fvPatchField, volMesh>& psi
)
:
interpolation<Type>(psi),
fieldInterpolation<Type, interpolationCellPoint<Type>>(psi),
psip_
(
volPointInterpolation::New(psi.mesh()).interpolate
@ -57,7 +57,7 @@ Foam::interpolationCellPoint<Type>::interpolationCellPoint
tmp<GeometricField<Type, pointPatchField, pointMesh>> psip
)
:
interpolation<Type>(psi),
fieldInterpolation<Type, interpolationCellPoint<Type>>(psi),
psip_(psip)
{
// Uses cellPointWeight to do interpolation which needs tet decomposition

4
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -48,7 +48,7 @@ namespace Foam
template<class Type>
class interpolationCellPoint
:
public interpolation<Type>
public fieldInterpolation<Type, interpolationCellPoint<Type>>
{
protected:

6
src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPointI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -51,7 +51,7 @@ inline Type Foam::interpolationCellPoint<Type>::interpolate
const label facei
) const
{
return interpolate(cellPointWeight(this->pMesh_, position, celli, facei));
return interpolate(cellPointWeight(this->mesh_, position, celli, facei));
}
@ -79,7 +79,7 @@ inline Type Foam::interpolationCellPoint<Type>::interpolate
}
}
const triFace triIs = tetIs.faceTriIs(this->pMesh_);
const triFace triIs = tetIs.faceTriIs(this->mesh_);
return
this->psi_[tetIs.cell()]*coordinates[0]

10
src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTet.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,8 +53,8 @@ bool interpolationCellPointFace<Type>::findTet
{
bool foundTet = false;
const labelList& thisFacePoints = this->pMeshFaces_[nFace];
tetPoints[2] = this->pMeshFaceCentres_[nFace];
const labelList& thisFacePoints = this->mesh_.faces()[nFace];
tetPoints[2] = this->mesh_.faceCentres()[nFace];
label pointi = 0;
@ -65,8 +65,8 @@ bool interpolationCellPointFace<Type>::findTet
tetPointLabels[0] = thisFacePoints[pointi];
tetPointLabels[1] = thisFacePoints[nextPointLabel];
tetPoints[0] = this->pMeshPoints_[tetPointLabels[0]];
tetPoints[1] = this->pMeshPoints_[tetPointLabels[1]];
tetPoints[0] = this->mesh_.points()[tetPointLabels[0]];
tetPoints[1] = this->mesh_.points()[tetPointLabels[1]];
bool inside = true;
scalar dist = 0.0;

8
src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTriangle.H
View File

@ -21,8 +21,8 @@ bool interpolationCellPointFace<Type>::findTriangle
{
bool foundTriangle = false;
vector tetPoints[3];
const labelList& facePoints = this->pMeshFaces_[nFace];
tetPoints[2] = this->pMeshFaceCentres_[nFace];
const labelList& facePoints = this->mesh_.faces()[nFace];
tetPoints[2] = this->mesh_.faceCentres()[nFace];
label pointi = 0;
@ -35,8 +35,8 @@ bool interpolationCellPointFace<Type>::findTriangle
tetPointLabels[0] = facePoints[pointi];
tetPointLabels[1] = facePoints[nextPointLabel];
tetPoints[0] = this->pMeshPoints_[tetPointLabels[0]];
tetPoints[1] = this->pMeshPoints_[tetPointLabels[1]];
tetPoints[0] = this->mesh_.points()[tetPointLabels[0]];
tetPoints[1] = this->mesh_.points()[tetPointLabels[1]];
vector fc = (tetPoints[0] + tetPoints[1] + tetPoints[2])/3.0;

30
src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -39,7 +39,7 @@ Foam::interpolationCellPointFace<Type>::interpolationCellPointFace
const GeometricField<Type, fvPatchField, volMesh>& psi
)
:
interpolation<Type>(psi),
fieldInterpolation<Type, interpolationCellPointFace<Type>>(psi),
psip_
(
volPointInterpolation::New(psi.mesh()).interpolate
@ -73,8 +73,8 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
// only use face information when the position is on a face
if (facei < 0)
{
const vector& cellCentre = this->pMesh_.cellCentres()[celli];
const labelList& cellFaces = this->pMesh_.cells()[celli];
const vector& cellCentre = this->mesh_.cellCentres()[celli];
const labelList& cellFaces = this->mesh_.cells()[celli];
vector projection = position - cellCentre;
tetPoints[3] = cellCentre;
@ -92,10 +92,10 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
{
label nFace = cellFaces[facei];
vector normal = this->pMeshFaceAreas_[nFace];
vector normal = this->mesh_.faceAreas()[nFace];
normal /= mag(normal);
const vector& faceCentreTmp = this->pMeshFaceCentres_[nFace];
const vector& faceCentreTmp = this->mesh_.faceCentres()[nFace];
scalar multiplierNumerator = (faceCentreTmp - cellCentre) & normal;
scalar multiplierDenominator = projection & normal;
@ -167,7 +167,7 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
while (facei < cellFaces.size() && !foundTet)
{
label nFace = cellFaces[facei];
if (nFace < this->pMeshFaceAreas_.size())
if (nFace < this->mesh_.faceAreas().size())
{
foundTet = findTet
(
@ -221,7 +221,7 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
else
{
label patchi =
this->pMesh_.boundaryMesh().whichPatch(closestFace);
this->mesh_.boundaryMesh().whichPatch(closestFace);
// If the boundary patch is not empty use the face value
// else use the cell value
@ -229,7 +229,7 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
{
ts[2] = this->psi_.boundaryField()[patchi]
[
this->pMesh_.boundaryMesh()[patchi].whichFace
this->mesh_.boundaryMesh()[patchi].whichFace
(
closestFace
)
@ -265,7 +265,7 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
else
{
label patchi =
this->pMesh_.boundaryMesh().whichPatch(closestFace);
this->mesh_.boundaryMesh().whichPatch(closestFace);
// If the boundary patch is not empty use the face value
// else use the cell value
@ -273,7 +273,7 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
{
t = this->psi_.boundaryField()[patchi]
[
this->pMesh_.boundaryMesh()[patchi].whichFace
this->mesh_.boundaryMesh()[patchi].whichFace
(
closestFace
)
@ -312,14 +312,14 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
}
else
{
label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
label patchi = this->mesh_.boundaryMesh().whichPatch(facei);
// If the boundary patch is not empty use the face value
// else use the cell value
if (this->psi_.boundaryField()[patchi].size())
{
t += phi[2]*this->psi_.boundaryField()[patchi]
[this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
[this->mesh_.boundaryMesh()[patchi].whichFace(facei)];
}
else
{
@ -336,14 +336,14 @@ Type Foam::interpolationCellPointFace<Type>::interpolate
}
else
{
label patchi = this->pMesh_.boundaryMesh().whichPatch(facei);
label patchi = this->mesh_.boundaryMesh().whichPatch(facei);
// If the boundary patch is not empty use the face value
// else use the cell value
if (this->psi_.boundaryField()[patchi].size())
{
t = this->psi_.boundaryField()[patchi]
[this->pMesh_.boundaryMesh()[patchi].whichFace(facei)];
[this->mesh_.boundaryMesh()[patchi].whichFace(facei)];
}
else
{

5
src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.H
View File

@ -46,7 +46,7 @@ namespace Foam
template<class Type>
class interpolationCellPointFace
:
public interpolation<Type>
public fieldInterpolation<Type, interpolationCellPointFace<Type>>
{
// Private Data
@ -95,6 +95,9 @@ public:
// Member Functions
//- Inherit interpolate from interpolation
using interpolation<Type>::interpolate;
//- Interpolate field to the given point in the given cell
Type interpolate
(

4
src/finiteVolume/interpolation/interpolation/interpolationPointMVC/interpolationPointMVC.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,7 +34,7 @@ Foam::interpolationPointMVC<Type>::interpolationPointMVC
const GeometricField<Type, fvPatchField, volMesh>& psi
)
:
interpolation<Type>(psi),
fieldInterpolation<Type, interpolationPointMVC<Type>>(psi),
psip_
(
volPointInterpolation::New(psi.mesh()).interpolate

2
src/finiteVolume/interpolation/interpolation/interpolationPointMVC/interpolationPointMVC.H
View File

@ -48,7 +48,7 @@ namespace Foam
template<class Type>
class interpolationPointMVC
:
public interpolation<Type>
public fieldInterpolation<Type, interpolationPointMVC<Type>>
{
protected:

4
src/finiteVolume/interpolation/interpolation/interpolationPointMVC/interpolationPointMVCI.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -45,7 +45,7 @@ inline Type Foam::interpolationPointMVC<Type>::interpolate
{
return interpolate
(
pointMVCWeight(this->pMesh_, position, celli, facei)
pointMVCWeight(this->mesh_, position, celli, facei)
);
}