Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

src/Allwmake

@@ -56,9 +56,7 @@ parallel/Allwmake $targetType $*
 randomProcesses/Allwmake $targetType $*
 
 wmake $targetType overset
-
 wmake $targetType ODE
-wmake $targetType fvMotionSolver
 
 transportModels/Allwmake $targetType $*
 thermophysicalModels/Allwmake $targetType $*
@@ -66,13 +64,16 @@ TurbulenceModels/Allwmake $targetType $*
 wmake $targetType combustionModels
 regionModels/Allwmake $targetType $*
 lagrangian/Allwmake $targetType $*
+
+wmake $targetType fvOptions
+wmake $targetType fvMotionSolver
+
 mesh/Allwmake $targetType $*
 renumber/Allwmake $targetType $*
 fvAgglomerationMethods/Allwmake $targetType $*
 wmake $targetType waveModels
 
 wmake $targetType engine
-wmake $targetType fvOptions
 wmake $targetType regionCoupled
 
 functionObjects/Allwmake $targetType $*

src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C

@@ -240,7 +240,8 @@ void Foam::displacementComponentLaplacianFvMotionSolver::solve()
     (
         fvm::laplacian
         (
-            diffusivityPtr_->operator()(),
+            dimensionedScalar("viscosity", dimViscosity, 1.0)
+           *diffusivityPtr_->operator()(),
             cellDisplacement_,
             "laplacian(diffusivity,cellDisplacement)"
         )
@@ -249,8 +250,8 @@ void Foam::displacementComponentLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-
     TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    fvOptions.correct(cellDisplacement_);
 }
 
 

src/fvMotionSolver/fvMotionSolvers/componentVelocity/componentLaplacian/velocityComponentLaplacianFvMotionSolver.C

@@ -138,7 +138,8 @@ void Foam::velocityComponentLaplacianFvMotionSolver::solve()
     (
         fvm::laplacian
         (
-            diffusivityPtr_->operator()(),
+            dimensionedScalar("viscosity", dimViscosity, 1.0)
+           *diffusivityPtr_->operator()(),
             cellMotionU_,
             "laplacian(diffusivity,cellMotionU)"
         )
@@ -147,8 +148,8 @@ void Foam::velocityComponentLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-
     TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    fvOptions.correct(cellMotionU_);
 }
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C

@@ -183,7 +183,11 @@ void Foam::displacementSBRStressFvMotionSolver::solve()
     diffusivityPtr_->correct();
     pointDisplacement_.boundaryFieldRef().updateCoeffs();
 
-    surfaceScalarField Df(diffusivityPtr_->operator()());
+    const surfaceScalarField Df
+    (
+        dimensionedScalar("viscosity", dimViscosity, 1.0)
+       *diffusivityPtr_->operator()()
+    );
 
     volTensorField gradCd("gradCd", fvc::grad(cellDisplacement_));
 
@@ -241,8 +245,8 @@ void Foam::displacementSBRStressFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-
     TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    fvOptions.correct(cellDisplacement_);
 }
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C

@@ -342,7 +342,8 @@ void Foam::displacementLaplacianFvMotionSolver::solve()
     (
         fvm::laplacian
         (
-            diffusivity().operator()(),
+            dimensionedScalar("viscosity", dimViscosity, 1.0)
+           *diffusivity().operator()(),
             cellDisplacement_,
             "laplacian(diffusivity,cellDisplacement)"
         )
@@ -351,8 +352,8 @@ void Foam::displacementLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-
     TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    fvOptions.correct(cellDisplacement_);
 }
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/solidBodyDisplacementLaplacian/solidBodyDisplacementLaplacianFvMotionSolver.C

@@ -351,7 +351,8 @@ void Foam::solidBodyDisplacementLaplacianFvMotionSolver::solve()
     (
         fvm::laplacian
         (
-            diffusivity().operator()(),
+            dimensionedScalar("viscosity", dimViscosity, 1.0)
+           *diffusivity().operator()(),
             cellDisplacement_,
             "laplacian(diffusivity,cellDisplacement)"
         )
@@ -360,8 +361,8 @@ void Foam::solidBodyDisplacementLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(TEqn);
-
     TEqn.solveSegregatedOrCoupled(TEqn.solverDict());
+    fvOptions.correct(cellDisplacement_);
 }
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/surfaceAlignedSBRStress/surfaceAlignedSBRStressFvMotionSolver.C

@@ -381,7 +381,12 @@ void Foam::surfaceAlignedSBRStressFvMotionSolver::solve()
         sigmaD_ = magNewSigmaD;
     }
 
-    const surfaceScalarField Df(diffusivity().operator()());
+    const surfaceScalarField Df
+    (
+        dimensionedScalar("viscosity", dimViscosity, 1.0)
+       *diffusivity().operator()()
+    );
+
     pointDisplacement_.boundaryFieldRef().updateCoeffs();
 
     fv::options& fvOptions(fv::options::New(fvMesh_));
