Skip to content
Permalink
Browse files

interpolationCellPointWallModified: Restored interpolation method

This interpolation method was previously removed by commit fbf0020.

The intention of this method is to provide a slip-like wall boundary
condition for the velocity when interpolated to the location of a
Lagrangian element. This is difficult because any velocity which points
through the wall can cause a drag model and a rebound wall interaction
to "fight"; i.e., the drag pushes the particle to the wall, the wall
bounces it back. This can result in the program hanging.

This method extrapolates a vector field to the wall points and then
modifies the result so that it does not point through the wall. It does
this by rotating the vectors towards the (reversed) point normal. The
result is also scaled so that is reduced to zero if the necessary
rotation exceeds 90 degrees.

This provides an alternate resolution to bug report
https://bugs.openfoam.org/view.php?id=2826
  • Loading branch information...
Will Bainbridge
Will Bainbridge committed Oct 12, 2018
1 parent 7ead70d commit 63b641a068c7fcbb08fb0ccffa3e3cf2d9a9413f
@@ -253,6 +253,7 @@ $(interpolation)/interpolationCellPatchConstrained/makeInterpolationCellPatchCon
$(interpolation)/interpolationCellPoint/cellPointWeight/cellPointWeight.C
$(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C
$(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C
$(interpolation)/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
$(interpolation)/interpolationPointMVC/pointMVCWeight.C
$(interpolation)/interpolationPointMVC/makeInterpolationPointMVC.C

@@ -50,4 +50,19 @@ Foam::interpolationCellPoint<Type>::interpolationCellPoint
}


template<class Type>
Foam::interpolationCellPoint<Type>::interpolationCellPoint
(
const GeometricField<Type, fvPatchField, volMesh>& psi,
tmp<GeometricField<Type, pointPatchField, pointMesh>> psip
)
:
interpolation<Type>(psi),
psip_(psip)
{
// Uses cellPointWeight to do interpolation which needs tet decomposition
(void)psi.mesh().tetBasePtIs();
}


// ************************************************************************* //
@@ -72,6 +72,13 @@ public:
const GeometricField<Type, fvPatchField, volMesh>& psi
);

//- Construct from components
interpolationCellPoint
(
const GeometricField<Type, fvPatchField, volMesh>& psi,
tmp<GeometricField<Type, pointPatchField, pointMesh>> psip
);


// Member Functions

@@ -0,0 +1,339 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "interpolationCellPointWallModified.H"
#include "syncTools.H"
#include "volPointInterpolation.H"
#include "wallPolyPatch.H"
#include "zeroGradientFvPatchFields.H"

// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //

template<class Type>
template<class TYPE>
Foam::tmp<Foam::GeometricField<TYPE, Foam::pointPatchField, Foam::pointMesh>>
Foam::interpolationCellPointWallModified<Type>::calcPointField
(
const GeometricField<TYPE, fvPatchField, volMesh>& psi
) const
{
FatalErrorInFunction
<< typeName << " interpolation is only defined for vector fields"
<< exit(FatalError);

return tmp<GeometricField<TYPE, pointPatchField, pointMesh>>(nullptr);
}


