From be59f346019c0b86f903b372ded82cba42b79aba Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 8 Jun 2018 17:22:40 +0200
Subject: [PATCH] BUG: error with surfMesh transfer from triSurfaceMesh (closes
 #862)

- face/point ownership is transferred to the surfMesh, so use these
  directly when sampling the interior.
---
 .../sampledTriSurfaceMeshTemplates.C          |  4 +-
 .../triSurfaceMesh/surfMeshSampleDiscrete.C   |  2 +-
 .../triSurfaceMesh/surfMeshSampleDiscrete.H   |  9 ++-
 .../surfMeshSampleDiscreteTemplates.C         | 74 +++++++------------
 .../surface/triSurfaceMesh/discreteSurface.C  |  1 +
 .../triSurfaceMesh/discreteSurfaceTemplates.C | 16 ++--
 6 files changed, 48 insertions(+), 58 deletions(-)

diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
index 48bd03308e..b4f83b3736 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
@@ -58,14 +58,14 @@ Foam::sampledTriSurfaceMesh::sampleOnFaces
     auto& values = tvalues.ref();
 
     const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-
+    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
 
     // Create flat boundary field
-    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
 
     Field<Type> bVals(nBnd, Zero);
 
     const auto& bField = sampler.psi().boundaryField();
+
     forAll(bField, patchi)
     {
         const label bFacei = (pbm[patchi].start() - mesh().nInternalFaces());
diff --git a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.C b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.C
index 2e02fa1982..dce44e0e63 100644
--- a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.C
+++ b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
diff --git a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.H b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.H
index 4dd47e2c89..c7a9b90b9d 100644
--- a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.H
+++ b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscrete.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,13 @@ class surfMeshSampleDiscrete
         //- Transfer mesh content from SurfaceSource to surfMesh
         void transferContent();
 
+        //- Sample field (from volume field) onto surface faces
+        template<class Type>
+        tmp<Field<Type>> sampleOnFaces
+        (
+            const interpolation<Type>& sampler
+        ) const;
+
         //- Sample field on surface.
         template<class Type>
         bool sampleType
diff --git a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscreteTemplates.C b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscreteTemplates.C
index 4d56412539..78cf386c52 100644
--- a/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscreteTemplates.C
+++ b/src/sampling/surfMeshSample/triSurfaceMesh/surfMeshSampleDiscreteTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,11 +24,32 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "surfMeshSampleDiscrete.H"
-#include "dimensionedType.H"
-#include "error.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+Foam::tmp<Foam::Field<Type>>
+Foam::surfMeshSampleDiscrete::sampleOnFaces
+(
+    const interpolation<Type>& sampler
+) const
+{
+    if (onBoundary())
+    {
+        return SurfaceSource::sampleOnFaces(sampler);
+    }
+
+    // Sample cells
+    return surfMeshSample::sampleOnFaces
+    (
+        sampler,
+        sampleElements(),  //< sampleElements == meshCells
+        surface().faces(),
+        surface().points()
+    );
+}
+
+
 template<class Type>
 bool Foam::surfMeshSampleDiscrete::sampleType
 (
@@ -38,7 +59,7 @@ bool Foam::surfMeshSampleDiscrete::sampleType
 {
     typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
 
-    const auto volFldPtr =
+    const auto* volFldPtr =
         SurfaceSource::mesh().lookupObjectPtr<VolFieldType>(fieldName);
 
     if (!volFldPtr)
@@ -48,52 +69,11 @@ bool Foam::surfMeshSampleDiscrete::sampleType
 
     auto samplerPtr = interpolation<Type>::New(sampleScheme, *volFldPtr);
 
-    getOrCreateSurfField<Type>(*volFldPtr).field()
-        = SurfaceSource::sampleOnFaces(*samplerPtr);
+    getOrCreateSurfField<Type>(*volFldPtr).field() =
+        sampleOnFaces(*samplerPtr);
 
     return true;
 }
 
 
-// template<class Type>
-// Foam::tmp<Foam::DimensionedField<Type, Foam::surfGeoMesh>>
-// Foam::surfMeshSampleDiscrete::sampleOnFaces
-// (
-//     const GeometricField<Type, fvPatchField, volMesh>& fld
-// ) const
-// {
-//     typedef DimensionedField<Type, surfGeoMesh> SurfFieldType;
-//
-//     tmp<Field<Type>> tfield = SurfaceSource::sampleOnFaces(fld);
-//     SurfFieldType& result = getOrCreateSurfField<Type>(fld);
-//
-//     // need to verify if this will be needed (in the future)
-//     const surfMesh& s = surface();
-//     if (result.size() != s.size())
-//     {
-//         // maybe resampling changed the surfMesh,
-//         // but the existing surfField wasn't updated
-//
-//         result.resize(s.size(), Zero);
-//     }
-//
-//     if (result.size() != sampleElements().size())
-//     {
-//         FatalErrorInFunction
-//             << "mismatch in field/mesh sizes "
-//             << result.name() << nl
-//             << " field has " << result.size()
-//             << " but sampleElements has " << sampleElements().size() << nl
-//             << endl
-//             << exit(FatalError);
-//
-//         Info<< "WARNING: "
-//             << endl;
-//     }
-//
-//     result = tfield;
-//     return result;
-// }
-
-
 // ************************************************************************* //
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurface.C b/src/sampling/surface/triSurfaceMesh/discreteSurface.C
index ad7ea43131..2c67b787c4 100644
--- a/src/sampling/surface/triSurfaceMesh/discreteSurface.C
+++ b/src/sampling/surface/triSurfaceMesh/discreteSurface.C
@@ -744,6 +744,7 @@ bool Foam::discreteSurface::interpolate() const
     return interpolate_;
 }
 
+
 bool Foam::discreteSurface::interpolate(bool b)
 {
     if (interpolate_ == b)
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C b/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
index a1a7b8a2c9..74256346f4 100644
--- a/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
+++ b/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
@@ -175,24 +175,26 @@ Foam::discreteSurface::sampleOnFaces
 
         Field<Type> bVals(nBnd, Zero);
 
-        forAll(vField.boundaryField(), patchi)
+        const auto& bField = vField.boundaryField();
+
+        forAll(bField, patchi)
         {
-            label bFacei = pbm[patchi].start() - mesh().nInternalFaces();
+            const label bFacei = pbm[patchi].start() - mesh().nInternalFaces();
 
             SubList<Type>
             (
                 bVals,
-                vField.boundaryField()[patchi].size(),
+                bField[patchi].size(),
                 bFacei
-            ) = vField.boundaryField()[patchi];
+            ) = bField[patchi];
         }
 
         // Sample in flat boundary field
 
         for (label i=0; i < len; ++i)
         {
-            label facei = elements[i];
-            values[i] = bVals[facei-mesh().nInternalFaces()];
+            const label bFacei = (elements[i] - mesh().nInternalFaces());
+            values[i] = bVals[bFacei];
         }
     }
 
@@ -230,7 +232,7 @@ Foam::discreteSurface::sampleOnPoints
 
         forAll(samplePoints_, pointi)
         {
-            label facei = sampleElements_[pointi];
+            const label facei = sampleElements_[pointi];
 
             values[pointi] = interpolator.interpolate
             (
-- 
GitLab