Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

Allwmake

@@ -26,12 +26,12 @@ else
 fi
 
 # build OpenFOAM libraries and applications
-src/Allwmake
+src/Allwmake $*
-applications/Allwmake
+applications/Allwmake $*
 
 if [ "$1" = doc ]
 then
-    doc/Allwmake
+    doc/Allwmake $*
 fi
 
 # ----------------------------------------------------------------- end-of-file

applications/Allwmake

@@ -16,7 +16,7 @@ wmakeCheckPwd "$WM_PROJECT_DIR/applications" || {
 
 set -x
 
-wmake all utilities
+wmake all utilities $*
-wmake all solvers
+wmake all solvers $*
 
 # ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C

@@ -254,47 +254,85 @@ void Foam::DistributedDelaunayMesh<Triangulation>::findProcessorBoundaryCells
        /Pstream::nProcs()
     );
 
-    std::list<Cell_handle> infinite_cells;
+//    std::list<Cell_handle> infinite_cells;
-    Triangulation::incident_cells
+//    Triangulation::incident_cells
-    (
+//    (
-        Triangulation::infinite_vertex(),
+//        Triangulation::infinite_vertex(),
-        std::back_inserter(infinite_cells)
+//        std::back_inserter(infinite_cells)
-    );
+//    );
-
+//
-    for
+//    for
-    (
+//    (
-        typename std::list<Cell_handle>::iterator vcit = infinite_cells.begin();
+//        typename std::list<Cell_handle>::iterator vcit
-        vcit != infinite_cells.end();
+//            = infinite_cells.begin();
-        ++vcit
+//        vcit != infinite_cells.end();
-    )
+//        ++vcit
-    {
+//    )
-        Cell_handle cit = *vcit;
+//    {
-
+//        Cell_handle cit = *vcit;
-        // Index of infinite vertex in this cell.
+//
-        int i = cit->index(Triangulation::infinite_vertex());
+//        // Index of infinite vertex in this cell.
-
+//        int i = cit->index(Triangulation::infinite_vertex());
-        Cell_handle c = cit->neighbor(i);
+//
-
+//        Cell_handle c = cit->neighbor(i);
-        if (c->unassigned())
+//
-        {
+//        if (c->unassigned())
-            c->cellIndex() = this->getNewCellIndex();
+//        {
-
+//            c->cellIndex() = this->getNewCellIndex();
-            if (checkProcBoundaryCell(c, circumsphereOverlaps))
+//
-            {
+//            if (checkProcBoundaryCell(c, circumsphereOverlaps))
-                cellToCheck.insert(c->cellIndex());
+//            {
-            }
+//                cellToCheck.insert(c->cellIndex());
-        }
+//            }
-    }
+//        }
+//    }
+//
+//
+//    for
+//    (
+//        Finite_cells_iterator cit = Triangulation::finite_cells_begin();
+//        cit != Triangulation::finite_cells_end();
+//        ++cit
+//    )
+//    {
+//        if (cit->parallelDualVertex())
+//        {
+//            if (cit->unassigned())
+//            {
+//                if (checkProcBoundaryCell(cit, circumsphereOverlaps))
+//                {
+//                    cellToCheck.insert(cit->cellIndex());
+//                }
+//            }
+//        }
+//    }
 
 
     for
     (
-        Finite_cells_iterator cit = Triangulation::finite_cells_begin();
+        All_cells_iterator cit = Triangulation::all_cells_begin();
-        cit != Triangulation::finite_cells_end();
+        cit != Triangulation::all_cells_end();
         ++cit
     )
     {
-        if (cit->parallelDualVertex())
+        if (Triangulation::is_infinite(cit))
+        {
+            // Index of infinite vertex in this cell.
+            int i = cit->index(Triangulation::infinite_vertex());
+
+            Cell_handle c = cit->neighbor(i);
+
+            if (c->unassigned())
+            {
+                c->cellIndex() = this->getNewCellIndex();
+
+                if (checkProcBoundaryCell(c, circumsphereOverlaps))
+                {
+                    cellToCheck.insert(c->cellIndex());
+                }
+            }
+        }
+        else if (cit->parallelDualVertex())
         {
             if (cit->unassigned())
             {

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C

@@ -1046,6 +1046,18 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         )
     ),
     decomposition_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::conformalVoronoiMesh::~conformalVoronoiMesh()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::conformalVoronoiMesh::initialiseForMotion()
 {
     if (foamyHexMeshControls().objOutput())
     {
@@ -1061,7 +1073,10 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
                 runTime_,
                 rndGen_,
                 geometryToConformTo_,
-                foamyHexMeshDict.subDict("backgroundMeshDecomposition")
+                foamyHexMeshControls().foamyHexMeshDict().subDict
+                (
+                    "backgroundMeshDecomposition"
+                )
             )
         );
     }
