ENH: timeVaryingMappedFixedValue: add perturbation to breat ties

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C

@@ -29,8 +29,8 @@ License
 #include "triSurface.H"
 #include "vector2D.H"
 #include "OFstream.H"
-#include "long.H"
 #include "AverageIOField.H"
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,6 +50,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(p, iF),
     fieldTableName_(iF.name()),
     setAverage_(false),
+    perturb_(0),
     referenceCS_(NULL),
     nearestVertex_(0),
     nearestVertexWeight_(0),
@@ -76,6 +77,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
+    perturb_(ptf.perturb_),
     referenceCS_(NULL),
     nearestVertex_(0),
     nearestVertexWeight_(0),
@@ -101,6 +103,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(p, iF),
     fieldTableName_(iF.name()),
     setAverage_(readBool(dict.lookup("setAverage"))),
+    perturb_(dict.lookupOrDefault("perturb", 1E-5)),
     referenceCS_(NULL),
     nearestVertex_(0),
     nearestVertexWeight_(0),
@@ -135,6 +138,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(ptf),
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
+    perturb_(ptf.perturb_),
     referenceCS_(ptf.referenceCS_),
     nearestVertex_(ptf.nearestVertex_),
     nearestVertexWeight_(ptf.nearestVertexWeight_),
@@ -160,6 +164,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(ptf, iF),
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
+    perturb_(ptf.perturb_),
     referenceCS_(ptf.referenceCS_),
     nearestVertex_(ptf.nearestVertex_),
     nearestVertexWeight_(ptf.nearestVertexWeight_),
@@ -334,16 +339,24 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()
     );
     vectorField& localVertices = tlocalVertices();
 
-    // Shear to avoid degenerate cases
+    const boundBox bb(localVertices, true);
+    const point bbMid(bb.midpoint());
+
+    if (debug)
+    {
+        Info<< "timeVaryingMappedFixedValueFvPatchField :"
+            << " Perturbing points with " << perturb_
+            << " fraction of a random position inside " << bb
+            << " to break any ties on regular meshes."
+            << nl << endl;
+    }
+
+    Random rndGen(123456);
     forAll(localVertices, i)
     {
-        point& pt = localVertices[i];
-        const scalar magPt = mag(pt);
-        const point nptDir = pt/magPt;
-        if (magPt > ROOTVSMALL)
-        {
-            pt += pow(magPt, 1.1 + Foam::sqrt(SMALL))*nptDir;
-        }
+        localVertices[i] +=
+            perturb_
+           *(rndGen.position(bb.min(), bb.max())-bbMid);
     }
 
     // Determine triangulation
@@ -354,6 +367,8 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()
         localVertices2D[i][1] = localVertices[i][1];
     }
 
+    triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
+
     tmp<pointField> tlocalFaceCentres
     (
         referenceCS().localPosition
@@ -361,40 +376,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()
             this->patch().patch().faceCentres()
         )
     );
-
-    pointField& localFaceCentres = tlocalFaceCentres();
-
-    // Shear to avoid degenerate cases
-    forAll(localFaceCentres, i)
-    {
-        point& pt = localFaceCentres[i];
-        const scalar magPt = mag(pt);
-        const point nptDir = pt/magPt;
-        if (magPt > ROOTVSMALL)
-        {
-            pt += pow(magPt, 1.1 + Foam::sqrt(SMALL))*nptDir;
-        }
-    }
-
-    if (debug)
-    {
-        OFstream str
-        (
-            this->db().time().path()/this->patch().name()
-          + "_localFaceCentres.obj"
-        );
-        Pout<< "readSamplePoints :"
-            << " Dumping face centres to " << str.name() << endl;
-
-        forAll(localFaceCentres, i)
-        {
-            const point& p = localFaceCentres[i];
-            str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
-        }
-    }
-
-
-    triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
+    const pointField& localFaceCentres = tlocalFaceCentres();
 
     if (debug)
     {
@@ -786,6 +768,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
     os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl;
+    os.writeKeyword("peturb") << perturb_ << token::END_STATEMENT << nl;
 
     if (fieldTableName_ != this->dimensionedInternalField().name())
     {

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,36 @@ Class
 Description
     Foam::timeVaryingMappedFixedValueFvPatchField
 
+    Interpolates from a set of supplied points in space and time. Supplied
+    data in constant/boundaryData/<patchname>:
+        - points : pointField with locations
+        - ddd    : supplied values at time ddd
+    Points need to be more or less on a plane since get triangulated in 2D.
+
+    At startup this bc does the triangulation and determines linear
+    interpolation (triangle it is in and weights to the 3 vertices)
+    for every face centre. Interpolates linearly inbetween times.
+
+    @verbatim
+        inlet
+        {
+            type            timeVaryingMappedFixedValue;
+
+            // Maintain average to that of the supplied values
+            setAverage      false;
+
+            // Optional: change perturbation (default 1E-5) to avoid any ties
+            // in triangulating regular geometries.
+            //perturb       0.0;
+
+            // Optional: use name instead of patchname for location of data
+            //fieldTableName samples;
+        }
+    @endverbatim
+
+    Switch on debug flag to have it dump the triangulation (in transformed
+    space) and transform face centres.
+
 SourceFiles
     timeVaryingMappedFixedValueFvPatchField.C
 
@@ -62,6 +92,9 @@ class timeVaryingMappedFixedValueFvPatchField
         //- If true adjust the mapped field to maintain average value
         bool setAverage_;
 
+        //- Fraction of perturbation (fraction of bounding box) to add
+        scalar perturb_;
+
         //- Coordinate system
         autoPtr<coordinateSystem> referenceCS_;