cfdTools: Added a number of functions for performing volume averages of
discontinuous fields, with the discontinuity defined by a level set. The functions do a proper integration of the discontinuous fields by tet- and tri-cutting along the plane of the level set.
This commit is contained in:
Will Bainbridge
@ -404,6 +404,7 @@ $(general)/adjustPhi/adjustPhi.C
|
||||
$(general)/bound/bound.C
|
||||
$(general)/CorrectPhi/correctUphiBCs.C
|
||||
$(general)/pressureControl/pressureControl.C
|
||||
$(general)/levelSet/levelSet.C
|
||||
|
||||
solutionControl = $(general)/solutionControl
|
||||
$(solutionControl)/solutionControl/solutionControl.C
|
||||
|
||||
147
src/finiteVolume/cfdTools/general/levelSet/levelSet.C
Normal file
147
src/finiteVolume/cfdTools/general/levelSet/levelSet.C
Normal file
@ -0,0 +1,147 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2017 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 "levelSet.H"
|
||||
#include "cut.H"
|
||||
#include "polyMeshTetDecomposition.H"
|
||||
#include "tetIndices.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh>>
|
||||
Foam::levelSetFraction
|
||||
(
|
||||
const fvMesh& mesh,
|
||||
const scalarField& levelC,
|
||||
const scalarField& levelP,
|
||||
const bool above
|
||||
)
|
||||
{
|
||||
DimensionedField<scalar, volMesh> sum
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"levelSetIntegral",
|
||||
mesh.time().timeName(),
|
||||
mesh
|
||||
),
|
||||
mesh,
|
||||
dimensionedScalar("0", dimVolume, 0)
|
||||
);
|
||||
|
||||
forAll(sum, cI)
|
||||
{
|
||||
const List<tetIndices> cellTetIs =
|
||||
polyMeshTetDecomposition::cellTetIndices(mesh, cI);
|
||||
|
||||
forAll(cellTetIs, cellTetI)
|
||||
{
|
||||
const tetIndices& tetIs = cellTetIs[cellTetI];
|
||||
const face& f = mesh.faces()[tetIs.face()];
|
||||
|
||||
const label pI0 = f[tetIs.faceBasePt()];
|
||||
const label pIA = f[tetIs.facePtA()];
|
||||
const label pIB = f[tetIs.facePtB()];
|
||||
|
||||
const FixedList<point, 4>
|
||||
tet =
|
||||
{
|
||||
mesh.cellCentres()[cI],
|
||||
mesh.points()[pI0],
|
||||
mesh.points()[pIA],
|
||||
mesh.points()[pIB]
|
||||
};
|
||||
const FixedList<scalar, 4>
|
||||
level =
|
||||
{
|
||||
levelC[cI],
|
||||
levelP[pI0],
|
||||
levelP[pIA],
|
||||
levelP[pIA]
|
||||
};
|
||||
|
||||
if (above)
|
||||
{
|
||||
sum[cI] += tetCut(tet, level, cut::volumeOp(), cut::noOp());
|
||||
}
|
||||
else
|
||||
{
|
||||
sum[cI] += tetCut(tet, level, cut::noOp(), cut::volumeOp());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return sum/mesh.V();
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::levelSetFraction
|
||||
(
|
||||
const fvPatch& patch,
|
||||
const scalarField& levelF,
|
||||
const scalarField& levelP,
|
||||
const bool above
|
||||
)
|
||||
{
|
||||
vectorField sum(patch.size(), vector::zero);
|
||||
|
||||
forAll(sum, fI)
|
||||
{
|
||||
const face& f = patch.patch().localFaces()[fI];
|
||||
|
||||
for(label eI = 0; eI < f.size(); ++ eI)
|
||||
{
|
||||
const edge e = f.faceEdge(eI);
|
||||
|
||||
const FixedList<point, 3>
|
||||
tri =
|
||||
{
|
||||
patch.patch().faceCentres()[fI],
|
||||
patch.patch().localPoints()[e[0]],
|
||||
patch.patch().localPoints()[e[1]]
|
||||
};
|
||||
const FixedList<scalar, 3>
|
||||
level =
|
||||
{
|
||||
levelF[fI],
|
||||
levelP[e[0]],
|
||||
levelP[e[1]]
|
||||
};
|
||||
|
||||
if (above)
|
||||
{
|
||||
sum[fI] += triCut(tri, level, cut::areaOp(), cut::noOp());
|
||||
}
|
||||
else
|
||||
{
|
||||
sum[fI] += triCut(tri, level, cut::noOp(), cut::areaOp());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return sum & patch.Sf()/sqr(patch.magSf());
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
110
src/finiteVolume/cfdTools/general/levelSet/levelSet.H
Normal file
110
src/finiteVolume/cfdTools/general/levelSet/levelSet.H
Normal file
@ -0,0 +1,110 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2017 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/>.
|
||||
|
||||
Description:
|
||||
Functions to create an averaged field from a discontinuous field defined by
|
||||
a level-set.
|
||||
|
||||
SourceFiles:
|
||||
levelSet.C
|
||||
levelSetTemplates.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef levelSet_H
|
||||
#define levelSet_H
|
||||
|
||||
#include "DimensionedField.H"
|
||||
#include "fvMesh.H"
|
||||
#include "pointMesh.H"
|
||||
#include "volMesh.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
//- Calculate the average value of two fields, one on each side of a level set
|
||||
// The level set and the fields are given both on the points and cell-centres.
|
||||
template<class Type>
|
||||
tmp<DimensionedField<Type, volMesh>> levelSetAverage
|
||||
(
|
||||
const fvMesh& mesh,
|
||||
const scalarField& levelC,
|
||||
const scalarField& levelP,
|
||||
const DimensionedField<Type, volMesh>& positiveC,
|
||||
const DimensionedField<Type, pointMesh>& positiveP,
|
||||
const DimensionedField<Type, volMesh>& negativeC,
|
||||
const DimensionedField<Type, pointMesh>& negativeP
|
||||
);
|
||||
|
||||
//- As the above overload, but on the faces of a patch
|
||||
template<class Type>
|
||||
tmp<Field<Type>> levelSetAverage
|
||||
(
|
||||
const fvPatch& patch,
|
||||
const scalarField& levelF,
|
||||
const scalarField& levelP,
|
||||
const Field<Type>& positiveF,
|
||||
const Field<Type>& positiveP,
|
||||
const Field<Type>& negativeF,
|
||||
const Field<Type>& negativeP
|
||||
);
|
||||
|
||||
//- Calculate the volume-fraction that a level set occupies. This gives the the
|
||||
// same result as levelSetAverage if the fields passed to the latter are
|
||||
// uniformly 0 and 1. The above flag flips the direction.
|
||||
tmp<DimensionedField<scalar, volMesh>> levelSetFraction
|
||||
(
|
||||
const fvMesh& mesh,
|
||||
const scalarField& levelC,
|
||||
const scalarField& levelP,
|
||||
const bool above
|
||||
);
|
||||
|
||||
//- As the above oveload, but on the faces of a patch
|
||||
tmp<scalarField> levelSetFraction
|
||||
(
|
||||
const fvPatch& patch,
|
||||
const scalarField& levelF,
|
||||
const scalarField& levelP,
|
||||
const bool above
|
||||
);
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#ifdef NoRepository
|
||||
#include "levelSetTemplates.C"
|
||||
#endif
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
178
src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
Normal file
178
src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C
Normal file
@ -0,0 +1,178 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2017 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 "levelSet.H"
|
||||
#include "cut.H"
|
||||
#include "polyMeshTetDecomposition.H"
|
||||
#include "tetIndices.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
template<class Type>
|
||||
Foam::tmp<Foam::DimensionedField<Type, Foam::volMesh>> Foam::levelSetAverage
|
||||
(
|
||||
const fvMesh& mesh,
|
||||
const scalarField& levelC,
|
||||
const scalarField& levelP,
|
||||
const DimensionedField<Type, volMesh>& positiveC,
|
||||
const DimensionedField<Type, pointMesh>& positiveP,
|
||||
const DimensionedField<Type, volMesh>& negativeC,
|
||||
const DimensionedField<Type, pointMesh>& negativeP
|
||||
)
|
||||
{
|
||||
DimensionedField<Type, volMesh> sum
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"levelSetIntegral",
|
||||
mesh.time().timeName(),
|
||||
mesh
|
||||
),
|
||||
mesh,
|
||||
dimensioned<Type>
|
||||
(
|
||||
"0",
|
||||
positiveC.dimensions()*dimVolume,
|
||||
pTraits<Type>::zero
|
||||
)
|
||||
);
|
||||
|
||||
forAll(sum, cI)
|
||||
{
|
||||
const List<tetIndices> cellTetIs =
|
||||
polyMeshTetDecomposition::cellTetIndices(mesh, cI);
|
||||
|
||||
forAll(cellTetIs, cellTetI)
|
||||
{
|
||||
const tetIndices& tetIs = cellTetIs[cellTetI];
|
||||
const face& f = mesh.faces()[tetIs.face()];
|
||||
|
||||
const label pI0 = f[tetIs.faceBasePt()];
|
||||
const label pIA = f[tetIs.facePtA()];
|
||||
const label pIB = f[tetIs.facePtB()];
|
||||
|
||||
const FixedList<point, 4>
|
||||
tet =
|
||||
{
|
||||
mesh.cellCentres()[cI],
|
||||
mesh.points()[pI0],
|
||||
mesh.points()[pIA],
|
||||
mesh.points()[pIB]
|
||||
};
|
||||
const FixedList<scalar, 4>
|
||||
level =
|
||||
{
|
||||
levelC[cI],
|
||||
levelP[pI0],
|
||||
levelP[pIA],
|
||||
levelP[pIA]
|
||||
};
|
||||
const cut::volumeIntegrateOp<Type>
|
||||
positive =
|
||||
{
|
||||
positiveC[cI],
|
||||
positiveP[pI0],
|
||||
positiveP[pIA],
|
||||
positiveP[pIB]
|
||||
};
|
||||
const cut::volumeIntegrateOp<Type>
|
||||
negative =
|
||||
{
|
||||
negativeC[cI],
|
||||
negativeP[pI0],
|
||||
negativeP[pIA],
|
||||
negativeP[pIB]
|
||||
};
|
||||
|
||||
sum[cI] += tetCut(tet, level, positive, negative);
|
||||
}
|
||||
}
|
||||
|
||||
return sum/mesh.V();
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
Foam::tmp<Foam::Field<Type>> Foam::levelSetAverage
|
||||
(
|
||||
const fvPatch& patch,
|
||||
const scalarField& levelF,
|
||||
const scalarField& levelP,
|
||||
const Field<Type>& positiveF,
|
||||
const Field<Type>& positiveP,
|
||||
const Field<Type>& negativeF,
|
||||
const Field<Type>& negativeP
|
||||
)
|
||||
{
|
||||
typedef typename outerProduct<Type, vector>::type sumType;
|
||||
|
||||
Field<sumType> sum(patch.size(), pTraits<sumType>::zero);
|
||||
|
||||
forAll(sum, fI)
|
||||
{
|
||||
const face& f = patch.patch().localFaces()[fI];
|
||||
|
||||
for(label eI = 0; eI < f.size(); ++ eI)
|
||||
{
|
||||
const edge e = f.faceEdge(eI);
|
||||
|
||||
const FixedList<point, 3>
|
||||
tri =
|
||||
{
|
||||
patch.patch().faceCentres()[fI],
|
||||
patch.patch().localPoints()[e[0]],
|
||||
patch.patch().localPoints()[e[1]]
|
||||
};
|
||||
const FixedList<scalar, 3>
|
||||
level =
|
||||
{
|
||||
levelF[fI],
|
||||
levelP[e[0]],
|
||||
levelP[e[1]]
|
||||
};
|
||||
const cut::areaIntegrateOp<Type>
|
||||
positive =
|
||||
{
|
||||
positiveF[fI],
|
||||
positiveP[e[0]],
|
||||
positiveP[e[1]]
|
||||
};
|
||||
const cut::areaIntegrateOp<Type>
|
||||
negative =
|
||||
{
|
||||
negativeF[fI],
|
||||
negativeP[e[0]],
|
||||
negativeP[e[1]]
|
||||
};
|
||||
|
||||
sum[fI] += triCut(tri, level, positive, negative);
|
||||
}
|
||||
}
|
||||
|
||||
return sum & patch.Sf()/sqr(patch.magSf());
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user