@@ -1125,13 +1140,53 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
 }
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+void Foam::conformalVoronoiMesh::initialiseForConformation()
+{
+    if (Pstream::parRun())
+    {
+        decomposition_.reset
+        (
+            new backgroundMeshDecomposition
+            (
+                runTime_,
+                rndGen_,
+                geometryToConformTo_,
+                foamyHexMeshControls().foamyHexMeshDict().subDict
+                (
+                    "backgroundMeshDecomposition"
+                )
+            )
+        );
+    }
 
-Foam::conformalVoronoiMesh::~conformalVoronoiMesh()
+    insertInitialPoints();
-{}
 
+    insertFeaturePoints();
+
+    // Improve the guess that the backgroundMeshDecomposition makes with the
+    // initial positions.  Use before building the surface conformation to
+    // better balance the surface conformation load.
+    distributeBackground(*this);
+
+    buildSurfaceConformation();
+
+    // The introduction of the surface conformation may have distorted the
+    // balance of vertices, distribute if necessary.
+    distributeBackground(*this);
+
+    if (Pstream::parRun())
+    {
+        sync(decomposition_().procBounds());
+    }
+
+    cellSizeMeshOverlapsBackground();
+
+    if (foamyHexMeshControls().printVertexInfo())
+    {
+        printVertexInfo(Info);
+    }
+}
 
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void Foam::conformalVoronoiMesh::move()
 {

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -1038,6 +1038,10 @@ public:
 
     // Member Functions
 
+        void initialiseForMotion();
+
+        void initialiseForConformation();
+
         //- Move the vertices according to the controller, re-conforming to the
         //  surface as required
         void move();

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -1787,6 +1787,11 @@ void Foam::conformalVoronoiMesh::indexDualVertices
         }
     }
 
+    OBJstream snapping1("snapToSurface1.obj");
+    OBJstream snapping2("snapToSurface2.obj");
+    OFstream tetToSnapTo("tetsToSnapTo.obj");
+    label offset = 0;
+
     for
     (
         Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -1892,6 +1897,88 @@ void Foam::conformalVoronoiMesh::indexDualVertices
                 }
             }
 
