mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
ENH: interpolationPoint : interpolation from point values only (using Mean Value Coordinates)
This commit is contained in:
mattijs
@ -196,6 +196,8 @@ $(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C
|
||||
$(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C
|
||||
$(interpolation)/interpolationCellPointWallModified/cellPointWeightWallModified/cellPointWeightWallModified.C
|
||||
$(interpolation)/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
|
||||
$(interpolation)/interpolationPoint/pointMVCWeight.C
|
||||
$(interpolation)/interpolationPoint/makeInterpolationPoint.C
|
||||
|
||||
volPointInterpolation = interpolation/volPointInterpolation
|
||||
/*
|
||||
|
||||
@ -0,0 +1,255 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
|
||||
\\/ 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 "interpolationPoint.H"
|
||||
#include "volPointInterpolation.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
|
||||
|
||||
template<class Type>
|
||||
Foam::interpolationPoint<Type>::interpolationPoint
|
||||
(
|
||||
const GeometricField<Type, fvPatchField, volMesh>& psi
|
||||
)
|
||||
:
|
||||
interpolation<Type>(psi),
|
||||
psip_(volPointInterpolation::New(psi.mesh()).interpolate(psi))
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
//template<class Type>
|
||||
//void Foam::interpolationPoint<Type>::calcWeights
|
||||
//(
|
||||
// const vector& position,
|
||||
// const label cellI,
|
||||
// const label faceI,
|
||||
// scalarField& weights
|
||||
//) const
|
||||
//{
|
||||
// const polyMesh& mesh = this->pMesh_;
|
||||
// const pointField& points = mesh.points();
|
||||
//
|
||||
//
|
||||
// const scalar eps = 0.00000001;
|
||||
//
|
||||
//
|
||||
// // Addressing - face vertices to local points
|
||||
// const labelList& toGlobal = mesh.cellPoints()[cellI];
|
||||
// Map<label> toLocal(2*toGlobal.size());
|
||||
// forAll(toGlobal, i)
|
||||
// {
|
||||
// toLocal.insert(toGlobal[i], i);
|
||||
// }
|
||||
//
|
||||
// // Initialise weights
|
||||
// weights.setSize(toGlobal.size());
|
||||
// weights = 0.0;
|
||||
//
|
||||
// // Point-to-vertex vectors and distances
|
||||
// scalarField dist(toGlobal.size());
|
||||
// vectorField uVec(toGlobal.size());
|
||||
// forAll(toGlobal, pid)
|
||||
// {
|
||||
// const point& pt = points[toGlobal[pid]];//-cc;
|
||||
// uVec[pid] = pt-position;
|
||||
// dist[pid] = mag(uVec[pid]);
|
||||
//
|
||||
// // Special case: point is close to vertex
|
||||
// if (dist[pid] < eps)
|
||||
// {
|
||||
// weights[pid] = 1.0;
|
||||
// return;
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// // Project onto unit sphere
|
||||
// uVec /= dist;
|
||||
//
|
||||
//
|
||||
// // Loop over all triangles of all polygons of cell to compute weights
|
||||
// DynamicList<scalar> alpha(100);
|
||||
// DynamicList<scalar> theta(100);
|
||||
//
|
||||
// const cell& cFaces = mesh.cells()[cellI];
|
||||
//
|
||||
// forAll(cFaces, iter)
|
||||
// {
|
||||
// label faceI = cFaces[iter];
|
||||
// const face& f = mesh.faces()[faceI];
|
||||
//
|
||||
// Pout<< "face:" << faceI << " at:"
|
||||
// << pointField(mesh.points(), f)
|
||||
// << endl;
|
||||
//
|
||||
// vector v(point::zero);
|
||||
// forAll(f, j)
|
||||
// {
|
||||
// label jPlus1 = f.fcIndex(j);
|
||||
// const point& uj = points[f[j]];//-cc;
|
||||
// const point& ujPlus1 = points[f[jPlus1]];//-cc;
|
||||
// Pout<< " uj:" << uj << " ujPlus1:" << ujPlus1 << endl;
|
||||
//
|
||||
// vector temp = uj ^ ujPlus1;
|
||||
// temp /= mag(temp);
|
||||
//
|
||||
// scalar l = mag(uj-ujPlus1);
|
||||
// scalar angle = 2.0*Foam::asin(l/2.0);
|
||||
//
|
||||
// v += 0.5*angle*temp;
|
||||
// }
|
||||
//
|
||||
// scalar vNorm = mag(v);
|
||||
// v /= vNorm;
|
||||
//
|
||||
// // Make sure v points towards the polygon
|
||||
// if ((v&points[f[0]]) < 0)
|
||||
// {
|
||||
// v = -v;
|
||||
// }
|
||||
//
|
||||
// Pout<< " v:" << v << endl;
|
||||
//
|
||||
// // angles between edges
|
||||
// forAll(f, j)
|
||||
// {
|
||||
// label jPlus1 = f.fcIndex(j);
|
||||
// const point& uj = points[f[j]];//-cc;
|
||||
// const point& ujPlus1 = points[f[jPlus1]];//-cc;
|
||||
// Pout<< " uj:" << uj << " ujPlus1:" << ujPlus1 << endl;
|
||||
//
|
||||
// vector n0 = uj ^ v;
|
||||
// n0 /= mag(n0);
|
||||
// vector n1 = ujPlus1 ^ v;
|
||||
// n1 /= mag(n1);
|
||||
//
|
||||
// scalar l = mag(n0-n1);
|
||||
// Pout<< " l:" << l << endl;
|
||||
// alpha(j) = 2.0*Foam::asin(l/2.0);
|
||||
//
|
||||
// vector temp = n0 ^ n1;
|
||||
// if ((temp&v) < 0.0)
|
||||
// {
|
||||
// alpha(j) = -alpha(j);
|
||||
// }
|
||||
//
|
||||
// l = mag(uj-v);
|
||||
// Pout<< " l:" << l << endl;
|
||||
// theta(j) = 2.0*Foam::asin(l/2.0);
|
||||
// }
|
||||
//
|
||||
//
|
||||
// bool outlierFlag = false;
|
||||
// forAll(f, j)
|
||||
// {
|
||||
// if (mag(theta(j)) < eps)
|
||||
// {
|
||||
// outlierFlag = true;
|
||||
//
|
||||
// label pid = toLocal[f[j]];
|
||||
// weights[pid] += vNorm / dist[pid];
|
||||
// break;
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// if (outlierFlag)
|
||||
// {
|
||||
// continue;
|
||||
// }
|
||||
//
|
||||
// scalar sum = 0.0;
|
||||
// forAll(f, j)
|
||||
// {
|
||||
// label jMin1 = f.rcIndex(j);
|
||||
// sum +=
|
||||
// 1.0
|
||||
// / Foam::tan(theta(j))
|
||||
// * (Foam::tan(alpha(j)/2.0)+Foam::tan(alpha(jMin1)/2.0));
|
||||
// }
|
||||
//
|
||||
// // The special case when x lies on the polygon, handle it using 2D mvc.
|
||||
// // In the 2D case, alpha = theta
|
||||
// if (mag(sum) < eps)
|
||||
// {
|
||||
// weights = 0.0;
|
||||
//
|
||||
// // recompute theta, the theta computed previously are not robust
|
||||
// forAll(f, j)
|
||||
// {
|
||||
// label jPlus1 = f.fcIndex(j);
|
||||
// const point& uj = points[f[j]];//-cc;
|
||||
// const point& ujPlus1 = points[f[jPlus1]];//-cc;
|
||||
// scalar l = mag(uj-ujPlus1);
|
||||
// theta(j) = 2.0*Foam::asin(l/2.0);
|
||||
// }
|
||||
//
|
||||
// scalar sumWeight = 0;
|
||||
// forAll(f, j)
|
||||
// {
|
||||
// label pid = toLocal[f[j]];
|
||||
// label jMin1 = f.rcIndex(j);
|
||||
// weights[pid] =
|
||||
// 1.0
|
||||
// / dist[pid]
|
||||
// * (Foam::tan(theta(jMin1)/2.0)+Foam::tan(theta(j)/2.0));
|
||||
// sumWeight += weights[pid];
|
||||
// }
|
||||
//
|
||||
// if (sumWeight < eps)
|
||||
// {
|
||||
// return;
|
||||
// }
|
||||
// weights /= sumWeight;
|
||||
// return;
|
||||
// }
|
||||
//
|
||||
//
|
||||
// // Normal 3D case
|
||||
// forAll(f, j)
|
||||
// {
|
||||
// label pid = toLocal[f[j]];
|
||||
// label jMin1 = f.rcIndex(j);
|
||||
// weights[pid] +=
|
||||
// vNorm
|
||||
// / sum
|
||||
// / dist[pid]
|
||||
// / Foam::sin(theta(j))
|
||||
// * (Foam::tan(alpha(j)/2.0)+Foam::tan(alpha(jMin1)/2.0));
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// // normalise weights
|
||||
// scalar sumWeight = sum(weights);
|
||||
//
|
||||
// if (mag(sumWeight) < eps)
|
||||
// {
|
||||
// return;
|
||||
// }
|
||||
// weights /= sumWeight;
|
||||
//}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,108 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
|
||||
\\/ 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/>.
|
||||
|
||||
Class
|
||||
Foam::interpolationPoint
|
||||
|
||||
Description
|
||||
Given cell centre values interpolates to vertices and uses these to
|
||||
do a Mean Value Coordinates interpolation.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef interpolationPoint_H
|
||||
#define interpolationPoint_H
|
||||
|
||||
#include "interpolation.H"
|
||||
#include "pointMVCWeight.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class interpolationPoint Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
template<class Type>
|
||||
class interpolationPoint
|
||||
:
|
||||
public interpolation<Type>
|
||||
{
|
||||
protected:
|
||||
|
||||
// Protected data
|
||||
|
||||
//- Interpolated volfield
|
||||
const GeometricField<Type, pointPatchField, pointMesh> psip_;
|
||||
|
||||
public:
|
||||
|
||||
//- Runtime type information
|
||||
TypeName("cellPoint");
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from components
|
||||
interpolationPoint
|
||||
(
|
||||
const GeometricField<Type, fvPatchField, volMesh>& psi
|
||||
);
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Interpolate field for the given cellPointWeight
|
||||
inline Type interpolate(const pointMVCWeight& cpw) const;
|
||||
|
||||
//- Interpolate field to the given point in the given cell
|
||||
inline Type interpolate
|
||||
(
|
||||
const vector& position,
|
||||
const label nCell,
|
||||
const label facei = -1
|
||||
) const;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#include "interpolationPointI.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#ifdef NoRepository
|
||||
# include "interpolationPoint.C"
|
||||
#endif
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,53 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
|
||||
\\/ 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/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
template<class Type>
|
||||
inline Type Foam::interpolationPoint<Type>::interpolate
|
||||
(
|
||||
const pointMVCWeight& cpw
|
||||
) const
|
||||
{
|
||||
return cpw.interpolate(psip_);
|
||||
}
|
||||
|
||||
|
||||
template<class Type>
|
||||
inline Type Foam::interpolationPoint<Type>::interpolate
|
||||
(
|
||||
const vector& position,
|
||||
const label celli,
|
||||
const label facei
|
||||
) const
|
||||
{
|
||||
return interpolate
|
||||
(
|
||||
pointMVCWeight(this->pMesh_, position, celli, facei)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,35 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
|
||||
\\/ 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 "interpolationPoint.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
makeInterpolation(interpolationPoint);
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,324 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
|
||||
\\/ 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 "pointMVCWeight.H"
|
||||
#include "polyMesh.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
int Foam::pointMVCWeight::debug
|
||||
(
|
||||
debug::debugSwitch("pointMVCWeight", 0)
|
||||
);
|
||||
|
||||
Foam::scalar Foam::pointMVCWeight::tol(SMALL);
|
||||
|
||||
|
||||
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
|
||||
|
||||
void Foam::pointMVCWeight::calcWeights
|
||||
(
|
||||
const Map<label>& toLocal,
|
||||
const face& f,
|
||||
const DynamicList<point>& u,
|
||||
const scalarField& dist,
|
||||
scalarField& weights
|
||||
) const
|
||||
{
|
||||
weights.setSize(toLocal.size());
|
||||
weights = 0.0;
|
||||
|
||||
scalarField theta(f.size());
|
||||
|
||||
// recompute theta, the theta computed previously are not robust
|
||||
forAll(f, j)
|
||||
{
|
||||
label jPlus1 = f.fcIndex(j);
|
||||
scalar l = mag(u[j]-u[jPlus1]);
|
||||
theta[j] = 2.0*Foam::asin(l/2.0);
|
||||
}
|
||||
|
||||
scalar sumWeight = 0;
|
||||
forAll(f, j)
|
||||
{
|
||||
label pid = toLocal[f[j]];
|
||||
label jMin1 = f.rcIndex(j);
|
||||
weights[pid] =
|
||||
1.0
|
||||
/ dist[pid]
|
||||
* (Foam::tan(theta[jMin1]/2.0)+Foam::tan(theta[j]/2.0));
|
||||
sumWeight += weights[pid];
|
||||
}
|
||||
|
||||
if (sumWeight >= tol)
|
||||
{
|
||||
weights /= sumWeight;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::pointMVCWeight::calcWeights
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const labelList& toGlobal,
|
||||
const Map<label>& toLocal,
|
||||
const vector& position,
|
||||
const vectorField& uVec,
|
||||
const scalarField& dist,
|
||||
scalarField& weights
|
||||
) const
|
||||
{
|
||||
// Loop over all triangles of all polygons of cell to compute weights
|
||||
DynamicList<scalar> alpha(100);
|
||||
DynamicList<scalar> theta(100);
|
||||
DynamicList<point> u(100);
|
||||
|
||||
const Foam::cell& cFaces = mesh.cells()[cellIndex_];
|
||||
|
||||
forAll(cFaces, iter)
|
||||
{
|
||||
label faceI = cFaces[iter];
|
||||
const face& f = mesh.faces()[faceI];
|
||||
|
||||
//Pout<< "face:" << faceI << " at:"
|
||||
// << pointField(mesh.points(), f)
|
||||
// << endl;
|
||||
|
||||
// Collect the uVec for the face
|
||||
forAll(f, j)
|
||||
{
|
||||
u(j) = uVec[toLocal[f[j]]];
|
||||
}
|
||||
|
||||
vector v(point::zero);
|
||||
forAll(f, j)
|
||||
{
|
||||
label jPlus1 = f.fcIndex(j);
|
||||
//Pout<< " uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl;
|
||||
|
||||
vector temp = u[j] ^ u[jPlus1];
|
||||
temp /= mag(temp);
|
||||
|
||||
//Pout<< " uj:" << u[j] << " ujPlus1:" << u[jPlus1]
|
||||
// << " temp:" << temp << endl;
|
||||
|
||||
scalar l = mag(u[j]-u[jPlus1]);
|
||||
scalar angle = 2.0*Foam::asin(l/2.0);
|
||||
|
||||
//Pout<< " j:" << j << " l:" << l
|
||||
// << " angle:" << angle << endl;
|
||||
|
||||
v += 0.5*angle*temp;
|
||||
}
|
||||
|
||||
scalar vNorm = mag(v);
|
||||
v /= vNorm;
|
||||
|
||||
// Make sure v points towards the polygon
|
||||
//if (((v&u[0]) < 0) != (mesh.faceOwner()[faceI] != cellIndex_))
|
||||
//{
|
||||
// FatalErrorIn("pointMVCWeight::calcWeights(..)")
|
||||
// << "v:" << v << " u[0]:" << u[0]
|
||||
// << exit(FatalError);
|
||||
//}
|
||||
|
||||
if ((v&u[0]) < 0)
|
||||
{
|
||||
v = -v;
|
||||
}
|
||||
|
||||
//Pout<< " v:" << v << endl;
|
||||
|
||||
// angles between edges
|
||||
forAll(f, j)
|
||||
{
|
||||
label jPlus1 = f.fcIndex(j);
|
||||
//Pout<< " uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl;
|
||||
|
||||
vector n0 = u[j] ^ v;
|
||||
n0 /= mag(n0);
|
||||
vector n1 = u[jPlus1] ^ v;
|
||||
n1 /= mag(n1);
|
||||
|
||||
scalar l = mag(n0-n1);
|
||||
//Pout<< " l:" << l << endl;
|
||||
alpha(j) = 2.0*Foam::asin(l/2.0);
|
||||
|
||||
vector temp = n0 ^ n1;
|
||||
if ((temp&v) < 0.0)
|
||||
{
|
||||
alpha[j] = -alpha[j];
|
||||
}
|
||||
|
||||
l = mag(u[j]-v);
|
||||
//Pout<< " l:" << l << endl;
|
||||
theta(j) = 2.0*Foam::asin(l/2.0);
|
||||
}
|
||||
|
||||
|
||||
bool outlierFlag = false;
|
||||
forAll(f, j)
|
||||
{
|
||||
if (mag(theta[j]) < tol)
|
||||
{
|
||||
outlierFlag = true;
|
||||
|
||||
label pid = toLocal[f[j]];
|
||||
weights[pid] += vNorm / dist[pid];
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (outlierFlag)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
scalar sum = 0.0;
|
||||
forAll(f, j)
|
||||
{
|
||||
label jMin1 = f.rcIndex(j);
|
||||
sum +=
|
||||
1.0
|
||||
/ Foam::tan(theta[j])
|
||||
* (Foam::tan(alpha[j]/2.0)+Foam::tan(alpha[jMin1]/2.0));
|
||||
}
|
||||
|
||||
// The special case when x lies on the polygon, handle it using 2D mvc.
|
||||
// In the 2D case, alpha = theta
|
||||
if (mag(sum) < tol)
|
||||
{
|
||||
// Calculate weights using face vertices only
|
||||
calcWeights(toLocal, f, u, dist, weights);
|
||||
return;
|
||||
}
|
||||
|
||||
|
||||
// Normal 3D case
|
||||
forAll(f, j)
|
||||
{
|
||||
label pid = toLocal[f[j]];
|
||||
label jMin1 = f.rcIndex(j);
|
||||
weights[pid] +=
|
||||
vNorm
|
||||
/ sum
|
||||
/ dist[pid]
|
||||
/ Foam::sin(theta[j])
|
||||
* (Foam::tan(alpha[j]/2.0)+Foam::tan(alpha[jMin1]/2.0));
|
||||
}
|
||||
}
|
||||
|
||||
// normalise weights
|
||||
scalar sumWeight = sum(weights);
|
||||
|
||||
if (mag(sumWeight) < tol)
|
||||
{
|
||||
return;
|
||||
}
|
||||
weights /= sumWeight;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::pointMVCWeight::pointMVCWeight
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const vector& position,
|
||||
const label cellIndex,
|
||||
const label faceIndex
|
||||
)
|
||||
:
|
||||
cellIndex_((cellIndex != -1) ? cellIndex : mesh.faceOwner()[faceIndex])
|
||||
{
|
||||
// Addressing - face vertices to local points and vice versa
|
||||
const labelList& toGlobal = mesh.cellPoints()[cellIndex_];
|
||||
Map<label> toLocal(2*toGlobal.size());
|
||||
forAll(toGlobal, i)
|
||||
{
|
||||
toLocal.insert(toGlobal[i], i);
|
||||
}
|
||||
|
||||
|
||||
// Initialise weights
|
||||
weights_.setSize(toGlobal.size());
|
||||
weights_ = 0.0;
|
||||
|
||||
|
||||
// Point-to-vertex vectors and distances
|
||||
vectorField uVec(toGlobal.size());
|
||||
scalarField dist(toGlobal.size());
|
||||
|
||||
forAll(toGlobal, pid)
|
||||
{
|
||||
const point& pt = mesh.points()[toGlobal[pid]];
|
||||
|
||||
uVec[pid] = pt-position;
|
||||
dist[pid] = mag(uVec[pid]);
|
||||
|
||||
// Special case: point is close to vertex
|
||||
if (dist[pid] < tol)
|
||||
{
|
||||
weights_[pid] = 1.0;
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
||||
// Project onto unit sphere
|
||||
uVec /= dist;
|
||||
|
||||
|
||||
if (faceIndex < 0)
|
||||
{
|
||||
// Face data not supplied
|
||||
calcWeights
|
||||
(
|
||||
mesh,
|
||||
toGlobal,
|
||||
toLocal,
|
||||
position,
|
||||
uVec,
|
||||
dist,
|
||||
|
||||
weights_
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
DynamicList<point> u(100);
|
||||
const face& f = mesh.faces()[faceIndex];
|
||||
// Collect the uVec for the face
|
||||
forAll(f, j)
|
||||
{
|
||||
u(j) = uVec[toLocal[f[j]]];
|
||||
}
|
||||
|
||||
// Calculate weights for face only
|
||||
calcWeights(toLocal, f, u, dist, weights_);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,162 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
|
||||
\\/ 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/>.
|
||||
|
||||
Class
|
||||
Foam::pointMVCWeight
|
||||
|
||||
Description
|
||||
Container to calculate weights for interpolating directly from vertices
|
||||
of cell using Mean Value Coordinates.
|
||||
|
||||
Based on (VTK's vtkMeanValueCoordinatesInterpolator's) implementation
|
||||
of "Spherical Barycentric Coordinates"
|
||||
2006 paper Eurographics Symposium on Geometry Processing
|
||||
by Torsten Langer, Alexander Belyaev and Hans-Peter Seide
|
||||
|
||||
|
||||
|
||||
SourceFiles
|
||||
pointMVCWeight.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef pointMVCWeight_H
|
||||
#define pointMVCWeight_H
|
||||
|
||||
#include "scalarField.H"
|
||||
#include "vectorField.H"
|
||||
#include "Map.H"
|
||||
#include "DynamicList.H"
|
||||
#include "point.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
class polyMesh;
|
||||
class pointMesh;
|
||||
template<class T> class pointPatchField;
|
||||
template<class Type, template<class> class PatchField, class GeoMesh>
|
||||
class GeometricField;
|
||||
class face;
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class pointMVCWeight Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class pointMVCWeight
|
||||
{
|
||||
protected:
|
||||
|
||||
// Protected data
|
||||
|
||||
//- Cell index
|
||||
const label cellIndex_;
|
||||
|
||||
//- Weights applied to cell vertices
|
||||
scalarField weights_;
|
||||
|
||||
|
||||
// Protected Member Functions
|
||||
|
||||
//- Calculate weights from single face's vertices only
|
||||
void calcWeights
|
||||
(
|
||||
const Map<label>& toLocal,
|
||||
const face& f,
|
||||
const DynamicList<point>& u,
|
||||
const scalarField& dist,
|
||||
scalarField& weights
|
||||
) const;
|
||||
|
||||
//- Calculate weights from all cell's vertices
|
||||
void calcWeights
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const labelList& toGlobal,
|
||||
const Map<label>& toLocal,
|
||||
const vector& position,
|
||||
const vectorField& uVec,
|
||||
const scalarField& dist,
|
||||
scalarField& weights
|
||||
) const;
|
||||
|
||||
public:
|
||||
|
||||
//- Debug switch
|
||||
static int debug;
|
||||
|
||||
//- Tolerance used in calculating barycentric co-ordinates
|
||||
// (applied to normalised values)
|
||||
static scalar tol;
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from components
|
||||
pointMVCWeight
|
||||
(
|
||||
const polyMesh& mesh,
|
||||
const vector& position,
|
||||
const label nCell,
|
||||
const label facei = -1
|
||||
);
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
//- Cell index
|
||||
inline label cell() const
|
||||
{
|
||||
return cellIndex_;
|
||||
}
|
||||
|
||||
//- interpolation weights (in order of cellPoints)
|
||||
inline const scalarField& weights() const
|
||||
{
|
||||
return weights_;
|
||||
}
|
||||
|
||||
//- Interpolate field
|
||||
template<class Type>
|
||||
inline Type interpolate
|
||||
(
|
||||
const GeometricField<Type, pointPatchField, pointMesh>& psip
|
||||
) const;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#include "pointMVCWeightI.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -0,0 +1,48 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
|
||||
\\/ 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 "pointFields.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
template<class Type>
|
||||
inline Type Foam::pointMVCWeight::interpolate
|
||||
(
|
||||
const GeometricField<Type, pointPatchField, pointMesh>& psip
|
||||
) const
|
||||
{
|
||||
const labelList& vertices = psip.mesh()().cellPoints()[cellIndex_];
|
||||
|
||||
Type t = pTraits<Type>::zero;
|
||||
forAll(vertices, i)
|
||||
{
|
||||
t += psip[vertices[i]]*weights_[i];
|
||||
}
|
||||
|
||||
return t;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user