iso surfaces on coupled patches

src/sampling/sampledSurface/distanceSurface/distanceSurfaceTemplates.C

@@ -62,13 +62,7 @@ Foam::distanceSurface::interpolateField
     );
 
     // Sample.
-    return surface().interpolate
-    (
-        cellDistancePtr_(),
-        pointDistance_,
-        volFld,
-        pointFld()
-    );
+    return surface().interpolate(volFld, pointFld());
 }
 
 

src/sampling/sampledSurface/isoSurface/isoSurface.H

@@ -56,6 +56,7 @@ SourceFiles
 #include "pointIndexHit.H"
 #include "PackedBoolList.H"
 #include "volFields.H"
+#include "slicedVolFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -92,6 +93,11 @@ class isoSurface
         //- Reference to mesh
         const fvMesh& mesh_;
 
+        const scalarField& pVals_;
+
+        //- Input volScalarField with separated coupled patches rewritten
+        autoPtr<slicedVolScalarField> cValsPtr_;
+
         //- Isosurface value
         const scalar iso_;
 
@@ -117,6 +123,38 @@ class isoSurface
 
     // Private Member Functions
 
+        // Point synchronisation
+
+            //- Does tensor differ (to within mergeTolerance) from identity
+            bool noTransform(const tensor& tt) const;
+
+            //- Is patch a collocated (i.e. non-separated) coupled patch?
+            static bool collocatedPatch(const polyPatch&);
+
+            //- Per face whether is collocated
+            PackedBoolList collocatedFaces(const coupledPolyPatch&) const;
+
+            //- Take value at local point patchPointI and assign it to its
+            //  correct place in patchValues (for transferral) and sharedValues
+            //  (for reduction)
+            void insertPointData
+            (
+                const processorPolyPatch& pp,
+                const Map<label>& meshToShared,
+                const pointField& pointValues,
+                const label patchPointI,
+                pointField& patchValues,
+                pointField& sharedValues
+            ) const;
+
+            //- Synchonise points on all non-separated coupled patches
+            void syncUnseparatedPoints
+            (
+                pointField& collapsedPoint,
+                const point& nullValue
+            ) const;
+
+
         //- Get location of iso value as fraction inbetween s0,s1
         scalar isoFraction
         (
@@ -136,6 +174,7 @@ class isoSurface
         void getNeighbour
         (
             const labelList& boundaryRegion,
+            const volVectorField& meshC,
             const volScalarField& cVals,
             const label cellI,
             const label faceI,
@@ -147,6 +186,7 @@ class isoSurface
         void calcCutTypes
         (
             const labelList& boundaryRegion,
+            const volVectorField& meshC,
             const volScalarField& cVals,
             const scalarField& pVals
         );
@@ -170,6 +210,7 @@ class isoSurface
         void calcSnappedCc
         (
             const labelList& boundaryRegion,
+            const volVectorField& meshC,
             const volScalarField& cVals,
             const scalarField& pVals,
             DynamicList<point>& snappedPoints,
@@ -182,12 +223,23 @@ class isoSurface
         (
             const PackedBoolList& isBoundaryPoint,
             const labelList& boundaryRegion,
+            const volVectorField& meshC,
             const volScalarField& cVals,
             const scalarField& pVals,
             DynamicList<point>& snappedPoints,
             labelList& snappedPoint
         ) const;
 
+
+        //- Return input field with coupled (and empty) patch values rewritten
+        template<class Type>
+        tmp<SlicedGeometricField
+        <Type, fvPatchField, slicedFvPatchField, volMesh> >
+        adaptPatchFields
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& fld
+        ) const;
+
         //- Generate single point by interpolation or snapping
         template<class Type>
         Type generatePoint
@@ -345,8 +397,8 @@ public:
     // Constructors
 
         //- Construct from cell values and point values. Uses boundaryField
-        //  for boundary values. Requires on coupled patchfields to be set
-        //  to the opposite cell value.
+        //  for boundary values. Holds reference to cellIsoVals and
+        //  pointIsoVals.
         isoSurface
         (
             const volScalarField& cellIsoVals,
@@ -376,8 +428,6 @@ public:
         template <class Type>
         tmp<Field<Type> > interpolate
         (
-            const volScalarField& cellIsoVals,
-            const scalarField& pointIsoVals,
             const GeometricField<Type, fvPatchField, volMesh>& cCoords,
             const Field<Type>& pCoords
         ) const;

src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C

@@ -27,9 +27,114 @@ License
 #include "isoSurface.H"
 #include "polyMesh.H"
 #include "syncTools.H"
+#include "surfaceFields.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+template<class Type>
+Foam::tmp<Foam::SlicedGeometricField
+<
+    Type,
+    Foam::fvPatchField,
+    Foam::slicedFvPatchField,
+    Foam::volMesh
+> >
+Foam::isoSurface::adaptPatchFields
+(
+    const GeometricField<Type, fvPatchField, volMesh>& fld
+) const
+{
+    typedef SlicedGeometricField
+    <
+        Type,
+        fvPatchField,
+        slicedFvPatchField,
+        volMesh
+    > FieldType;
+
+    tmp<FieldType> tsliceFld
+    (
+        new FieldType
+        (
+            IOobject
+            (
+                fld.name(),
+                fld.instance(),
+                fld.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fld,        // internal field
+            true        // preserveCouples
+        )
+    );
+    FieldType& sliceFld = tsliceFld();
+
+    const fvMesh& mesh = fld.mesh();
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (isA<emptyPolyPatch>(pp))
+        {
+            // Clear old value. Cannot resize it since is a slice.
+            sliceFld.boundaryField().set(patchI, NULL);
+
+            // Set new value we can change
+            sliceFld.boundaryField().set
+            (
+                patchI,
+                new calculatedFvPatchField<Type>
+                (
+                    mesh.boundary()[patchI],
+                    sliceFld
+                )
+            );
+            sliceFld.boundaryField()[patchI] ==
+                mesh.boundary()[patchI].patchInternalField
+                (
+                    sliceFld
+                );
+        }
+        else if (isA<cyclicPolyPatch>(pp))
+        {
+            // Already has interpolate as value
+        }
+        else if (isA<processorPolyPatch>(pp) && !collocatedPatch(pp))
+        {
+            fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
+            (
+                fld.boundaryField()[patchI]
+            );
+
+            const scalarField& w = mesh.weights().boundaryField()[patchI];
+
+            tmp<Field<Type> > f =
+                w*pfld.patchInternalField()
+              + (1.0-w)*pfld.patchNeighbourField();
+
+            PackedBoolList isCollocated
+            (
+                collocatedFaces(refCast<const processorPolyPatch>(pp))
+            );
+
+            forAll(isCollocated, i)
+            {
+                if (!isCollocated[i])
+                {
+                    pfld[i] = f()[i];
+                }
+            }
+        }
+    }
+    return tsliceFld;
+}
+
+
 template<class Type>
 Type Foam::isoSurface::generatePoint
 (
@@ -389,7 +494,6 @@ void Foam::isoSurface::generateTriPoints
     }
 
 
-
     // Generate triangle points
 
     triPoints.clear();
@@ -456,16 +560,16 @@ void Foam::isoSurface::generateTriPoints
     syncTools::swapBoundaryFaceList(mesh_, neiSnappedPoint, false);
 
 
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
 
-        if
-        (
-            isA<processorPolyPatch>(pp)
-         && !refCast<const processorPolyPatch>(pp).separated()
-        )
+        if (isA<processorPolyPatch>(pp) && collocatedPatch(pp))
         {
+            // Coincident processor patches. Boundary field of 
+            // cVals and cCoords is opposite cell.
+
             //if (refCast<const processorPolyPatch>(pp).owner())
             {
                 label faceI = pp.start();
@@ -474,18 +578,6 @@ void Foam::isoSurface::generateTriPoints
                 {
                     if (faceCutType_[faceI] != NOTCUT)
                     {
-                        label bFaceI = faceI-mesh_.nInternalFaces();
-                        if
-                        (
-                            neiSnapped[bFaceI]
-                         && (neiSnappedPoint[bFaceI]==pTraits<Type>::zero)
-                        )
-                        {
-                            FatalErrorIn("isoSurface::generateTriPoints(..)")
-                                << "problem:" << abort(FatalError);
-                        }
-
-
                         generateFaceTriPoints
                         (
                             cVals,
@@ -512,84 +604,6 @@ void Foam::isoSurface::generateTriPoints
                 }
             }
         }
-        else if (isA<emptyPolyPatch>(pp))
-        {
-            // Check if field is there (when generating geometry the
-            // empty patches have been rewritten to be the face centres),
-            // otherwise use zero-gradient.
-
-            label faceI = pp.start();
-
-
-            const fvPatchScalarField& fvp = cVals.boundaryField()[patchI];
-
-            // Owner value of cVals
-            scalarField internalVals;
-            if (fvp.size() == 0)
-            {
-                internalVals.setSize(pp.size());
-                forAll(pp, i)
-                {
-                    internalVals[i] = cVals[own[pp.start()+i]];
-                }
-            }
-            const scalarField& bVals =
-            (
-                fvp.size() > 0
-              ? static_cast<const scalarField&>(fvp)
-              : internalVals
-            );
-
-
-            const fvPatchField<Type>& pc = cCoords.boundaryField()[patchI];
-
-            // Owner value of cCoords
-            Field<Type> internalCoords;
-            if (pc.size() == 0)
-            {
-                internalCoords.setSize(pp.size());
-                forAll(pp, i)
-                {
-                    internalCoords[i] = cCoords[own[pp.start()+i]];
-                }
-            }
-            const Field<Type>& bCoords =
-            (
-                pc.size() > 0
-              ? static_cast<const Field<Type>&>(pc)
-              : internalCoords
-            );
-
-
-            forAll(pp, i)
-            {
-                if (faceCutType_[faceI] != NOTCUT)
-                {
-                    generateFaceTriPoints
-                    (
-                        cVals,
-                        pVals,
-
-                        cCoords,
-                        pCoords,
-
-                        snappedPoints,
-                        snappedCc,
-                        snappedPoint,
-                        faceI,
-
-                        bVals[i],
-                        bCoords[i],
-                        false,  // fc not snapped
-                        pTraits<Type>::zero,
-
-                        triPoints,
-                        triMeshCells
-                    );
-                }
-                faceI++;
-            }
-        }
         else
         {
             label faceI = pp.start();
@@ -642,12 +656,20 @@ template <class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::isoSurface::interpolate
 (
-    const volScalarField& cVals,
-    const scalarField& pVals,
     const GeometricField<Type, fvPatchField, volMesh>& cCoords,
     const Field<Type>& pCoords
 ) const
 {
+    // Recalculate boundary values
+    tmp<SlicedGeometricField
+    <
+        Type,
+        fvPatchField,
+        slicedFvPatchField,
+        volMesh
+    > > c2(adaptPatchFields(cCoords));
+
+
     DynamicList<Type> triPoints(nCutCells_);
     DynamicList<label> triMeshCells(nCutCells_);
 
@@ -658,10 +680,10 @@ Foam::isoSurface::interpolate
 
     generateTriPoints
     (
-        cVals,
-        pVals,
+        cValsPtr_(),
+        pVals_,
 
-        cCoords,
+        c2(),
         pCoords,
 
         snappedPoints,

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceTemplates.C

@@ -71,13 +71,7 @@ Foam::sampledIsoSurface::interpolateField
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
         // Sample.
-        return surface().interpolate
-        (
-            *volSubFieldPtr_,
-            *pointSubFieldPtr_,
-            volSubFld,
-            tpointSubFld()
-        );
+        return surface().interpolate(volSubFld, tpointSubFld());
     }
     else
     {
@@ -85,13 +79,7 @@ Foam::sampledIsoSurface::interpolateField
             volPointInterpolation::New(volFld.mesh()).interpolate(volFld);
 
         // Sample.
-        return surface().interpolate
-        (
-            *volFieldPtr_,
-            *pointFieldPtr_,
-            volFld,
-            tpointFld()
-        );
+        return surface().interpolate(volFld, tpointFld());
     }
 }
 

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlaneTemplates.C

@@ -67,13 +67,7 @@ Foam::sampledCuttingPlane::interpolateField
             volPointInterpolation::New(volSubFld.mesh()).interpolate(volSubFld);
 
         // Sample.
-        return surface().interpolate
-        (
-            cellDistancePtr_(),
-            pointDistance_,
-            volSubFld,
-            tpointSubFld()
-        );
+        return surface().interpolate(volSubFld, tpointSubFld());
     }
     else
     {
@@ -83,13 +77,7 @@ Foam::sampledCuttingPlane::interpolateField
         );
 
         // Sample.
-        return surface().interpolate
-        (
-            cellDistancePtr_(),
-            pointDistance_,
-            volFld,
-            tpointFld()
-        );
+        return surface().interpolate(volFld, tpointFld());
     }
 }