+            {
+                // Snapping points far outside
+                if (cit->boundaryDualVertex() && !cit->parallelDualVertex())
+                {
+                    pointFromPoint dual = cit->dual();
+
+                    pointIndexHit hitInfo;
+                    label surfHit;
+
+                    // Find nearest surface point
+                    geometryToConformTo_.findSurfaceNearest
+                    (
+                        dual,
+                        sqr(targetCellSize(dual)),
+                        hitInfo,
+                        surfHit
+                    );
+
+                    if (!hitInfo.hit())
+                    {
+                        // Project dual to nearest point on tet
+
+                        tetPointRef tet
+                        (
+                            topoint(cit->vertex(0)->point()),
+                            topoint(cit->vertex(1)->point()),
+                            topoint(cit->vertex(2)->point()),
+                            topoint(cit->vertex(3)->point())
+                        );
+
+                        pointFromPoint nearestPointOnTet =
+                            tet.nearestPoint(dual).rawPoint();
+
+                        // Get nearest point on surface from tet.
+                        geometryToConformTo_.findSurfaceNearest
+                        (
+                            nearestPointOnTet,
+                            sqr(targetCellSize(nearestPointOnTet)),
+                            hitInfo,
+                            surfHit
+                        );
+
+                        vector snapDir = nearestPointOnTet - dual;
+                        snapDir /= mag(snapDir) + SMALL;
+
+                        drawDelaunayCell(tetToSnapTo, cit, offset);
+                        offset += 1;
+
+                        vectorField norm(1);
+                        allGeometry_[surfHit].getNormal
+                        (
+                            List<pointIndexHit>(1, hitInfo),
+                            norm
+                        );
+                        norm[0] /= mag(norm[0]) + SMALL;
+
+                        if
+                        (
+                            hitInfo.hit()
+                         && (mag(snapDir & norm[0]) > 0.5)
+                        )
+                        {
+                            snapping1.write
+                            (
+                                linePointRef(dual, nearestPointOnTet)
+                            );
+
+                            snapping2.write
+                            (
+                                linePointRef
+                                (
+                                    nearestPointOnTet,
+                                    hitInfo.hitPoint()
+                                )
+                            );
+
+                            pts[cit->cellIndex()] = hitInfo.hitPoint();
+                        }
+                    }
+                }
+            }
+
             if (cit->boundaryDualVertex())
             {
                 if (cit->featureEdgeDualVertex())

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C

@@ -58,7 +58,7 @@ List<Vb::Point> pointFile::initialPoints() const
         IOobject
         (
             pointFileName_.name(),
-            foamyHexMesh_.time().constant(),
+            foamyHexMesh_.time().timeName(),
             foamyHexMesh_.time(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE

applications/utilities/mesh/generation/foamyHexMesh/foamyHexMesh.C

@@ -45,12 +45,22 @@ int main(int argc, char *argv[])
         "check all surface geometry for quality"
     );
 
+    Foam::argList::addBoolOption
+    (
+        "conformationOnly",
+        "conform to the initial points without any point motion"
+    );
+
+    #include "addOverwriteOption.H"
+
     #include "setRootCase.H"
     #include "createTime.H"
 
     runTime.functionObjects().off();
 
     const bool checkGeometry = args.optionFound("checkGeometry");
+    const bool conformationOnly = args.optionFound("conformationOnly");
+    const bool overwrite = args.optionFound("overwrite");
 
     IOdictionary foamyHexMeshDict
     (
@@ -104,16 +114,33 @@ int main(int argc, char *argv[])
     conformalVoronoiMesh mesh(runTime, foamyHexMeshDict);
 
 
-    while (runTime.loop())
+    if (conformationOnly)
     {
-        Info<< nl << "Time = " << runTime.timeName() << endl;
+        mesh.initialiseForConformation();
 
-        mesh.move();
+        if (!overwrite)
+        {
+            runTime++;
+        }
 
-        Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        mesh.writeMesh(runTime.timeName());
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << endl;
     }
+    else
+    {
+        mesh.initialiseForMotion();
+
+        while (runTime.loop())
+        {
+            Info<< nl << "Time = " << runTime.timeName() << endl;
+
+            mesh.move();
+
+            Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+                << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+                << endl;
+        }
+    }
+
 
     Info<< nl << "End" << nl << endl;
 

src/Allwmake

@@ -24,7 +24,7 @@ wmakeLnInclude OpenFOAM
 wmakeLnInclude OSspecific/${WM_OSTYPE:-POSIX}
 Pstream/Allwmake $*
 
-OSspecific/${WM_OSTYPE:-POSIX}/Allwmake
+OSspecific/${WM_OSTYPE:-POSIX}/Allwmake $*
 wmake $makeType OpenFOAM
 
 wmake $makeType fileFormats

src/OSspecific/POSIX/Allwmake

@@ -1,5 +1,6 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # run from this directory
+makeType=${1:-libo}
 
 unset COMP_FLAGS LINK_FLAGS
 
@@ -19,6 +20,6 @@ fi
 
 
 # make (non-shared) object
-wmake libo
+wmake $makeType 
 
 # ----------------------------------------------------------------- end-of-file

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomerate.C

@@ -26,177 +26,7 @@ License
 #include "pairGAMGAgglomeration.H"
 #include "lduAddressing.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
-
-Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
-(
-    label& nCoarseCells,
-    const lduAddressing& fineMatrixAddressing,
-    const scalarField& faceWeights
-)
-{
-    const label nFineCells = fineMatrixAddressing.size();
-
-    const labelUList& upperAddr = fineMatrixAddressing.upperAddr();
-    const labelUList& lowerAddr = fineMatrixAddressing.lowerAddr();
-
-    // For each cell calculate faces
-    labelList cellFaces(upperAddr.size() + lowerAddr.size());
-    labelList cellFaceOffsets(nFineCells + 1);
-
-    // memory management
-    {
-        labelList nNbrs(nFineCells, 0);
-
-        forAll(upperAddr, facei)
-        {
-            nNbrs[upperAddr[facei]]++;
-        }
-
-        forAll(lowerAddr, facei)
-        {
-            nNbrs[lowerAddr[facei]]++;
-        }
-
-        cellFaceOffsets[0] = 0;
-        forAll(nNbrs, celli)
-        {
-            cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
-        }
-
-        // reset the whole list to use as counter
-        nNbrs = 0;
-
-        forAll(upperAddr, facei)
-        {
-            cellFaces
-            [
-                cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
-            ] = facei;
-
-            nNbrs[upperAddr[facei]]++;
-        }
-
-        forAll(lowerAddr, facei)
-        {
-            cellFaces
-            [
-                cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
-            ] = facei;
-
-            nNbrs[lowerAddr[facei]]++;
-        }
-    }
-
-
-    // go through the faces and create clusters
-
-    tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
-    labelField& coarseCellMap = tcoarseCellMap();
-
-    nCoarseCells = 0;
-
-    for (label celli=0; celli<nFineCells; celli++)
-    {
-        if (coarseCellMap[celli] < 0)
-        {
-            label matchFaceNo = -1;
-            scalar maxFaceWeight = -GREAT;
-
-            // check faces to find ungrouped neighbour with largest face weight
-            for
-            (
-                label faceOs=cellFaceOffsets[celli];
-                faceOs<cellFaceOffsets[celli+1];
-                faceOs++
-            )
-            {
-                label facei = cellFaces[faceOs];
-
-                // I don't know whether the current cell is owner or neighbour.
-                // Therefore I'll check both sides
-                if
-                (
-                    coarseCellMap[upperAddr[facei]] < 0
-                 && coarseCellMap[lowerAddr[facei]] < 0
-                 && faceWeights[facei] > maxFaceWeight
-                )
-                {
-                    // Match found. Pick up all the necessary data
-                    matchFaceNo = facei;
-                    maxFaceWeight = faceWeights[facei];
-                }
-            }
-
-            if (matchFaceNo >= 0)
-            {
-                // Make a new group
-                coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
-                coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
-                nCoarseCells++;
-            }
-            else
-            {
-                // No match. Find the best neighbouring cluster and
-                // put the cell there
-                label clusterMatchFaceNo = -1;
-                scalar clusterMaxFaceCoeff = -GREAT;
-
-                for
-                (
-                    label faceOs=cellFaceOffsets[celli];
-                    faceOs<cellFaceOffsets[celli+1];
-                    faceOs++
-                )
-                {
-                    label facei = cellFaces[faceOs];
-
-                    if (faceWeights[facei] > clusterMaxFaceCoeff)
-                    {
-                        clusterMatchFaceNo = facei;
-                        clusterMaxFaceCoeff = faceWeights[facei];
-                    }
-                }
-
-                if (clusterMatchFaceNo >= 0)
-                {
-                    // Add the cell to the best cluster
-                    coarseCellMap[celli] = max
-                    (
-                        coarseCellMap[upperAddr[clusterMatchFaceNo]],
-                        coarseCellMap[lowerAddr[clusterMatchFaceNo]]
-                    );
-                }
-            }
-        }
-    }
-
-
-    // Check that all cells are part of clusters,
-    // if not create single-cell "clusters" for each
-    for (label celli=0; celli<nFineCells; celli++)
-    {
-        if (coarseCellMap[celli] < 0)
-        {
-            coarseCellMap[celli] = nCoarseCells;
-            nCoarseCells++;
-        }
-    }
-
-    // Reverse the map ordering to improve the next level of agglomeration
-    // (doesn't always help and is sometimes detrimental)
-    nCoarseCells--;
-
-    forAll(coarseCellMap, celli)
-    {
-        coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
-    }
-
-    nCoarseCells++;
-
-    return tcoarseCellMap;
-}
-
 
 void Foam::pairGAMGAgglomeration::agglomerate
 (
@@ -285,4 +115,187 @@ void Foam::pairGAMGAgglomeration::agglomerate
 }
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::labelField> Foam::pairGAMGAgglomeration::agglomerate
+(
+    label& nCoarseCells,
+    const lduAddressing& fineMatrixAddressing,
+    const scalarField& faceWeights
+)
+{
+    const label nFineCells = fineMatrixAddressing.size();
+
+    const labelUList& upperAddr = fineMatrixAddressing.upperAddr();
+    const labelUList& lowerAddr = fineMatrixAddressing.lowerAddr();
+
+    // For each cell calculate faces
+    labelList cellFaces(upperAddr.size() + lowerAddr.size());
+    labelList cellFaceOffsets(nFineCells + 1);
+
+    // memory management
+    {
+        labelList nNbrs(nFineCells, 0);
+
+        forAll(upperAddr, facei)
+        {
+            nNbrs[upperAddr[facei]]++;
+        }
+
+        forAll(lowerAddr, facei)
+        {
+            nNbrs[lowerAddr[facei]]++;
+        }
+
+        cellFaceOffsets[0] = 0;
+        forAll(nNbrs, celli)
+        {
+            cellFaceOffsets[celli+1] = cellFaceOffsets[celli] + nNbrs[celli];
+        }
+
+        // reset the whole list to use as counter
+        nNbrs = 0;
+
+        forAll(upperAddr, facei)
+        {
+            cellFaces
+            [
+                cellFaceOffsets[upperAddr[facei]] + nNbrs[upperAddr[facei]]
+            ] = facei;
+
+            nNbrs[upperAddr[facei]]++;
+        }
+
+        forAll(lowerAddr, facei)
+        {
+            cellFaces
+            [
+                cellFaceOffsets[lowerAddr[facei]] + nNbrs[lowerAddr[facei]]
+            ] = facei;
+
+            nNbrs[lowerAddr[facei]]++;
+        }
+    }
+
+
+    // go through the faces and create clusters
+
+    tmp<labelField> tcoarseCellMap(new labelField(nFineCells, -1));
+    labelField& coarseCellMap = tcoarseCellMap();
+
+    nCoarseCells = 0;
+    label celli;
+
+    for (label cellfi=0; cellfi<nFineCells; cellfi++)
+    {
+        // Change cell ordering depending on direction for this level
+        celli = forward_ ? cellfi : nFineCells - cellfi - 1;
+
+        if (coarseCellMap[celli] < 0)
+        {
+            label matchFaceNo = -1;
+            scalar maxFaceWeight = -GREAT;
+
+            // check faces to find ungrouped neighbour with largest face weight
+            for
+            (
+                label faceOs=cellFaceOffsets[celli];
+                faceOs<cellFaceOffsets[celli+1];
+                faceOs++
+            )
+            {
+                label facei = cellFaces[faceOs];
+
+                // I don't know whether the current cell is owner or neighbour.
+                // Therefore I'll check both sides
+                if
+                (
+                    coarseCellMap[upperAddr[facei]] < 0
+                 && coarseCellMap[lowerAddr[facei]] < 0
+                 && faceWeights[facei] > maxFaceWeight
+                )
+                {
+                    // Match found. Pick up all the necessary data
+                    matchFaceNo = facei;
+                    maxFaceWeight = faceWeights[facei];
+                }
+            }
+
+            if (matchFaceNo >= 0)
+            {
+                // Make a new group
+                coarseCellMap[upperAddr[matchFaceNo]] = nCoarseCells;
+                coarseCellMap[lowerAddr[matchFaceNo]] = nCoarseCells;
+                nCoarseCells++;
+            }
+            else
+            {
+                // No match. Find the best neighbouring cluster and
+                // put the cell there
+                label clusterMatchFaceNo = -1;
+                scalar clusterMaxFaceCoeff = -GREAT;
+
+                for
+                (
+                    label faceOs=cellFaceOffsets[celli];
+                    faceOs<cellFaceOffsets[celli+1];
+                    faceOs++
+                )
+                {
+                    label facei = cellFaces[faceOs];
+
+                    if (faceWeights[facei] > clusterMaxFaceCoeff)
+                    {
+                        clusterMatchFaceNo = facei;
+                        clusterMaxFaceCoeff = faceWeights[facei];
+                    }
+                }
+
+                if (clusterMatchFaceNo >= 0)
+                {
+                    // Add the cell to the best cluster
+                    coarseCellMap[celli] = max
+                    (
+                        coarseCellMap[upperAddr[clusterMatchFaceNo]],
+                        coarseCellMap[lowerAddr[clusterMatchFaceNo]]
+                    );
+                }
+            }
+        }
+    }
+
+    // Check that all cells are part of clusters,
+    // if not create single-cell "clusters" for each
+    for (label cellfi=0; cellfi<nFineCells; cellfi++)
+    {
+        // Change cell ordering depending on direction for this level
+        celli = forward_ ? cellfi : nFineCells - cellfi - 1;
+
+        if (coarseCellMap[celli] < 0)
+        {
+            coarseCellMap[celli] = nCoarseCells;
+            nCoarseCells++;
+        }
+    }
+
+    if (!forward_)
+    {
+        nCoarseCells--;
+
+        forAll(coarseCellMap, celli)
+        {
+            coarseCellMap[celli] = nCoarseCells - coarseCellMap[celli];
+        }
+
+        nCoarseCells++;
+    }
+
+    // Reverse the map ordering for the next level
+    // to improve the next level of agglomeration
+    forward_ = !forward_;
+
+    return tcoarseCellMap;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,7 +29,8 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(pairGAMGAgglomeration, 0);
+    defineTypeNameAndDebug(pairGAMGAgglomeration, 0);
+    bool pairGAMGAgglomeration::forward_(true);
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/pairGAMGAgglomeration/pairGAMGAgglomeration.H

@@ -56,6 +56,9 @@ class pairGAMGAgglomeration
         //- Number of levels to merge, 1 = don't merge, 2 = merge pairs etc.
         label mergeLevels_;
 
+        //- Direction of cell loop for the current level
+        static bool forward_;
+
 
 protected:
 
@@ -97,7 +100,6 @@ public:
             const lduAddressing& fineMatrixAddressing,
             const scalarField& faceWeights
         );
-
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H

@@ -230,7 +230,8 @@ class GAMGSolver
             const lduMatrix& m,
             const FieldField<Field, scalar>& interfaceBouCoeffs,
             const lduInterfaceFieldPtrsList& interfaces,
-            const scalarField& source,
+            const labelList& restrictAddressing,
+            const scalarField& psiC,
             const direction cmpt
         ) const;
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverInterpolate.C

@@ -34,7 +34,8 @@ void Foam::GAMGSolver::interpolate
     const lduMatrix& m,
     const FieldField<Field, scalar>& interfaceBouCoeffs,
     const lduInterfaceFieldPtrsList& interfaces,
-    const scalarField& source,
+    const labelList& restrictAddressing,
+    const scalarField& psiC,
     const direction cmpt
 ) const
 {
@@ -80,6 +81,30 @@ void Foam::GAMGSolver::interpolate
     {
         psiPtr[celli] = -ApsiPtr[celli]/(diagPtr[celli]);
     }
+
+    register const label nCCells = psiC.size();
+    scalarField corrC(nCCells, 0);
+    scalarField diagC(nCCells, 0);
+
+    for (register label celli=0; celli<nCells; celli++)
+    {
+        corrC[restrictAddressing[celli]] += diagPtr[celli]*psiPtr[celli];
+        diagC[restrictAddressing[celli]] += diagPtr[celli];
+        //corrC[restrictAddressing[celli]] += psiPtr[celli]/diagPtr[celli];
+        //diagC[restrictAddressing[celli]] += 1.0/diagPtr[celli];
+        //corrC[restrictAddressing[celli]] += psiPtr[celli];
+        //diagC[restrictAddressing[celli]] += 1.0;
+    }
+
+    for (register label ccelli=0; ccelli<nCCells; ccelli++)
+    {
+        corrC[ccelli] = psiC[ccelli] - corrC[ccelli]/diagC[ccelli];
+    }
+
+    for (register label celli=0; celli<nCells; celli++)
+    {
+        psiPtr[celli] += corrC[restrictAddressing[celli]];
+    }
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C

@@ -313,7 +313,7 @@ void Foam::GAMGSolver::Vcycle
             scalarField& ACfRef =
                 const_cast<scalarField&>(ACf.operator const scalarField&());
 
-            if (interpolateCorrection_)
+            if (interpolateCorrection_ && leveli < coarsestLevel - 2)
             {
                 interpolate
                 (
@@ -322,7 +322,8 @@ void Foam::GAMGSolver::Vcycle
                     matrixLevels_[leveli],
                     interfaceLevelsBouCoeffs_[leveli],
                     interfaceLevels_[leveli],
-                    coarseSources[leveli],
+                    agglomeration_.restrictAddressing(leveli + 1),
+                    coarseCorrFields[leveli + 1],
                     cmpt
                 );
             }
@@ -382,7 +383,8 @@ void Foam::GAMGSolver::Vcycle
             matrix_,
             interfaceBouCoeffs_,
             interfaces_,
-            finestResidual,
+            agglomeration_.restrictAddressing(0),
+            coarseCorrFields[0],
             cmpt
         );
     }

