Commit cdb0a8a8 authored by mattijs's avatar mattijs
Browse files

iso surface

parent 204f93e8
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "isoSurface.H"
#include "dictionary.H"
#include "polyMesh.H"
#include "volFields.H"
#include "mergePoints.H"
#include "tetMatcher.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(isoSurface, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::isoSurface::isoSurface
(
const polyMesh& mesh,
const scalarField& cellValues,
const scalarField& pointValues,
const scalar iso
)
:
mesh_(mesh),
cellValues_(cellValues),
pointValues_(pointValues),
iso_(iso)
{
const pointField& cellCentres = mesh.cellCentres();
tetMatcher tet;
DynamicList<point> triPoints;
DynamicList<label> triMeshCells;
forAll(mesh.cells(), cellI)
{
label oldNPoints = triPoints.size();
const cell& cFaces = mesh.cells()[cellI];
if (tet.isA(mesh, cellI))
{
// For tets don't do cell-centre decomposition, just use the
// tet points and values
const face& f0 = mesh.faces()[cFaces[0]];
// Get the other point
const face& f1 = mesh.faces()[cFaces[1]];
label oppositeI = -1;
forAll(f1, fp)
{
oppositeI = f1[fp];
if (findIndex(f0, oppositeI) == -1)
{
break;
}
}
vertexInterp
(
iso,
pointValues[f0[0]],
pointValues[f0[1]],
pointValues[f0[2]],
pointValues[oppositeI],
mesh.points()[f0[0]],
mesh.points()[f0[1]],
mesh.points()[f0[2]],
mesh.points()[oppositeI],
triPoints
);
}
else
{
const cell& cFaces = mesh.cells()[cellI];
forAll(cFaces, cFaceI)
{
label faceI = cFaces[cFaceI];
const face& f = mesh.faces()[faceI];
// Do a tetrahedrisation. Each face to cc becomes pyr.
// Each pyr gets split into two tets by diagionalisation
// of face. So
// - f[0], f[1], f[2], cc
// - f[0], f[2], f[3], cc
for(label fp = 1; fp < f.size() - 1; fp++)
{
vertexInterp
(
iso,
pointValues[f[0]],
pointValues[f[fp]],
pointValues[f[f.fcIndex(fp)]],
cellValues[cellI],
mesh.points()[f[0]],
mesh.points()[f[fp]],
mesh.points()[f[f.fcIndex(fp)]],
cellCentres[cellI],
triPoints
);
}
}
}
// Every three triPoints is a cell
label nCells = (triPoints.size()-oldNPoints)/3;
for (label i = 0; i < nCells; i++)
{
triMeshCells.append(cellI);
}
}
triPoints.shrink();
triMeshCells.shrink();
meshCells_.transfer(triMeshCells);
pointField newPoints;
mergePoints(triPoints, SMALL, false, triPointMergeMap_, newPoints);
DynamicList<labelledTri> tris(meshCells_.size());
forAll(meshCells_, triI)
{
tris.append
(
labelledTri
(
triPointMergeMap_[3*triI],
triPointMergeMap_[3*triI+1],
triPointMergeMap_[3*triI+2],
0
)
);
}
//- 1.5.x:
//tris.shrink();
//triSurface::operator=
//(
// triSurface(tris, geometricSurfacePatchList(0), newPoints)
//);
triSurface::operator=
(
triSurface(tris, geometricSurfacePatchList(0), newPoints, true)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
//Foam::tmp<Foam::scalarField>
//Foam::isoSurface::sample
//(
// const volScalarField& vField
//) const
//{
// return sampleField(vField);
//}
//
//
//Foam::tmp<Foam::vectorField>
//Foam::isoSurface::sample
//(
// const volVectorField& vField
//) const
//{
// return sampleField(vField);
//}
//
//Foam::tmp<Foam::sphericalTensorField>
//Foam::isoSurface::sample
//(
// const volSphericalTensorField& vField
//) const
//{
// return sampleField(vField);
//}
//
//
//Foam::tmp<Foam::symmTensorField>
//Foam::isoSurface::sample
//(
// const volSymmTensorField& vField
//) const
//{
// return sampleField(vField);
//}
//
//
//Foam::tmp<Foam::tensorField>
//Foam::isoSurface::sample
//(
// const volTensorField& vField
//) const
//{
// return sampleField(vField);
//}
//
//
//Foam::tmp<Foam::scalarField>
//Foam::isoSurface::interpolate
//(
// const interpolation<scalar>& interpolator
//) const
//{
// return interpolateField(interpolator);
//}
//
//
//Foam::tmp<Foam::vectorField>
//Foam::isoSurface::interpolate
//(
// const interpolation<vector>& interpolator
//) const
//{
// return interpolateField(interpolator);
//}
//
//Foam::tmp<Foam::sphericalTensorField>
//Foam::isoSurface::interpolate
//(
// const interpolation<sphericalTensor>& interpolator
//) const
//{
// return interpolateField(interpolator);
//}
//
//
//Foam::tmp<Foam::symmTensorField>
//Foam::isoSurface::interpolate
//(
// const interpolation<symmTensor>& interpolator
//) const
//{
// return interpolateField(interpolator);
//}
//
//
//Foam::tmp<Foam::tensorField>
//Foam::isoSurface::interpolate
//(
// const interpolation<tensor>& interpolator
//) const
//{
// return interpolateField(interpolator);
//}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
Foam::isoSurface
Description
A surface formed by the iso value.
After "Polygonising A Scalar Field Using Tetrahedrons", Paul Bourke
SourceFiles
isoSurface.C
\*---------------------------------------------------------------------------*/
#ifndef isoSurface_H
#define isoSurface_H
#include "triSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class polyMesh;
/*---------------------------------------------------------------------------*\
Class isoSurface Declaration
\*---------------------------------------------------------------------------*/
class isoSurface
:
public triSurface
{
// Private data
//- Reference to mesh
const polyMesh& mesh_;
//- Reference to cell values
const scalarField& cellValues_;
//- Reference to point values
const scalarField& pointValues_;
//- Isosurface value
const scalar iso_;
//- For every triangle the original cell in mesh
labelList meshCells_;
//- For every unmerged triangle point the point in the triSurface
labelList triPointMergeMap_;
// Private Member Functions
template<class T>
static T vertexInterp
(
const scalar iso,
const T& p0,
const T& p1,
const scalar s0,
const scalar s1
);
template<class T>
static void vertexInterp
(
const scalar iso,
const scalar s0,
const scalar s1,
const scalar s2,
const scalar s3,
const T& p0,
const T& p1,
const T& p2,
const T& p3,
DynamicList<T>& points
);
public:
//- Runtime type information
TypeName("isoSurface");
// Constructors
//- Construct from dictionary
isoSurface
(
const polyMesh& mesh,
const scalarField& cellValues,
const scalarField& pointValues,
const scalar iso
);
// Member Functions
//- For every face original cell in mesh
const labelList& meshCells() const
{
return meshCells_;
}
//- For every unmerged triangle point the point in the triSurface
const labelList triPointMergeMap() const
{
return triPointMergeMap_;
}
//- sample field on faces
template <class Type>
tmp<Field<Type> > sample
(
const Field<Type>& sampleCellValues
) const;
//- interpolate field to points
template <class Type>
tmp<Field<Type> >
interpolate
(
const Field<Type>& sampleCellValues,
const Field<Type>& samplePointValues
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "isoSurfaceTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "isoSurface.H"
#include "polyMesh.H"
#include "tetMatcher.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
Type Foam::isoSurface::vertexInterp
(
const scalar iso,
const Type& p0,
const Type& p1,
const scalar s0,
const scalar s1
)
{
scalar d = s1-s0;
if (mag(d) > VSMALL)
{
return (iso-s0)/d*p1 + (s1-iso)/d*p0;
}
else
{
return 0.5*(p0+p1);
}
}
// After "Polygonising A Scalar Field Using Tetrahedrons"
// by Paul Bourke
// Get value consistent with uncompacted triangle points.
// Given tet corner sample values s0..s3 interpolate the corresponding
// values p0..p3 to construct the surface corresponding to sample value iso.
template<class Type>
void Foam::isoSurface::vertexInterp
(
const scalar iso,
const scalar s0,
const scalar s1,
const scalar s2,
const scalar s3,
const Type& p0,
const Type& p1,
const Type& p2,
const Type& p3,
DynamicList<Type>& points
)
{
int triIndex = 0;
if (s0 < iso)
{
triIndex |= 1;
}
if (s1 < iso)
{
triIndex |= 2;
}
if (s2 < iso)
{
triIndex |= 4;
}
if (s3 < iso)
{
triIndex |= 8;
}
/* Form the vertices of the triangles for each case */
switch (triIndex)
{
case 0x00:
case 0x0F:
break;
case 0x0E:
case 0x01:
points.append(vertexInterp(iso,p0,p1,s0,s1));
points.append(vertexInterp(iso,p0,p2,s0,s2));
points.append(vertexInterp(iso,p0,p3,s0,s3));
break;
case 0x0D:
case 0x02:
points.append(vertexInterp(iso,p1,p0,s1,s0));
points.append(vertexInterp(iso,p1,p3,s1,s3));
points.append(vertexInterp(iso,p1,p2,s1,s2));
break;
case 0x0C:
case 0x03:
{
const Type tp1 = vertexInterp(iso,p0,p2,s0,s2);
const Type tp2 = vertexInterp(iso,p1,p3,s1,s3);
points.append(vertexInterp(iso,p0,p3,s0,s3));
points.append(tp1);
points.append(tp2);
points.append(tp2);
points.append(vertexInterp(iso,p1,p2,s1,s2));
points.append(tp1);
}
break;
case 0x0B:
case 0x04:
{
points.append(vertexInterp(iso,p2,p0,s2,s0));
points.append(vertexInterp(iso,p2,p1,s2,s1));
points.append(vertexInterp(iso,p2,p3,s2,s3));
}
break;
case 0x0A:
case 0x05:
{
const Type tp0 = vertexInterp(iso,p0,p1,s0,s1);
const Type tp1 = vertexInterp(iso,p2,p3,s2,s3);
points.append(tp0);
points.append(tp1);
points.append(vertexInterp(iso,p0,p3,s0,s3));
points.append(tp0);