diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 3084f3b55..49b0a6b72 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -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 diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSet.C b/src/finiteVolume/cfdTools/general/levelSet/levelSet.C new file mode 100644 index 000000000..83057c56a --- /dev/null +++ b/src/finiteVolume/cfdTools/general/levelSet/levelSet.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#include "levelSet.H" +#include "cut.H" +#include "polyMeshTetDecomposition.H" +#include "tetIndices.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::tmp> +Foam::levelSetFraction +( + const fvMesh& mesh, + const scalarField& levelC, + const scalarField& levelP, + const bool above +) +{ + DimensionedField sum + ( + IOobject + ( + "levelSetIntegral", + mesh.time().timeName(), + mesh + ), + mesh, + dimensionedScalar("0", dimVolume, 0) + ); + + forAll(sum, cI) + { + const List 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 + tet = + { + mesh.cellCentres()[cI], + mesh.points()[pI0], + mesh.points()[pIA], + mesh.points()[pIB] + }; + const FixedList + 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::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 + tri = + { + patch.patch().faceCentres()[fI], + patch.patch().localPoints()[e[0]], + patch.patch().localPoints()[e[1]] + }; + const FixedList + 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()); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSet.H b/src/finiteVolume/cfdTools/general/levelSet/levelSet.H new file mode 100644 index 000000000..f7b7004a4 --- /dev/null +++ b/src/finiteVolume/cfdTools/general/levelSet/levelSet.H @@ -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 . + +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 +tmp> levelSetAverage +( + const fvMesh& mesh, + const scalarField& levelC, + const scalarField& levelP, + const DimensionedField& positiveC, + const DimensionedField& positiveP, + const DimensionedField& negativeC, + const DimensionedField& negativeP +); + +//- As the above overload, but on the faces of a patch +template +tmp> levelSetAverage +( + const fvPatch& patch, + const scalarField& levelF, + const scalarField& levelP, + const Field& positiveF, + const Field& positiveP, + const Field& negativeF, + const Field& 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> 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 levelSetFraction +( + const fvPatch& patch, + const scalarField& levelF, + const scalarField& levelP, + const bool above +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "levelSetTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C b/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C new file mode 100644 index 000000000..9ac490efb --- /dev/null +++ b/src/finiteVolume/cfdTools/general/levelSet/levelSetTemplates.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#include "levelSet.H" +#include "cut.H" +#include "polyMeshTetDecomposition.H" +#include "tetIndices.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +Foam::tmp> Foam::levelSetAverage +( + const fvMesh& mesh, + const scalarField& levelC, + const scalarField& levelP, + const DimensionedField& positiveC, + const DimensionedField& positiveP, + const DimensionedField& negativeC, + const DimensionedField& negativeP +) +{ + DimensionedField sum + ( + IOobject + ( + "levelSetIntegral", + mesh.time().timeName(), + mesh + ), + mesh, + dimensioned + ( + "0", + positiveC.dimensions()*dimVolume, + pTraits::zero + ) + ); + + forAll(sum, cI) + { + const List 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 + tet = + { + mesh.cellCentres()[cI], + mesh.points()[pI0], + mesh.points()[pIA], + mesh.points()[pIB] + }; + const FixedList + level = + { + levelC[cI], + levelP[pI0], + levelP[pIA], + levelP[pIA] + }; + const cut::volumeIntegrateOp + positive = + { + positiveC[cI], + positiveP[pI0], + positiveP[pIA], + positiveP[pIB] + }; + const cut::volumeIntegrateOp + negative = + { + negativeC[cI], + negativeP[pI0], + negativeP[pIA], + negativeP[pIB] + }; + + sum[cI] += tetCut(tet, level, positive, negative); + } + } + + return sum/mesh.V(); +} + + +template +Foam::tmp> Foam::levelSetAverage +( + const fvPatch& patch, + const scalarField& levelF, + const scalarField& levelP, + const Field& positiveF, + const Field& positiveP, + const Field& negativeF, + const Field& negativeP +) +{ + typedef typename outerProduct::type sumType; + + Field sum(patch.size(), pTraits::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 + tri = + { + patch.patch().faceCentres()[fI], + patch.patch().localPoints()[e[0]], + patch.patch().localPoints()[e[1]] + }; + const FixedList + level = + { + levelF[fI], + levelP[e[0]], + levelP[e[1]] + }; + const cut::areaIntegrateOp + positive = + { + positiveF[fI], + positiveP[e[0]], + positiveP[e[1]] + }; + const cut::areaIntegrateOp + negative = + { + negativeF[fI], + negativeP[e[0]], + negativeP[e[1]] + }; + + sum[fI] += triCut(tri, level, positive, negative); + } + } + + return sum & patch.Sf()/sqr(patch.magSf()); +} + + +// ************************************************************************* //