src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomerate.C

@@ -30,11 +30,6 @@ Description
 #include "fvMesh.H"
 #include "syncTools.H"
 
-//extern "C"
-//{
-//#   include "mgridgen.h"
-//}
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::MGridGenGAMGAgglomeration::

src/fvAgglomerationMethods/MGridGenGamgAgglomeration/MGridGenGAMGAgglomeration.C

@@ -59,7 +59,11 @@ Foam::MGridGenGAMGAgglomeration::MGridGenGAMGAgglomeration
 
     // Start geometric agglomeration from the cell volumes and areas of the mesh
     scalarField* VPtr = const_cast<scalarField*>(&fvMesh_.cellVolumes());
-    SubField<vector> Sf(fvMesh_.faceAreas(), fvMesh_.nInternalFaces());
+
+    vectorField magFaceAreas(vector::one*mag(fvMesh_.faceAreas()));
+    SubField<vector> Sf(magFaceAreas, fvMesh_.nInternalFaces());
+    //SubField<vector> Sf(fvMesh_.faceAreas(), fvMesh_.nInternalFaces());
+
     vectorField* SfPtr = const_cast<vectorField*>
     (
         &Sf.operator const vectorField&()

src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationMotionSolver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,9 +36,6 @@ Description
 
     Tables are in the \a constant/tables directory.
 
-Note
-    could be a motionSolver - does not use any fvMesh structure.
-
 SourceFiles
     displacementInterpolationMotionSolver.C
 

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -253,6 +253,18 @@ Foam::displacementLayeredMotionMotionSolver::faceZoneEvaluate
         // Only on boundary faces - follow boundary conditions
         fld = vectorField(pointDisplacement_, meshPoints);
     }
