From f050ca8248001dcbcb112779e074516e07c0d0a6 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 21 Oct 2008 19:52:10 +0100
Subject: [PATCH] read fields

---
 .../isoSurface/sampledIsoSurface.C            | 122 ++++++++++++++++++
 .../isoSurface/sampledIsoSurface.H            |   3 +
 .../isoSurface/sampledIsoSurfaceTemplates.C   |  63 +--------
 3 files changed, 129 insertions(+), 59 deletions(-)

diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index 53344e6282..a5ad5dc0ed 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -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  * * * * * * * * * * * * * * //
 
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index 57445082f6..a55cb6641a 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -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
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
index 755fc95939..dc2abbd74b 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C
@@ -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()));