@@ -417,6 +422,8 @@ void Foam::surfaceAlignedSBRStressFvMotionSolver::solve()
 
         // Note: solve uncoupled
         DEqn.solveSegregatedOrCoupled(DEqn.solverDict());
+
+        fvOptions.correct(cellDisp);
     }
 }
 

src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C

@@ -131,7 +131,8 @@ void Foam::velocityLaplacianFvMotionSolver::solve()
     (
         fvm::laplacian
         (
-            diffusivityPtr_->operator()(),
+            dimensionedScalar("viscosity", dimViscosity, 1.0)
+           *diffusivityPtr_->operator()(),
             cellMotionU_,
             "laplacian(diffusivity,cellMotionU)"
         )
@@ -140,8 +141,8 @@ void Foam::velocityLaplacianFvMotionSolver::solve()
     );
 
     fvOptions.constrain(UEqn);
-
     UEqn.solveSegregatedOrCoupled(UEqn.solverDict());
+    fvOptions.correct(cellMotionU_);
 }
 
 

src/meshTools/Make/files

@@ -76,6 +76,7 @@ meshSearch/meshSearch.C
 meshSearch/meshSearchFACE_CENTRE_TRISMeshObject.C
 meshSearch/meshSearchMeshObject.C
 
+fields/polyFields.C
 meshTools/meshTools.C
 
 algorithms = algorithms