+    else if (type == "uniformFollow")
+    {
+        // Reads name of name of patch. Then get average point dislacement on
+        // patch. That becomes the value of fld.
+        const word patchName(dict.lookup("patch"));
+        label patchID = mesh().boundaryMesh().findPatchID(patchName);
+        pointField pdf
+        (
+            pointDisplacement_.boundaryField()[patchID].patchInternalField()
+        );
+        fld = gAverage(pdf);
+    }
     else
     {
         FatalIOErrorIn
@@ -399,22 +411,6 @@ Info<< "For cellZone:" << cellZoneI
         // Implement real bc.
         patchDisp[patchI].correctBoundaryConditions();
 
-
-//Info<< "Writing displacement for faceZone " << fz.name()
-//    << " to " << patchDisp[patchI].name() << endl;
-//patchDisp[patchI].write();
-
-//        // Copy into pointDisplacement for other fields to use
-//        forAll(isZonePoint, pointI)
-//        {
-//            if (isZonePoint[pointI])
-//            {
-//                pointDisplacement_[pointI] = patchDisp[patchI][pointI];
-//            }
-//        }
-//        pointDisplacement_.correctBoundaryConditions();
-
-
         patchI++;
     }
 
@@ -423,37 +419,40 @@ Info<< "For cellZone:" << cellZoneI
     // ~~~~~
     // solving the interior is just interpolating
 
