multiple isosurfaces from unregistered samplefield

src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C

@@ -498,8 +498,12 @@ Foam::isoSurface::interpolate
 
 
     // One value per point
-    tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
+    tmp<Field<Type> > tvalues
+    (
+        new Field<Type>(points().size(), pTraits<Type>::zero)
+    );
     Field<Type>& values = tvalues();
+    labelList nValues(values.size(), 0);
 
     forAll(triPoints, i)
     {
@@ -507,8 +511,34 @@ Foam::isoSurface::interpolate
 
         if (mergedPointI >= 0)
         {
-            values[mergedPointI] = triPoints[i];
+            values[mergedPointI] += triPoints[i];
+            nValues[mergedPointI]++;
+        }
+    }
+
+    if (debug)
+    {
+        Pout<< "nValues:" << values.size() << endl;
+        label nMult = 0;
+        forAll(nValues, i)
+        {
+            if (nValues[i] == 0)
+            {
+                FatalErrorIn("isoSurface::interpolate(..)")
+                    << "point:" << i << " nValues:" << nValues[i]
+                    << abort(FatalError);
+            }
+            else if (nValues[i] > 1)
+            {
+                nMult++;
+            }
         }
+        Pout<< "Of which mult:" << nMult << endl;
+    }
+
+    forAll(values, i)
+    {
+        values[i] /= scalar(nValues[i]);
     }
 
     return tvalues;

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C

@@ -65,7 +65,7 @@ void Foam::sampledIsoSurface::getIsoFields() const
         if (debug)
         {
             Info<< "sampledIsoSurface::getIsoField() : reading "
-                << isoField_ << " from time " <<fvm.time().timeName()
+                << isoField_ << " from time " << fvm.time().timeName()
                 << endl;
         }
 
@@ -124,10 +124,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
 
     if (debug)
     {
-        Info<< "sampledIsoSurface::getIsoField() : obtained volField "
+        Info<< "sampledIsoSurface::getIsoField() : volField "
             << volFieldPtr_->name() << " min:" << min(*volFieldPtr_).value()
             << " max:" << max(*volFieldPtr_).value() << endl;
-        Info<< "sampledIsoSurface::getIsoField() : obtained pointField "
+        Info<< "sampledIsoSurface::getIsoField() : pointField "
             << pointFieldPtr_->name()
             << " min:" << gMin(pointFieldPtr_->internalField())
             << " max:" << gMax(pointFieldPtr_->internalField()) << endl;

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C

@@ -57,11 +57,16 @@ Foam::sampledIsoSurface::interpolateField
     const GeometricField<Type, fvPatchField, volMesh>& volFld =
         interpolator.psi();
 
-    tmp<GeometricField<Type, pointPatchField, pointMesh> > pointFld
+    tmp<GeometricField<Type, pointPatchField, pointMesh> > tpointFld
     (
         volPointInterpolation::New(fvm).interpolate(volFld)
     );
 
+    const GeometricField<Type, pointPatchField, pointMesh>& pointFld =
+        tpointFld();
+
+    // Get pointers to sampling field (both original and interpolated one)
+    getIsoFields();
     // Recreate geometry if time has changed
     createGeometry();
 
@@ -71,7 +76,7 @@ Foam::sampledIsoSurface::interpolateField
         *volFieldPtr_,
         *pointFieldPtr_,
         volFld,
-        pointFld()
+        pointFld
     );
 }