ec6c0d7e0f
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.
160 lines
4.4 KiB
C
160 lines
4.4 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / 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
|
|
)
|
|
{
|
|
tmp<DimensionedField<scalar, volMesh>> tResult
|
|
(
|
|
new DimensionedField<scalar, volMesh>
|
|
(
|
|
IOobject
|
|
(
|
|
"levelSetFraction",
|
|
mesh.time().timeName(),
|
|
mesh
|
|
),
|
|
mesh,
|
|
dimensionedScalar("0", dimless, 0)
|
|
)
|
|
);
|
|
DimensionedField<scalar, volMesh>& result = tResult.ref();
|
|
|
|
forAll(result, cI)
|
|
{
|
|
const List<tetIndices> cellTetIs =
|
|
polyMeshTetDecomposition::cellTetIndices(mesh, cI);
|
|
|
|
scalar v = 0, r = 0;
|
|
|
|
forAll(cellTetIs, cellTetI)
|
|
{
|
|
const triFace triIs = cellTetIs[cellTetI].faceTriIs(mesh);
|
|
|
|
const FixedList<point, 4>
|
|
tet =
|
|
{
|
|
mesh.cellCentres()[cI],
|
|
mesh.points()[triIs[0]],
|
|
mesh.points()[triIs[1]],
|
|
mesh.points()[triIs[2]]
|
|
};
|
|
const FixedList<scalar, 4>
|
|
level =
|
|
{
|
|
levelC[cI],
|
|
levelP[triIs[0]],
|
|
levelP[triIs[1]],
|
|
levelP[triIs[2]]
|
|
};
|
|
|
|
v += cut::volumeOp()(tet);
|
|
|
|
if (above)
|
|
{
|
|
r += tetCut(tet, level, cut::volumeOp(), cut::noOp());
|
|
}
|
|
else
|
|
{
|
|
r += tetCut(tet, level, cut::noOp(), cut::volumeOp());
|
|
}
|
|
}
|
|
|
|
result[cI] = r/v;
|
|
}
|
|
|
|
return tResult;
|
|
}
|
|
|
|
|
|
Foam::tmp<Foam::scalarField> Foam::levelSetFraction
|
|
(
|
|
const fvPatch& patch,
|
|
const scalarField& levelF,
|
|
const scalarField& levelP,
|
|
const bool above
|
|
)
|
|
{
|
|
tmp<scalarField> tResult(new scalarField(patch.size(), 0));
|
|
scalarField& result = tResult.ref();
|
|
|
|
forAll(result, fI)
|
|
{
|
|
const face& f = patch.patch().localFaces()[fI];
|
|
|
|
vector a = vector::zero, r = vector::zero;
|
|
|
|
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]]
|
|
};
|
|
|
|
a += cut::areaOp()(tri);
|
|
|
|
if (above)
|
|
{
|
|
r += triCut(tri, level, cut::areaOp(), cut::noOp());
|
|
}
|
|
else
|
|
{
|
|
r += triCut(tri, level, cut::noOp(), cut::areaOp());
|
|
}
|
|
}
|
|
|
|
result[fI] = a/magSqr(a) & r;
|
|
}
|
|
|
|
return tResult;
|
|
}
|
|
|
|
// ************************************************************************* //
|