src/meshTools/fields/polyFields.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "polyFields.H"
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+// Naming to shadow volScalarField::Internal etc.
+// keep synchronized with finiteVolume volFields.C
+
+template<>
+const word DimensionedField<scalar, polyGeoMesh>::typeName
+(
+    "volScalarField::Internal"
+);
+
+template<>
+const word DimensionedField<vector, polyGeoMesh>::typeName
+(
+    "volVectorField::Internal"
+);
+
+template<>
+const word DimensionedField<sphericalTensor, polyGeoMesh>::typeName
+(
+    "volSphericalTensorField::Internal"
+);
+
+template<>
+const word DimensionedField<symmTensor, polyGeoMesh>::typeName
+(
+    "volSymmTensorField::Internal"
+);
+
+template<>
+const word DimensionedField<tensor, polyGeoMesh>::typeName
+(
+    "volTensorField::Internal"
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/meshTools/fields/polyFields.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+InClass
+    Foam::polyFields
+
+Description
+    A polyMesh-based naming and storage for internal volume fields when a
+    Foam::fvMesh is unavailable.  Use sparingly.
+
+SourceFields
+    polyFields.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyFields_H
+#define polyFields_H
+
+#include "DimensionedField.H"
+#include "polyGeoMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef DimensionedField<scalar,polyGeoMesh> polyScalarField;
+typedef DimensionedField<vector,polyGeoMesh> polyVectorField;
+typedef DimensionedField<sphericalTensor,polyGeoMesh> polySphericalTensorField;
+typedef DimensionedField<symmTensor,polyGeoMesh> polySymmTensorField;
+typedef DimensionedField<tensor,polyGeoMesh> polyTensorField;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/fields/polyGeoMesh.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::polyGeoMesh
+
+Description
+    The polyMesh GeoMesh for holding internal fields without an fvMesh.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyGeoMesh_H
+#define polyGeoMesh_H
+
+#include "GeoMesh.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class polyGeoMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyGeoMesh
+:
+    public GeoMesh<polyMesh>
+{
+public:
+
+    // Constructors
+
+        //- Construct from polyMesh reference
+        explicit polyGeoMesh(const polyMesh& mesh)
+        :
+            GeoMesh<polyMesh>(mesh)
+        {}
+
+
+    // Member Functions
+
+        //- Return size
+        static label size(const polyMesh& mesh)
+        {
+            return mesh.nCells();
+        }
+
+        //- Return size
+        label size() const
+        {
+            return size(mesh_);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/regionSplit/regionSplit.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,36 +57,85 @@ void Foam::regionSplit::calcNonCompactRegionSplit
 
     // Take over blockedFaces by seeding a negative number
     // (so is always less than the decomposition)
-    label nUnblocked = 0;
-    forAll(faceData, facei)
+
+    forAll(blockedFace, facei)
     {
-        if (blockedFace.size() && blockedFace[facei])
+        if (blockedFace[facei])
         {
             faceData[facei] = minData(-2);
         }
-        else
+    }
+
+    // Seed all faces on (real) boundaries and faces on cells next to blockFace
+    // (since regions can only occur because of boundaries (or blocked faces))
+
+    PackedBoolList isSeed(mesh().nFaces());
+
+    // Get internal or coupled faces
+    PackedBoolList isConnection(syncTools::getInternalOrCoupledFaces(mesh()));
+
+
+    // 1. Seed (real) boundaries
+    for
+    (
+        label facei = mesh().nInternalFaces();
+        facei < mesh().nFaces();
+        facei++
+    )
+    {
+        if (!isConnection[facei])
         {
-            nUnblocked++;
+            isSeed.set(facei);
         }
     }
 
-    // Seed unblocked faces
-    labelList seedFaces(nUnblocked);
-    List<minData> seedData(nUnblocked);
-    nUnblocked = 0;
-
+    // 2. Seed (internal) faces on cells next to blocked (internal) faces
+    //    (since we can't seed the blocked faces themselves)
 
-    forAll(faceData, facei)
+    if (blockedFace.size())
     {
-        if (blockedFace.empty() || !blockedFace[facei])
+        for (const cell& cFaces : mesh().cells())
         {
-            seedFaces[nUnblocked] = facei;
-            // Seed face with globally unique number
-            seedData[nUnblocked] = minData(globalFaces.toGlobal(facei));
-            nUnblocked++;
+            bool blockedCell = false;
+            label connectedFacei = -1;
+
+            for (const label facei : cFaces)
+            {
+                if (blockedFace[facei])
+                {
+                    blockedCell = true;
+                }
+                else if (isConnection[facei])
+                {
+                    connectedFacei = facei;
+                }
+            }
+
+            if (blockedCell)
+            {
+                isSeed.set(connectedFacei);  // silently ignores -1
+            }
         }
     }
 
+    List<label> seedFaces(isSeed.used());
+    List<minData> seedData(seedFaces.size());
+
+    // Seed face with globally unique number
+    forAll(seedFaces, i)
+    {
+        const label facei = seedFaces[i];
+        seedData[i] = minData(globalFaces.toGlobal(facei));
+    }
+
+    if (debug)
+    {
+        Pout<< "Seeded " << seedFaces.size()
+            << " boundary and internal faces out of"
+            << " total:" << mesh().nInternalFaces()
+            << " boundary:" << mesh().nFaces()-mesh().nInternalFaces()
+            << endl;
+    }
 
     // Propagate information inwards
     FaceCellWave<minData> deltaCalc
@@ -112,7 +161,10 @@ void Foam::regionSplit::calcNonCompactRegionSplit
         }
         else
         {
-            // Unvisited cell -> only possible if surrounded by blocked faces.
+            // Unvisited cell -> only possible if surrounded by blocked faces
+            // since:
+            // - all walls become seed faces
+            // - all faces on cells with blocked faces become seed faces
             // If so make up region from any of the faces
             const cell& cFaces = mesh().cells()[celli];
             label facei = cFaces[0];
@@ -343,7 +395,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
 
 
     // Get the wanted region labels into recvNonLocal
-    labelListList recvNonLocal(Pstream::nProcs());
+    labelListList recvNonLocal;
     Pstream::exchange<labelList, label>(sendNonLocal, recvNonLocal);
 
     // Now we have the wanted compact region labels that proci wants in

src/sampling/Make/files

@@ -16,6 +16,7 @@ sampledSet/sampledSets/sampledSetsGrouping.C
 sampledSet/triSurfaceMeshPointSet/triSurfaceMeshPointSet.C
 sampledSet/uniform/uniformSet.C
 sampledSet/array/arraySet.C
+sampledSet/shortestPath/shortestPathSet.C
 
 surface/cuttingPlane/cuttingPlane.C
 surface/isoSurface/isoSurface.C

src/sampling/sampledSet/sampledSet/sampledSet.C

@@ -29,6 +29,7 @@ License
 #include "meshSearch.H"
 #include "writer.H"
 #include "particle.H"
+#include "SortableList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -402,6 +403,77 @@ void Foam::sampledSet::setSamples
 }
 
 
+Foam::autoPtr<Foam::coordSet> Foam::sampledSet::gather
+(
+    labelList& indexSet
+) const
+{
+    // Combine sampleSet from processors. Sort by curveDist. Return
+    // ordering in indexSet.
+    // Note: only master results are valid
+
+    // Collect data from all processors
+    List<List<point>> gatheredPts(Pstream::nProcs());
+    gatheredPts[Pstream::myProcNo()] = *this;
+    Pstream::gatherList(gatheredPts);
+
+    List<labelList> gatheredSegments(Pstream::nProcs());
+    gatheredSegments[Pstream::myProcNo()] = segments();
+    Pstream::gatherList(gatheredSegments);
+
+    List<scalarList> gatheredDist(Pstream::nProcs());
+    gatheredDist[Pstream::myProcNo()] = curveDist();
+    Pstream::gatherList(gatheredDist);
+
+
+    // Combine processor lists into one big list.
+    List<point> allPts
+    (
+        ListListOps::combine<List<point>>
+        (
+            gatheredPts, accessOp<List<point>>()
+        )
+    );
+    labelList allSegments
+    (
+        ListListOps::combine<labelList>
+        (
+            gatheredSegments, accessOp<labelList>()
+        )
+    );
+    scalarList allCurveDist
+    (
+        ListListOps::combine<scalarList>
+        (
+            gatheredDist, accessOp<scalarList>()
+        )
+    );
+
+
+    if (Pstream::master() && allCurveDist.size() == 0)
+    {
+        WarningInFunction