-//    // Get normalised distance
+    if (debug)
-//    pointScalarField distance
+    {
-//    (
+        // Get normalised distance
-//        IOobject
+        pointScalarField distance
-//        (
+        (
-//            "distance",
+            IOobject
-//            mesh().time().timeName(),
+            (
-//            mesh(),
+                "distance",
-//            IOobject::NO_READ,
+                mesh().time().timeName(),
-//            IOobject::NO_WRITE,
+                mesh(),
-//            false
+                IOobject::NO_READ,
-//        ),
+                IOobject::NO_WRITE,
-//        pointMesh::New(mesh()),
+                false
-//        dimensionedScalar("distance", dimLength, 0.0)
+            ),
-//    );
+            pointMesh::New(mesh()),
-//    forAll(distance, pointI)
+            dimensionedScalar("distance", dimLength, 0.0)
-//    {
+        );
-//        if (isZonePoint[pointI])
+        forAll(distance, pointI)
-//        {
+        {
-//            scalar d1 = patchDist[0][pointI];
+            if (isZonePoint[pointI])
-//            scalar d2 = patchDist[1][pointI];
+            {
-//            if (d1+d2 > SMALL)
+                scalar d1 = patchDist[0][pointI];
-//            {
+                scalar d2 = patchDist[1][pointI];
-//                scalar s = d1/(d1+d2);
+                if (d1+d2 > SMALL)
-//                distance[pointI] = s;
+                {
-//            }
+                    scalar s = d1/(d1+d2);
-//        }
+                    distance[pointI] = s;
-//    }
+                }
-//    Info<< "Writing distance pointScalarField to " << mesh().time().timeName()
+            }
-//        << endl;
+        }
-//    distance.write();
+        Info<< "Writing distance pointScalarField to "
+            << mesh().time().timeName() << endl;
+        distance.write();
+    }
 
     // Average
     forAll(pointDisplacement_, pointI)
@@ -470,7 +469,6 @@ Info<< "For cellZone:" << cellZoneI
               + s*patchDisp[1][pointI];
         }
     }
-    pointDisplacement_.correctBoundaryConditions();
 }
 
 
@@ -484,9 +482,7 @@ displacementLayeredMotionMotionSolver
 )
 :
     displacementMotionSolver(mesh, dict, typeName)
-{
+{}
-    pointDisplacement_.correctBoundaryConditions();
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -518,6 +514,9 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
     // the motionSolver accordingly
     movePoints(mesh().points());
 
+    // Apply boundary conditions
+    pointDisplacement_.boundaryField().updateCoeffs();
+
     // Apply all regions (=cellZones)
 
     const dictionary& regionDicts = coeffDict().subDict("regions");
@@ -544,6 +543,9 @@ void Foam::displacementLayeredMotionMotionSolver::solve()
 
         cellZoneSolve(zoneI, regionDict);
     }
+
+    // Update pointDisplacement for solved values
+    pointDisplacement_.correctBoundaryConditions();
 }
 
 

src/fvMotionSolver/fvMotionSolvers/displacement/layeredSolver/displacementLayeredMotionMotionSolver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,11 +28,6 @@ Description
     Mesh motion solver for an (multi-block) extruded fvMesh. Gets given the
     structure of the mesh blocks and boundary conditions on these blocks.
 
-    Note: should not be an fvMotionSolver but just a motionSolver. Only here
-    so we can reuse displacementFvMotionSolver functionality (e.g. surface
-    following boundary conditions)
-
-
     The displacementLayeredMotionCoeffs subdict of dynamicMeshDict specifies
     per region (=cellZone) the boundary conditions on two opposing patches
     (=faceZones). It then interpolates the boundary values using topological
@@ -44,6 +39,9 @@ Description
             Use this for faceZones on boundary faces (so it uses the
             proper boundary conditions on the pointDisplacement).
 
+    uniformFollow: like 'follow' but takes the average value of
+            a specified 'patch' (which is not necessarily colocated)
+
     fixedValue: fixed value.
 
     timeVaryingUniformFixedValue: table-driven fixed value.

src/postProcessing/Allwmake

@@ -3,7 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 makeType=${1:-libso}
 set -x
 
-wmake libo postCalc
+wmake ${1:-libo} postCalc
 wmake $makeType foamCalcFunctions
 
 functionObjects/Allwmake $*

src/sampling/probes/patchProbesTemplates.C

@@ -43,7 +43,9 @@ void Foam::patchProbes::sampleAndWrite
         unsigned int w = IOstream::defaultPrecision() + 7;
         OFstream& probeStream = *probeFilePtrs_[vField.name()];
 
-        probeStream << setw(w) << vField.time().value();
+        probeStream
+            << setw(w)
+            << vField.time().timeToUserTime(vField.time().value());
 
         forAll(values, probeI)
         {
@@ -67,7 +69,9 @@ void Foam::patchProbes::sampleAndWrite
         unsigned int w = IOstream::defaultPrecision() + 7;
         OFstream& probeStream = *probeFilePtrs_[sField.name()];
 
-        probeStream << setw(w) << sField.time().value();
+        probeStream
+            << setw(w)
+            << sField.time().timeToUserTime(sField.time().value());
 
         forAll(values, probeI)
         {

src/sampling/probes/probesTemplates.C

@@ -76,7 +76,7 @@ void Foam::probes::sampleAndWrite
         unsigned int w = IOstream::defaultPrecision() + 7;
         OFstream& os = *probeFilePtrs_[vField.name()];
 
-        os  << setw(w) << vField.time().value();
+        os  << setw(w) << vField.time().timeToUserTime(vField.time().value());
 
         forAll(values, probeI)
         {
@@ -90,17 +90,17 @@ void Foam::probes::sampleAndWrite
 template<class Type>
 void Foam::probes::sampleAndWrite
 (
-    const GeometricField<Type, fvsPatchField, surfaceMesh>& vField
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
 )
 {
-    Field<Type> values(sample(vField));
+    Field<Type> values(sample(sField));
 
     if (Pstream::master())
     {
         unsigned int w = IOstream::defaultPrecision() + 7;
-        OFstream& os = *probeFilePtrs_[vField.name()];
+        OFstream& os = *probeFilePtrs_[sField.name()];
 
-        os  << setw(w) << vField.time().value();
+        os  << setw(w) << sField.time().timeToUserTime(sField.time().value());
 
         forAll(values, probeI)
         {
@@ -259,7 +259,7 @@ template<class Type>
 Foam::tmp<Foam::Field<Type> >
 Foam::probes::sample
 (
-    const GeometricField<Type, fvsPatchField, surfaceMesh>& vField
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
 ) const
 {
     const Type unsetVal(-VGREAT*pTraits<Type>::one);
@@ -275,7 +275,7 @@ Foam::probes::sample
     {
         if (faceList_[probeI] >= 0)
         {
-            values[probeI] = vField[faceList_[probeI]];
+            values[probeI] = sField[faceList_[probeI]];
         }
     }
 

src/triSurface/triSurface/triSurface.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -494,6 +494,7 @@ void Foam::triSurface::write
         (
             "triSurface::write(const fileName&, const word&, const bool)"
         )   << "unknown file extension " << ext
+            << " for file " << name
             << ". Supported extensions are '.ftr', '.stl', '.stlb', "
             << "'.gts', '.obj', '.vtk'"
             << ", '.off', '.dx', '.smesh', '.ac' and '.tri'"

tutorials/multiphase/twoPhaseEulerFoam/bed/system/fvSolution

@@ -53,6 +53,12 @@ solvers
         relTol          0;
     }
 
+    alpha1
+    {
+        nAlphaCorr      2;
+        nAlphaSubCycles 3;
+    }
+
     alpha
     {
         solver          PCG;
@@ -72,12 +78,9 @@ PIMPLE
 {
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
-    nAlphaCorr      2;
-    nAlphaSubCycles 3;
     correctAlpha    yes;
     pRefCell        0;
     pRefValue       0;
 }
 
-
 // ************************************************************************* //

wmake/wmake

@@ -57,6 +57,7 @@ or a special target:
   libo              build statically linked lib (.o)
   libso             build dynamically linked lib (.so)
   jar               build Java jar
+  dep               build lnInclude and dependencies only
 
 USAGE
     exit 1
@@ -290,7 +291,7 @@ OBJECTS_DIR=$MakeDir/$WM_OPTIONS
 touch $OBJECTS_DIR/dontIncludeDeps
 
 case "$makeType" in
-lib | libo | libso )
+lib | libo | libso | dep )
     $make -s -f $WM_DIR/Makefile MAKE_DIR=$MakeDir INCLUDE_DEPS=$OBJECTS_DIR/dontIncludeDeps lnInclude/uptodate
     ;;
 esac
@@ -305,8 +306,11 @@ rc=$?
 # Make the object files and link
 #------------------------------------------------------------------------------
 
-cmd="$make -f $WM_DIR/Makefile MAKE_DIR=$MakeDir INCLUDE_DEPS=$OBJECTS_DIR/includeDeps $makeType"
+if [ "$makeType" != dep ]
-# echo "cmd=$cmd"
+then
-exec $cmd
+    cmd="$make -f $WM_DIR/Makefile MAKE_DIR=$MakeDir INCLUDE_DEPS=$OBJECTS_DIR/includeDeps $makeType"
+    # echo "cmd=$cmd"
+    exec $cmd
+fi
 
 #------------------------------------------------------------------------------