template<class Type>
Foam::tmp<Foam::pointVectorField>
Foam::interpolationCellPointWallModified<Type>::calcPointField
(
const volVectorField& psi
) const
{
using namespace constant::mathematical;

const fvMesh& mesh = psi.mesh();


// Create the point field
tmp<pointVectorField> tPsip
(
new pointVectorField
(
IOobject
(
"volPointInterpolateWallModified(" + psi.name() + ')',
mesh.time().timeName(),
mesh
),
pointMesh::New(mesh),
dimensioned<Type>("zero", psi.dimensions(), Zero)
)
);
pointVectorField& psip = tPsip.ref();


// Interpolate to the points with wall patches extrapolated
{
wordList patchTypes(psi.boundaryField().size());
forAll(patchTypes, patchi)
{
if (isA<wallPolyPatch>(mesh.boundaryMesh()[patchi]))
{
patchTypes[patchi] = zeroGradientFvPatchVectorField::typeName;
}
else
{
patchTypes[patchi] = calculatedFvPatchVectorField::typeName;
}
}
volVectorField psiExtrapolated
(
IOobject
(
psi.name() + "Extrapolated",
mesh.time().timeName(),
mesh
),
psi,
patchTypes
);
psiExtrapolated.correctBoundaryConditions();
volPointInterpolation::New(mesh).interpolate(psiExtrapolated, psip);
}


// Generate point normals across the entire boundary
pointField pointNormals(mesh.nPoints(), vector::zero);
{
scalarField pointCount(mesh.nPoints(), 0);

forAll(mesh.boundaryMesh(), patchi)
{
const polyPatch& patch = mesh.boundaryMesh()[patchi];

forAll(patch, patchFacei)
{
const face& f = patch[patchFacei];
const vector& n = patch.faceNormals()[patchFacei];

forAll(f, i)
{
pointNormals[f[i]] += n;
pointCount[f[i]] += 1;
}
}
}

syncTools::syncPointList
(
mesh,
pointNormals,
plusEqOp<vector>(),
vector::zero
);

syncTools::syncPointList
(
mesh,
pointCount,
plusEqOp<scalar>(),
scalar(0)
);

pointNormals /= max(pointCount, small);
}


// Calculate the rotation necessary from the vector to the (negative) point
// normal needed to make the interpolated field point into the mesh
scalarField theta0(mesh.nPoints(), -vGreat), theta1(mesh.nPoints(), vGreat);
scalar maxVHatDotN = - vGreat;
forAll(mesh.boundaryMesh(), patchi)
{
const polyPatch& patch = mesh.boundaryMesh()[patchi];
if (isA<wallPolyPatch>(patch))
{
forAll(patch, patchFacei)
{
const label facei = patch.start() + patchFacei;
const face& f = patch[patchFacei];

for (label i = 1; i < f.size() - 1; ++ i)
{
const triFace tri
(
tetIndices(mesh.faceOwner()[facei], facei, i)
.faceTriIs(mesh)
);

const vector n = tri.normal(mesh.points());

forAll(tri, triPointI)
{
const label pointi = tri[triPointI];

const vector& v = psip[pointi];

const scalar vHatDotN = normalised(v) & n;
maxVHatDotN = max(maxVHatDotN, vHatDotN);

const vector a =
normalised(v ^ pointNormals[pointi]);

const scalar C = v & n, S = (v ^ a) & n;

const scalar theta = atan2(C, - S);

theta0[pointi] = max(theta0[pointi], theta);
theta1[pointi] = min(theta1[pointi], theta + pi);
}
}
}
}
}

if (debug)
{
reduce(maxVHatDotN, maxOp<scalar>());
Info<< typeName << ": Maximum in-to-wall dot product before = "
<< maxVHatDotN << endl;
}

syncTools::syncPointList
(
mesh,
theta0,
maxEqOp<scalar>(),
scalar(0)
);

syncTools::syncPointList
(
mesh,
theta1,
minEqOp<scalar>(),
scalar(0)
);


if (debug > 1)
{
pointVectorField(psip.name() + "Before", psip).write();
}


// Apply the rotations so that the interpolated field points into the mesh
forAll(mesh.boundaryMesh(), patchi)
{
const polyPatch& patch = mesh.boundaryMesh()[patchi];
if (isA<wallPolyPatch>(patch))
{
forAll(patch.meshPoints(), patchPointi)
{
const label pointi = patch.meshPoints()[patchPointi];

vector& v = psip[pointi];

if (theta0[pointi] <= 0 && theta1[pointi] >= 0)
{
continue;
}

if (theta0[pointi] >= theta1[pointi])
{
v = Zero;
continue;
}

const scalar theta =
theta0[pointi] > 0 ? theta0[pointi] : theta1[pointi];

const scalar c = cos(theta), s = sin(theta);

const scalar scale = max(c, 0); // or mag(theta)/(pi/2) ...

const vector a = normalised(v ^ pointNormals[pointi]);

v = scale*(tensor::I*c - (*a)*s + sqr(a)*(1 - c)) & v;

theta0[pointi] = theta1[pointi] = 0;
}
}
}


// Report the field-normal dot products
if (debug)
{
maxVHatDotN = - vGreat;

forAll(mesh.boundaryMesh(), patchi)
{
const polyPatch& patch = mesh.boundaryMesh()[patchi];
if (isA<wallPolyPatch>(patch))
{
forAll(patch, patchFacei)
{
const label facei = patch.start() + patchFacei;
const face& f = patch[patchFacei];

for (label i = 1; i < f.size() - 1; ++ i)
{
const triFace tri
(
tetIndices(mesh.faceOwner()[facei], facei, i)
.faceTriIs(mesh)
);

const vector n = tri.normal(mesh.points());

forAll(tri, triPointI)
{
const label pointi = tri[triPointI];

const vector& v = psip[pointi];

const scalar vHatDotN = normalised(v) & n;
maxVHatDotN = max(maxVHatDotN, vHatDotN);
}
}
}
}
}

reduce(maxVHatDotN, maxOp<scalar>());
Info<< typeName << ": Maximum in-to-wall dot product after = "
<< maxVHatDotN << endl;
}


// Write the interpolated field
if (debug > 1)
{
psip.write();
}


return tPsip;
}


// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //

template<class Type>
Foam::interpolationCellPointWallModified<Type>::
interpolationCellPointWallModified
(
const GeometricField<Type, fvPatchField, volMesh>& psi
)
:
interpolationCellPoint<Type>(psi, calcPointField(psi))
{}


// ************************************************************************* //

0 comments on commit 63b641a

Please sign in to comment.
You can’t perform that action at this time.