diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index 8c4acea9967b39942148f148308d5ae83c0af641..f337d14253dd69b8e361200083f92899ef8530f5 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -30,7 +30,8 @@ License
#include "volPointInterpolation.H"
#include "addToRunTimeSelectionTable.H"
#include "fvMesh.H"
-#include "isoSurfaceCell.H"
+#include "isoSurface.H"
+//#include "isoSurfaceCell.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -42,19 +43,46 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-void Foam::distanceSurface::createGeometry() const
+void Foam::distanceSurface::createGeometry()
{
// Clear any stored topo
facesPtr_.clear();
+
+ const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
// Distance to cell centres
- scalarField cellDistance(mesh().nCells());
+ // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+ cellDistancePtr_.reset
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ "cellDistance",
+ fvm.time().timeName(),
+ fvm.time(),
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ fvm,
+ dimensionedScalar("zero", dimless/dimTime, 0)
+ )
+ );
+ volScalarField& cellDistance = cellDistancePtr_();
+
+ // Internal field
{
+ const pointField& cc = fvm.C();
+ scalarField& fld = cellDistance.internalField();
+
List<pointIndexHit> nearest;
surfPtr_().findNearest
(
- mesh().cellCentres(),
- scalarField(mesh().nCells(), GREAT),
+ cc,
+ scalarField(cc.size(), GREAT),
nearest
);
@@ -63,98 +91,109 @@ void Foam::distanceSurface::createGeometry() const
vectorField normal;
surfPtr_().getNormal(nearest, normal);
- forAll(cellDistance, cellI)
+ forAll(nearest, i)
{
- vector d(mesh().cellCentres()[cellI]-nearest[cellI].hitPoint());
+ vector d(cc[i]-nearest[i].hitPoint());
- if ((d&normal[cellI]) > 0)
+ if ((d&normal[i]) > 0)
{
- cellDistance[cellI] = Foam::mag(d);
+ fld[i] = Foam::mag(d);
}
else
{
- cellDistance[cellI] = -Foam::mag(d);
+ fld[i] = -Foam::mag(d);
}
}
}
else
{
- forAll(cellDistance, cellI)
+ forAll(nearest, i)
{
- cellDistance[cellI] = Foam::mag
- (
- nearest[cellI].hitPoint()
- - mesh().cellCentres()[cellI]
- );
+ fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
}
}
}
+ // Patch fields
+ {
+ forAll(fvm.C().boundaryField(), patchI)
+ {
+ const pointField& cc = fvm.C().boundaryField()[patchI];
+ fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
+
+ List<pointIndexHit> nearest;
+ surfPtr_().findNearest
+ (
+ cc,
+ scalarField(cc.size(), GREAT),
+ nearest
+ );
+
+ if (signed_)
+ {
+ vectorField normal;
+ surfPtr_().getNormal(nearest, normal);
+
+ forAll(nearest, i)
+ {
+ vector d(cc[i]-nearest[i].hitPoint());
+
+ if ((d&normal[i]) > 0)
+ {
+ fld[i] = Foam::mag(d);
+ }
+ else
+ {
+ fld[i] = -Foam::mag(d);
+ }
+ }
+ }
+ else
+ {
+ forAll(nearest, i)
+ {
+ fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
+ }
+ }
+ }
+ }
+
+ // On processor patches the mesh.C() will already be the cell centre
+ // on the opposite side so no need to swap cellDistance.
+
+
// Distance to points
- scalarField pointDistance(mesh().nPoints());
+ pointDistance_.setSize(fvm.nPoints());
{
List<pointIndexHit> nearest;
surfPtr_().findNearest
(
- mesh().points(),
- scalarField(mesh().nPoints(), GREAT),
+ fvm.points(),
+ scalarField(fvm.nPoints(), GREAT),
nearest
);
- forAll(pointDistance, pointI)
+ forAll(pointDistance_, pointI)
{
- pointDistance[pointI] = Foam::mag
+ pointDistance_[pointI] = Foam::mag
(
nearest[pointI].hitPoint()
- - mesh().points()[pointI]
+ - fvm.points()[pointI]
);
}
}
//- Direct from cell field and point field.
- const isoSurfaceCell iso
+ isoSurfPtr_.reset
(
- mesh(),
- cellDistance,
- pointDistance,
- distance_,
- regularise_
+ new isoSurface
+ (
+ cellDistance,
+ pointDistance_,
+ distance_,
+ regularise_
+ )
);
- ////- From point field and interpolated cell.
- //scalarField cellAvg(mesh().nCells(), scalar(0.0));
- //labelField nPointCells(mesh().nCells(), 0);
- //{
- // for (label pointI = 0; pointI < mesh().nPoints(); pointI++)
- // {
- // const labelList& pCells = mesh().pointCells(pointI);
- //
- // forAll(pCells, i)
- // {
- // label cellI = pCells[i];
- //
- // cellAvg[cellI] += pointDistance[pointI];
- // nPointCells[cellI]++;
- // }
- // }
- //}
- //forAll(cellAvg, cellI)
- //{
- // cellAvg[cellI] /= nPointCells[cellI];
- //}
- //
- //const isoSurface iso
- //(
- // mesh(),
- // cellAvg,
- // pointDistance,
- // distance_,
- // regularise_
- //);
-
-
- const_cast<distanceSurface&>(*this).triSurface::operator=(iso);
- meshCells_ = iso.meshCells();
-
if (debug)
{
print(Pout);
@@ -193,9 +232,8 @@ Foam::distanceSurface::distanceSurface
signed_(readBool(dict.lookup("signed"))),
regularise_(dict.lookupOrDefault("regularise", true)),
zoneName_(word::null),
- facesPtr_(NULL),
- storedTimeIndex_(-1),
- meshCells_(0)
+ isoSurfPtr_(NULL),
+ facesPtr_(NULL)
{
// label zoneId = -1;
// if (dict.readIfPresent("zone", zoneName_))
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
index 82fba8a83f3ff696122ffb23102d4838cb6f6c4d..64aebb1a9a0a588c9201dac2261b2727dc554871 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
@@ -37,8 +37,9 @@ SourceFiles
#define distanceSurface_H
#include "sampledSurface.H"
-#include "triSurface.H"
#include "searchableSurface.H"
+//#include "isoSurfaceCell.H"
+#include "isoSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -51,8 +52,7 @@ namespace Foam
class distanceSurface
:
- public sampledSurface,
- public triSurface
+ public sampledSurface
{
// Private data
@@ -71,23 +71,24 @@ class distanceSurface
//- zone name (if restricted to zones)
word zoneName_;
- //- triangles converted to faceList
- mutable autoPtr<faceList> facesPtr_;
+ //- Distance to cell centres
+ autoPtr<volScalarField> cellDistancePtr_;
- // Recreated for every isoSurface
+ //- Distance to points
+ scalarField pointDistance_;
- //- Time at last call
- mutable label storedTimeIndex_;
+ //- Constructed iso surface
+ autoPtr<isoSurface> isoSurfPtr_;
- //- For every triangle the original cell in mesh
- mutable labelList meshCells_;
+ //- triangles converted to faceList
+ mutable autoPtr<faceList> facesPtr_;
// Private Member Functions
- //- Create iso surface (if time has changed)
- void createGeometry() const;
+ //- Create iso surface
+ void createGeometry();
//- sample field on faces
template <class Type>
@@ -129,7 +130,7 @@ public:
//- Points of surface
virtual const pointField& points() const
{
- return triSurface::points();
+ return surface().points();
}
//- Faces of surface
@@ -137,7 +138,7 @@ public:
{
if (!facesPtr_.valid())
{
- const triSurface& s = *this;
+ const triSurface& s = surface();
facesPtr_.reset(new faceList(s.size()));
@@ -153,6 +154,10 @@ public:
//- Correct for mesh movement and/or field changes
virtual void correct(const bool meshChanged);
+ const isoSurface& surface() const
+ {
+ return isoSurfPtr_();
+ }
//- sample field on surface
virtual tmp<scalarField> sample
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
index ce983c7f5f08bffe29eda422a0f47dbff83d08cb..fd1447af78ceae51b06ee26529dfdc8deb1f2734 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C
@@ -39,7 +39,7 @@ Foam::distanceSurface::sampleField
const GeometricField<Type, fvPatchField, volMesh>& vField
) const
{
- return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
+ return tmp<Field<Type> >(new Field<Type>(vField, surface().meshCells()));
}
@@ -50,33 +50,25 @@ Foam::distanceSurface::interpolateField
const interpolation<Type>& interpolator
) const
{
- // One value per point
- tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
- Field<Type>& values = tvalues();
-
- boolList pointDone(points().size(), false);
-
- forAll(faces(), cutFaceI)
- {
- const face& f = faces()[cutFaceI];
-
- forAll(f, faceVertI)
- {
- label pointI = f[faceVertI];
-
- if (!pointDone[pointI])
- {
- values[pointI] = interpolator.interpolate
- (
- points()[pointI],
- meshCells_[cutFaceI]
- );
- pointDone[pointI] = true;
- }
- }
- }
-
- return tvalues;
+ const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
+
+ // Get fields to sample. Assume volPointInterpolation!
+ const GeometricField<Type, fvPatchField, volMesh>& volFld =
+ interpolator.psi();
+
+ tmp<GeometricField<Type, pointPatchField, pointMesh> > pointFld
+ (
+ volPointInterpolation::New(fvm).interpolate(volFld)
+ );
+
+ // Sample.
+ return surface().interpolate
+ (
+ cellDistancePtr_(),
+ pointDistance_,
+ volFld,
+ pointFld()
+ );
}