Commit f050ca82 authored by mattijs's avatar mattijs
Browse files

read fields

parent 755e8848
......@@ -41,6 +41,128 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::sampledIsoSurface::createGeometry() const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
if (fvm.time().timeIndex() != storedTimeIndex_)
{
storedTimeIndex_ = fvm.time().timeIndex();
// Clear any stored topo
facesPtr_.clear();
// Optionally read volScalarField
autoPtr<volScalarField> readFieldPtr_;
// 1. see if field in database
// 2. see if field can be read
const volScalarField* cellFldPtr = NULL;
if (fvm.foundObject<volScalarField>(isoField_))
{
if (debug)
{
Info<< "sampledIsoSurface::createGeometry() : lookup "
<< isoField_ << endl;
}
cellFldPtr = &fvm.lookupObject<volScalarField>(isoField_);
}
else
{
// Bit of a hack. Read field and store.
if (debug)
{
Info<< "sampledIsoSurface::createGeometry() : reading "
<< isoField_ << " from time " <<fvm.time().timeName()
<< endl;
}
readFieldPtr_.reset
(
new volScalarField
(
IOobject
(
isoField_,
fvm.time().timeName(),
fvm,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
),
fvm
)
);
cellFldPtr = readFieldPtr_.operator->();
}
const volScalarField& cellFld = *cellFldPtr;
tmp<pointScalarField> pointFld
(
volPointInterpolation::New(fvm).interpolate(cellFld)
);
//- Direct from cell field and point field. Gives bad continuity.
//const isoSurface iso
//(
// fvm,
// cellFld.internalField(),
// pointFld().internalField(),
// isoVal_
//);
//- From point field and interpolated cell.
scalarField cellAvg(fvm.nCells(), scalar(0.0));
labelField nPointCells(fvm.nCells(), 0);
{
for (label pointI = 0; pointI < fvm.nPoints(); pointI++)
{
const labelList& pCells = fvm.pointCells(pointI);
forAll(pCells, i)
{
label cellI = pCells[i];
cellAvg[cellI] += pointFld().internalField()[pointI];
nPointCells[cellI]++;
}
}
}
forAll(cellAvg, cellI)
{
cellAvg[cellI] /= nPointCells[cellI];
}
const isoSurface iso
(
fvm,
cellAvg,
pointFld().internalField(),
isoVal_
);
const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
meshCells_ = iso.meshCells();
triPointMergeMap_ = iso.triPointMergeMap();
if (debug)
{
Pout<< "sampledIsoSurface::createGeometry() : constructed iso:"
<< nl
<< " isoField : " << isoField_ << nl
<< " isoValue : " << isoVal_ << nl
<< " unmerged points: " << triPointMergeMap_.size() << nl
<< " points : " << points().size() << nl
<< " tris : " << triSurface::size() << nl
<< " cut cells : " << meshCells_.size() << endl;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......
......@@ -84,6 +84,9 @@ class sampledIsoSurface
// Private Member Functions
//- Create iso surface (if time has changed)
void createGeometry() const;
//- sample field on faces
template <class Type>
tmp<Field<Type> > sampleField
......
......@@ -39,35 +39,8 @@ Foam::sampledIsoSurface::sampleField
const GeometricField<Type, fvPatchField, volMesh>& vField
) const
{
const fvMesh& fvm = vField.mesh();
if (fvm.time().timeIndex() != storedTimeIndex_)
{
storedTimeIndex_ = fvm.time().timeIndex();
//- Clear any stored topo
facesPtr_.clear();
const volScalarField& cellFld =
fvm.lookupObject<volScalarField>(isoField_);
tmp<pointScalarField> pointFld
(
volPointInterpolation::New(fvm).interpolate(cellFld)
);
const isoSurface iso
(
fvm,
cellFld.internalField(),
pointFld().internalField(),
isoVal_
);
const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
meshCells_ = iso.meshCells();
triPointMergeMap_ = iso.triPointMergeMap();
}
// Recreate geometry if time has changed
createGeometry();
return tmp<Field<Type> >(new Field<Type>(vField, meshCells_));
}
......@@ -80,36 +53,8 @@ Foam::sampledIsoSurface::interpolateField
const interpolation<Type>& interpolator
) const
{
const fvMesh& fvm = static_cast<const fvMesh&>(mesh());
if (fvm.time().timeIndex() != storedTimeIndex_)
{
//- Clear any stored topo
facesPtr_.clear();
storedTimeIndex_ = fvm.time().timeIndex();
const volScalarField& cellFld =
fvm.lookupObject<volScalarField>(isoField_);
tmp<pointScalarField> pointFld
(
volPointInterpolation::New(fvm).interpolate(cellFld)
);
const isoSurface iso
(
fvm,
cellFld.internalField(),
pointFld().internalField(),
isoVal_
);
const_cast<sampledIsoSurface&>(*this).triSurface::operator=(iso);
meshCells_ = iso.meshCells();
triPointMergeMap_ = iso.triPointMergeMap();
}
// Recreate geometry if time has changed
createGeometry();
// One value per point
tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment