New fvMotionSolver. Small fix to foamDebugSwitches

applications/utilities/miscellaneous/foamDebugSwitches/Make/options

@@ -1,26 +1,43 @@
 EXE_LIBS = \
-    -lLESfilters \
+    -lautoMesh \
-    -lODE \
-    -lcompressibleLESmodels \
-    \
     -lbasicThermophysicalModels \
-    -lfiniteVolume \
+    -lchemistryModel \
     -lcombustionThermophysicalModels \
+    -lcompressibleLESmodels \
     -lcompressibleTurbulenceModels \
+    -ldecompositionMethods \
     -ldieselSpray \
+    -ldynamicFvMesh \
+    -ldynamicMesh \
     -ledgeMesh \
     -lengine \
     -lerrorEstimation \
-    -ldynamicMesh \
+    -lfiniteVolume \
+    -lforces \
+    -lfvMotionSolvers \
+    -lincompressibleLESmodels \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModels \
-    -lincompressibleLESmodels \
+    -linterfaceProperties \
+    -llagrangianIntermediate \
     -llagrangian \
     -llaminarFlameSpeedModels \
+    -lLESdeltas \
+    -lLESfilters \
+    -lliquidMixture \
     -lliquids \
+    -lmeshTools \
+    -lODE \
+    -lOpenFOAM \
+    -lpdf \
+    -lphaseModel \
+    -lradiation \
     -lrandomProcesses \
     -lsampling \
+    -lsolidMixture \
+    -lsolids \
     -lspecie \
     -lthermophysicalFunctions \
-    -ltriSurface \
+    -ltopoChangerFvMesh \
-    -lpdf
+    -ltriSurface
+

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C

@@ -2882,8 +2882,7 @@ void Foam::autoHexMeshDriver::addLayers
 
         if (debug_)
         {
-            Info<< "*** Writing layer mesh to "
+            Info<< "Writing layer mesh to " << mesh_.time().timeName() << endl;
-                << mesh_.time().timeName() << endl;
             newMesh.write();
             cellSet addedCellSet
             (

src/fvMotionSolver/Make/files

@@ -1,5 +1,6 @@
 fvMotionSolvers/fvMotionSolver/fvMotionSolver.C
 fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C
+fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C
 fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C
 fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C
 

src/fvMotionSolver/fvMotionSolvers/displacement/interpolation/displacementInterpolationFvMotionSolver.C

@@ -0,0 +1,465 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "displacementInterpolationFvMotionSolver.H"
+#include "addToRunTimeSelectionTable.H"
+#include "SortableList.H"
+#include "IOList.H"
+#include "Tuple2.H"
+#include "mapPolyMesh.H"
+#include "interpolateXY.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(displacementInterpolationFvMotionSolver, 0);
+
+    addToRunTimeSelectionTable
+    (
+        fvMotionSolver,
+        displacementInterpolationFvMotionSolver,
+        dictionary
+    );
+
+    template<>
+    const word IOList<Tuple2<scalar, vector> >::typeName("scalarVectorTable");
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::displacementInterpolationFvMotionSolver::
+displacementInterpolationFvMotionSolver
+(
+    const polyMesh& mesh,
+    Istream& msData
+)
+:
+    fvMotionSolver(mesh),
+    points0_
+    (
+        pointIOField
+        (
+            IOobject
+            (
+                "points",
+                mesh.time().constant(),
+                polyMesh::meshSubDir,
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        )
+    ),
+    dynamicMeshCoeffs_
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict(typeName + "Coeffs")
+    )
+{
+    // Get zones and their interpolation tables for displacement
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    List<Pair<word> > faceZoneToTable
+    (
+        dynamicMeshCoeffs_.lookup("interpolationTables")
+    );
+
+    const faceZoneMesh& fZones = mesh.faceZones();
+
+    times_.setSize(fZones.size());
+    displacements_.setSize(fZones.size());
+
+    forAll(faceZoneToTable, i)
+    {
+        const word& zoneName = faceZoneToTable[i][0];
+        label zoneI = fZones.findZoneID(zoneName);
+
+        if (zoneI == -1)
+        {
+            FatalErrorIn
+            (
+                "displacementInterpolationFvMotionSolver::"
+                "displacementInterpolationFvMotionSolver(const polyMesh&,"
+                "Istream&)"
+            )   << "Cannot find zone " << zoneName << endl
+                << "Valid zones are " << mesh.faceZones().names()
+                << exit(FatalError);
+        }
+
+        const word& tableName = faceZoneToTable[i][1];
+
+        IOList<Tuple2<scalar, vector> > table
+        (
+            IOobject
+            (
+                tableName,
+                mesh.time().constant(),
+                "tables",
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        );
+
+        // Copy table
+        times_[zoneI].setSize(table.size());
+        displacements_[zoneI].setSize(table.size());
+
+        forAll(table, j)
+        {
+            times_[zoneI][j] = table[j].first();
+            displacements_[zoneI][j] = table[j].second();
+        }
+    }
+
+
+
+    // Sort points into bins according to position relative to faceZones
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Done in all three directions.
+
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        // min and max coordinates of all faceZones
+        SortableList<scalar> zoneCoordinates(2*faceZoneToTable.size());
+
+        forAll(faceZoneToTable, i)
+        {
+            const word& zoneName = faceZoneToTable[i][0];
+            label zoneI = fZones.findZoneID(zoneName);
+            const faceZone& fz = fZones[zoneI];
+
+            scalarField fzCoords = fz().localPoints().component(dir);
+
+            zoneCoordinates[2*i] = gMin(fzCoords);
+            zoneCoordinates[2*i+1] = gMax(fzCoords);
+        }
+        zoneCoordinates.sort();
+
+        // Slightly tweak min and max face zone so points sort within
+        zoneCoordinates[0] -= SMALL;
+        zoneCoordinates[zoneCoordinates.size()-1] += SMALL;
+
+        // Check if we have static min and max mesh bounds
+        const scalarField meshCoords = points0_.component(dir);
+
+        scalar minCoord = gMin(meshCoords);
+        scalar maxCoord = gMax(meshCoords);
+
+        if (debug)
+        {
+            Pout<< "direction " << dir << " : "
+                << "mesh ranges from coordinate " << minCoord << " to "
+                << maxCoord << endl;
+        }
+
+        // Make copy of zoneCoordinates; include min and max of mesh
+        // if necessary. Mark min and max with zoneI=-1.
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        labelList& rangeZone = rangeToZone_[dir];
+        labelListList& rangePoints = rangeToPoints_[dir];
+        List<scalarField>& rangeWeights = rangeToWeights_[dir];
+
+        scalarField rangeToCoord(zoneCoordinates.size());
+        rangeZone.setSize(zoneCoordinates.size());
+        label rangeI = 0;
+
+        if (minCoord < zoneCoordinates[0])
+        {
+            label sz = rangeZone.size();
+            rangeToCoord.setSize(sz+1);
+            rangeZone.setSize(sz+1);
+            rangeToCoord[rangeI] = minCoord-SMALL;
+            rangeZone[rangeI] = -1;
+
+            if (debug)
+            {
+                Pout<< "direction " << dir << " : "
+                    << "range " << rangeI << " at coordinate "
+                    << rangeToCoord[rangeI] << " from min of mesh "
+                    << rangeZone[rangeI] << endl;
+            }
+            rangeI = 1;
+        }
+        forAll(zoneCoordinates, i)
+        {
+            rangeToCoord[rangeI] = zoneCoordinates[i];
+            rangeZone[rangeI] = zoneCoordinates.indices()[i]/2;
+
+            if (debug)
+            {
+                Pout<< "direction " << dir << " : "
+                    << "range " << rangeI << " at coordinate "
+                    << rangeToCoord[rangeI]
+                    << " from zone " << rangeZone[rangeI] << endl;
+            }
+            rangeI++;
+        }
+        if (maxCoord > zoneCoordinates[zoneCoordinates.size()-1])
+        {
+            label sz = rangeToCoord.size();
+            rangeToCoord.setSize(sz+1);
+            rangeZone.setSize(sz+1);
+            rangeToCoord[sz] = maxCoord+SMALL;
+            rangeZone[sz] = -1;
+
+            if (debug)
+            {
+                Pout<< "direction " << dir << " : "
+                    << "range " << rangeI << " at coordinate "
+                    << rangeToCoord[sz] << " from max of mesh "
+                    << rangeZone[sz] << endl;
+            }
+        }
+
+
+        // Sort the points
+        // ~~~~~~~~~~~~~~~
+
+        // Count all the points inbetween rangeI and rangeI+1
+        labelList nRangePoints(rangeToCoord.size(), 0);
+
+        forAll(meshCoords, pointI)
+        {
+            label rangeI = findLower(rangeToCoord, meshCoords[pointI]);
+
+            if (rangeI == -1 || rangeI == rangeToCoord.size()-1)
+            {
+                FatalErrorIn
+                (
+                    "displacementInterpolationFvMotionSolver::"
+                    "displacementInterpolationFvMotionSolver"
+                    "(const polyMesh&, Istream&)"
+                )   << "Did not find point " << points0_[pointI]
+                    << " coordinate " << meshCoords[pointI]
+                    << " in ranges " << rangeToCoord
+                    << abort(FatalError);
+            }
+            nRangePoints[rangeI]++;
+        }
+
+        if (debug)
+        {
+            for (label rangeI = 0; rangeI < rangeToCoord.size()-1; rangeI++)
+            {
+                // Get the two zones bounding the range
+                Pout<< "direction " << dir << " : "
+                    << "range from " << rangeToCoord[rangeI]
+                    << " to " << rangeToCoord[rangeI+1]
+                    << " contains " << nRangePoints[rangeI]
+                    << " points." << endl;
+            }
+        }
+
+        // Sort
+        rangePoints.setSize(nRangePoints.size());
+        rangeWeights.setSize(nRangePoints.size());
+        forAll(rangePoints, rangeI)
+        {
+            rangePoints[rangeI].setSize(nRangePoints[rangeI]);
+            rangeWeights[rangeI].setSize(nRangePoints[rangeI]);
+        }
+        nRangePoints = 0;
+        forAll(meshCoords, pointI)
+        {
+            label rangeI = findLower(rangeToCoord, meshCoords[pointI]);
+            label& nPoints = nRangePoints[rangeI];
+            rangePoints[rangeI][nPoints] = pointI;
+            rangeWeights[rangeI][nPoints] =
+                (meshCoords[pointI]-rangeToCoord[rangeI])
+              / (rangeToCoord[rangeI+1]-rangeToCoord[rangeI]);
+            nPoints++;
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::displacementInterpolationFvMotionSolver::~displacementInterpolationFvMotionSolver()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::pointField>
+Foam::displacementInterpolationFvMotionSolver::curPoints() const
+{
+    if (mesh().nPoints() != points0_.size())
+    {
+        FatalErrorIn
+        (
+            "displacementInterpolationFvMotionSolver::curPoints() const"
+        )   << "The number of points in the mesh seems to have changed." << endl
+            << "In constant/polyMesh there are " << points0_.size()
+            << " points; in the current mesh there are " << mesh().nPoints()
+            << " points." << exit(FatalError);
+    }
+
+    tmp<pointField> tcurPoints(new pointField(points0_));
+    pointField& curPoints = tcurPoints();
+
+    // Interpolate the diplacement of the face zones.
+    vectorField zoneDisp(displacements_.size(), vector::zero);
+    forAll(zoneDisp, zoneI)
+    {
+        if (times_[zoneI].size() > 0)
+        {
+            zoneDisp[zoneI] = interpolateXY
+            (
+                mesh().time().value(),
+                times_[zoneI],
+                displacements_[zoneI]
+            );
+        }
+    }
+    if (debug)
+    {
+        Pout<< "Zone displacements:" << zoneDisp << endl;
+    }
+
+
+    // Interpolate the point location
+    for (direction dir = 0; dir < vector::nComponents; dir++)
+    {
+        const labelList& rangeZone = rangeToZone_[dir];
+        const labelListList& rangePoints = rangeToPoints_[dir];
+        const List<scalarField>& rangeWeights = rangeToWeights_[dir];
+
+        for (label rangeI = 0; rangeI < rangeZone.size()-1; rangeI++)
+        {
+            const labelList& rPoints = rangePoints[rangeI];
+            const scalarField& rWeights = rangeWeights[rangeI];
+
+            // Get the two zones bounding the range
+            label minZoneI = rangeZone[rangeI];
+            //vector minDisp =
+            //    (minZoneI == -1 ? vector::zero : zoneDisp[minZoneI]);
+            scalar minDisp = (minZoneI == -1 ? 0.0 : zoneDisp[minZoneI][dir]);
+            label maxZoneI = rangeZone[rangeI+1];
+            //vector maxDisp =
+            //    (maxZoneI == -1 ? vector::zero : zoneDisp[maxZoneI]);
+            scalar maxDisp = (maxZoneI == -1 ? 0.0 : zoneDisp[maxZoneI][dir]);
+
+            forAll(rPoints, i)
+            {
+                label pointI = rPoints[i];
+                scalar w = rWeights[i];
+                //curPoints[pointI] += (1.0-w)*minDisp+w*maxDisp;
+                curPoints[pointI][dir] += (1.0-w)*minDisp+w*maxDisp;
+            }
+        }
+    }
+    return tcurPoints;
+}
+
+
+void Foam::displacementInterpolationFvMotionSolver::updateMesh
+(
+    const mapPolyMesh& mpm
+)
+{
+    fvMotionSolver::updateMesh(mpm);
+
+    // Map points0_. Bit special since we somehow have to come up with
+    // a sensible points0 position for introduced points.
+    // Find out scaling between points0 and current points
+
+    // Get the new points either from the map or the mesh
+    const pointField& points =
+    (
+        mpm.hasMotionPoints()
+      ? mpm.preMotionPoints()
+      : fvMesh_.points()
+    );
+
+    // Note: boundBox does reduce
+    const boundBox bb0(points0_, true);
+    const vector span0(bb0.max()-bb0.min());
+    const boundBox bb(points, true);
+    const vector span(bb.max()-bb.min());
+
+    vector scaleFactors(cmptDivide(span0, span));
+
+    pointField newPoints0(mpm.pointMap().size());
+
+    forAll(newPoints0, pointI)
+    {
+        label oldPointI = mpm.pointMap()[pointI];
+    
+        if (oldPointI >= 0)
+        {
+            label masterPointI = mpm.reversePointMap()[oldPointI];
+
+            if (masterPointI == pointI)
+            {
+                newPoints0[pointI] = points0_[oldPointI];
+            }
+            else
+            {
+                // New point. Assume motion is scaling.
+                newPoints0[pointI] =
+                    points0_[oldPointI]
+                  + cmptMultiply
+                    (
+                        scaleFactors,
+                        points[pointI]-points[masterPointI]
+                    );
+            }
+        }
+        else
+        {
+            FatalErrorIn
+            (
+                "displacementLaplacianFvMotionSolver::updateMesh"
+                "(const mapPolyMesh& mpm)"
+            )   << "Cannot work out coordinates of introduced vertices."
+                << " New vertex " << pointI << " at coordinate "
+                << points[pointI] << exit(FatalError);
+        }
+    }
+    points0_.transfer(newPoints0);
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::displacementInterpolationFvMotionSolver
+
+Description
+    Mesh motion solver for an fvMesh. Scales inbetween motion prescribed on
+    faceZones. Works out per point the distance between the bounding
+    face zones (in all three directions) at the start and then every time
+    step
+    - moves the faceZones based on tables
+    - interpolates the displacement of all points based on the
+    faceZone motion.
+
+    Tables are in the constant/tables directory.
+
+    Note: could be a motionSolver - does not use any fvMesh structure.
+
+SourceFiles
+    displacementInterpolationFvMotionSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef displacementInterpolationFvMotionSolver_H
+#define displacementInterpolationFvMotionSolver_H
+
+#include "fvMotionSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+        Class displacementInterpolationFvMotionSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class displacementInterpolationFvMotionSolver
+:
+    public fvMotionSolver
+{
+    // Private data
+
+        //- Reference point field
+        pointField points0_;
+
+        //- Additional settings for motion solver
+        dictionary dynamicMeshCoeffs_;
+
+
+        // Face zone information (note: could pack these to only contain
+        // used zones)
+
+            //- Interpolation table. From faceZone to times.
+            List<scalarField> times_;
+
+            //- Interpolation table. From faceZone to displacements.
+            List<vectorField> displacements_;
+
+        // Range information.
+
+            //- Per direction, per range the index of the lower
+            //  faceZone
+            FixedList<labelList, 3> rangeToZone_;
+
+            //- Per direction, per range the points that are in it
+            FixedList<labelListList, 3> rangeToPoints_;
+
+            //- Per direction, per range the weight of the points relative
+            //  to this and the next range.
+            FixedList<List<scalarField>, 3> rangeToWeights_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        displacementInterpolationFvMotionSolver
+        (
+            const displacementInterpolationFvMotionSolver&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const displacementInterpolationFvMotionSolver&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("displacementInterpolation");
+
+
+    // Constructors
+
+        displacementInterpolationFvMotionSolver
+        (
+            const polyMesh& mesh,
+            Istream& msData
+        );
+
+
+    // Destructor
+
+        ~displacementInterpolationFvMotionSolver();
+
+
+    // Member Functions
+
+        //- Return reference to the reference field
+        const pointField& points0() const
+        {
+            return points0_;
+        }
+
+        //- Return point location obtained from the current motion field
+        virtual tmp<pointField> curPoints() const;
+
+        //- Solve for motion
+        virtual void solve()
+        {}
+
+        //- Update topology
+        virtual void updateMesh(const mapPolyMesh&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //