fvOptions: verticalDamping: Allow multiple directions of ramping

The ramp function used to graduate the vertical damping force can now be
applied along a number of paths, rather than just one. The keywords
"origins" and "directions" can be used to define a list of paths.

    verticalDamping1
    {
        type            verticalDamping;

        origins         ((1200 0 0) (1200 100 0) (1200 -100 0));
        directions      ((1 0 0) (0 1 0) (0 -1 0));

        // ...
    }

The ramping function will be calculated along each of the paths defined
by the origin-direction pair, and the maximum of the calculated values
will be used.

The "origin" and "direction" keywords can still be used with non-list
values.

This work was supported by Jan Kaufmann and Jan Oberhagemann at DNV GL.

src/fvOptions/sources/derived/verticalDamping/verticalDamping.C

@@ -30,6 +30,7 @@ License
 #include "meshTools.H"
 #include "Function1.H"
 #include "uniformDimensionedFields.H"
+#include "zeroGradientFvPatchField.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -61,12 +62,13 @@ void Foam::fv::verticalDamping::add
     const DimensionedField<scalar, volMesh>& V = mesh_.V();
 
     // Calculate the scale
-    const scalarField s
-    (
-        ramp_.valid()
-      ? ramp_->value((mesh_.cellCentres() - origin_) & direction_)
-      : tmp<scalarField>(new scalarField(mesh_.nCells(), 1))
-    );
+    scalarField s(mesh_.nCells(), ramp_.valid() ? 0 : 1);
+    forAll(origins_, i)
+    {
+        const vectorField& c = mesh_.cellCentres();
+        const scalarField x((c - origins_[i]) & directions_[i]);
+        s = max(s, ramp_->value(x));
+    }
 
     // Check dimensions
     eqn.dimensions()
@@ -85,6 +87,25 @@ void Foam::fv::verticalDamping::add
         const label c = cells_[i];
         eqn.source()[c] += force[i];
     }
+
+    // Write out the force coefficient for debugging
+    if (debug && mesh_.time().writeTime())
+    {
+        volScalarField forceCoeff
+        (
+            IOobject
+            (
+                type() + ":forceCoeff",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            lambda_*mag(g),
+            zeroGradientFvPatchField<scalar>::typeName
+        );
+        forceCoeff.primitiveFieldRef() *= s;
+        forceCoeff.write();
+    }
 }
 
 
@@ -101,8 +122,8 @@ Foam::fv::verticalDamping::verticalDamping
     cellSetOption(name, modelType, dict, mesh),
     lambda_("lambda", dimless/dimTime, coeffs_.lookup("lambda")),
     ramp_(),
-    origin_(),
-    direction_()
+    origins_(),
+    directions_()
 {
     read(dict);
 }
@@ -156,22 +177,67 @@ bool Foam::fv::verticalDamping::read(const dictionary& dict)
             );
 
         const bool foundRamp = coeffs_.found("ramp");
-        const bool foundOrigin = coeffs_.found("origin");
-        const bool foundDirection = coeffs_.found("direction");
-        if (foundRamp && foundOrigin && foundDirection)
+        const bool foundOgn = coeffs_.found("origin");
+        const bool foundDir = coeffs_.found("direction");
+        const bool foundOgns = coeffs_.found("origins");
+        const bool foundDirs = coeffs_.found("directions");
+        const bool foundAll =
+            foundRamp
+         && (
+                (foundOgn && foundDir && !foundOgns && !foundDirs)
+             || (!foundOgn && !foundDir && foundOgns && foundDirs)
+            );
+         const bool foundAny =
+            foundRamp || foundOgn || foundDir || foundOgns || foundDirs;
+
+        if (!foundAll)
+        {
+            ramp_ = autoPtr<Function1<scalar>>();
+            origins_.clear();
+            directions_.clear();
+        }
+
+        if (foundAll)
         {
             ramp_ = Function1<scalar>::New("ramp", coeffs_);
-            coeffs_.lookup("origin") >> origin_;
-            coeffs_.lookup("direction") >> direction_;
-            direction_ /= mag(direction_);
+            if (foundOgn)
+            {
+                origins_.setSize(1);
+                directions_.setSize(1);
+                coeffs_.lookup("origin") >> origins_.last();
+                coeffs_.lookup("direction") >> directions_.last();
+            }
+            else
+            {
+                coeffs_.lookup("origins") >> origins_;
+                coeffs_.lookup("directions") >> directions_;
+
+                if
+                (
+                    origins_.size() == 0
+                 || directions_.size() == 0
+                 || origins_.size() != directions_.size()
+                )
+                {
+                    FatalErrorInFunction
+                        << "The same, non-zero number of origins and "
+                        << "directions must be provided" << exit(FatalError);
+                }
+            }
+            forAll(directions_, i)
+            {
+                directions_[i] /= mag(directions_[i]);
+            }
         }
-        else if (foundRamp || foundOrigin || foundDirection)
+
+        if (!foundAll && foundAny)
         {
             WarningInFunction
                 << "The ramping specification is incomplete. \"ramp\", "
-                << "\"origin\" and \"direction\", must all be specified in "
-                << "order to ramp the damping. The damping will be applied "
-                << "uniformly across the cell set." << endl;
+                << "\"origin\" and \"direction\" (or \"origins\" and "
+                << "\"directions\"), must all be specified in order to ramp "
+                << "the damping. The damping will be applied uniformly across "
+                << "the cell set." << endl << endl;
         }
 
         fieldNames_ = wordList(1, coeffs_.lookupOrDefault<word>("U", "U"));

src/fvOptions/sources/derived/verticalDamping/verticalDamping.H

@@ -75,8 +75,14 @@ Usage
 
         selectionMode   all;
 
+        // Define the line along which to apply the graduation
         origin          (1200 0 0);
         direction       (1 0 0);
+
+        // Or, define multiple lines
+        //origins         ((1200 0 0) (1200 -300 0) (1200 300 0));
+        //directions      ((1 0 0) (0 -1 0) (0 1 0));
+
         ramp
         {
             type        halfCosineRamp;
@@ -129,11 +135,11 @@ private:
         //- The ramping function
         autoPtr<Function1<scalar>> ramp_;
 
-        //- Origin of the ramping coordinate
-        vector origin_;
+        //- Origins of the ramping coordinate
+        vectorField origins_;
 
-        //- Direction of increasing ramping coordinate
-        vector direction_;
+        //- Directions of increasing ramping coordinate
+        vectorField directions_;
 
 
     // Private Member Functions