src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "FitData.H"
  27 #include "surfaceFields.H"
  28 #include "volFields.H"
  29 #include "SVD.H"
  30
  31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
  32
  33 template<class Form, class ExtendedStencil, class Polynomial>
  34 Foam::FitData<Form, ExtendedStencil, Polynomial>::FitData
  35 (
  36     const fvMesh& mesh,
  37     const ExtendedStencil& stencil,
  38     const bool linearCorrection,
  39     const scalar linearLimitFactor,
  40     const scalar centralWeight
  41 )
  42 :
  43     MeshObject<fvMesh, Form>(mesh),
  44     stencil_(stencil),
  45     linearCorrection_(linearCorrection),
  46     linearLimitFactor_(linearLimitFactor),
  47     centralWeight_(centralWeight),
  48 #   ifdef SPHERICAL_GEOMETRY
  49     dim_(2),
  50 #   else
  51     dim_(mesh.nGeometricD()),
  52 #   endif
  53     minSize_(Polynomial::nTerms(dim_))
  54 {
  55     // Check input
  56     if (linearLimitFactor <= SMALL || linearLimitFactor > 3)
  57     {
  58         FatalErrorIn("FitData<Polynomial>::FitData(..)")
  59             << "linearLimitFactor requested = " << linearLimitFactor
  60             << " should be between zero and 3"
  61             << exit(FatalError);
  62     }
  63 }
  64
  65
  66 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  67
  68 template<class FitDataType, class ExtendedStencil, class Polynomial>
  69 void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
  70 (
  71     vector& idir,        // value changed in return
  72     vector& jdir,        // value changed in return
  73     vector& kdir,        // value changed in return
  74     const label facei
  75 )
  76 {
  77     const fvMesh& mesh = this->mesh();
  78
  79     idir = mesh.faceAreas()[facei];
  80     idir /= mag(idir);
  81
  82 #   ifndef SPHERICAL_GEOMETRY
  83     if (mesh.nGeometricD() <= 2) // find the normal direction
  84     {
  85         if (mesh.geometricD()[0] == -1)
  86         {
  87             kdir = vector(1, 0, 0);
  88         }
  89         else if (mesh.geometricD()[1] == -1)
  90         {
  91             kdir = vector(0, 1, 0);
  92         }
  93         else
  94         {
  95             kdir = vector(0, 0, 1);
  96         }
  97     }
  98     else // 3D so find a direction in the plane of the face
  99     {
 100         const face& f = mesh.faces()[facei];
 101         kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei];
 102     }
 103 #   else
 104     // Spherical geometry so kdir is the radial direction
 105     kdir = mesh.faceCentres()[facei];
 106 #   endif
 107
 108     if (mesh.nGeometricD() == 3)
 109     {
 110         // Remove the idir component from kdir and normalise
 111         kdir -= (idir & kdir)*idir;
 112
 113         scalar magk = mag(kdir);
 114
 115         if (magk < SMALL)
 116         {
 117             FatalErrorIn("findFaceDirs(..)") << " calculated kdir = zero"
 118                 << exit(FatalError);
 119         }
 120         else
 121         {
 122             kdir /= magk;
 123         }
 124     }
 125
 126     jdir = kdir ^ idir;
 127 }
 128
 129
 130 template<class FitDataType, class ExtendedStencil, class Polynomial>
 131 void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
 132 (
 133     scalarList& coeffsi,
 134     const List<point>& C,
 135     const scalar wLin,
 136     const label facei
 137 )
 138 {
 139     vector idir(1,0,0);
 140     vector jdir(0,1,0);
 141     vector kdir(0,0,1);
 142     findFaceDirs(idir, jdir, kdir, facei);
 143
 144     // Setup the point weights
 145     scalarList wts(C.size(), scalar(1));
 146     wts[0] = centralWeight_;
 147     if (linearCorrection_)
 148     {
 149         wts[1] = centralWeight_;
 150     }
 151
 152     // Reference point
 153     point p0 = this->mesh().faceCentres()[facei];
 154
 155     // Info<< "Face " << facei << " at " << p0 << " stencil points at:\n"
 156     //     << C - p0 << endl;
 157
 158     // p0 -> p vector in the face-local coordinate system
 159     vector d;
 160
 161     // Local coordinate scaling
 162     scalar scale = 1;
 163
 164     // Matrix of the polynomial components
 165     scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
 166
 167     forAll(C, ip)
 168     {
 169         const point& p = C[ip];
 170
 171         d.x() = (p - p0)&idir;
 172         d.y() = (p - p0)&jdir;
 173 #       ifndef SPHERICAL_GEOMETRY
 174         d.z() = (p - p0)&kdir;
 175 #       else
 176         d.z() = mag(p) - mag(p0);
 177 #       endif
 178
 179         if (ip == 0)
 180         {
 181             scale = cmptMax(cmptMag((d)));
 182         }
 183
 184         // Scale the radius vector
 185         d /= scale;
 186
 187         Polynomial::addCoeffs
 188         (
 189             B[ip],
 190             d,
 191             wts[ip],
 192             dim_
 193         );
 194     }
 195
 196     // Additional weighting for constant and linear terms
 197     for (label i = 0; i < B.n(); i++)
 198     {
 199         B[i][0] *= wts[0];
 200         B[i][1] *= wts[0];
 201     }
 202
 203     // Set the fit
 204     label stencilSize = C.size();
 205     coeffsi.setSize(stencilSize);
 206
 207     bool goodFit = false;
 208     for (int iIt = 0; iIt < 8 && !goodFit; iIt++)
 209     {
 210         SVD svd(B, SMALL);
 211
 212         scalar maxCoeff = 0;
 213         label maxCoeffi = 0;
 214
 215         for (label i=0; i<stencilSize; i++)
 216         {
 217             coeffsi[i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
 218             if (mag(coeffsi[i]) > maxCoeff)
 219             {
 220                 maxCoeff = mag(coeffsi[i]);
 221                 maxCoeffi = i;
 222             }
 223         }
 224
 225         if (linearCorrection_)
 226         {
 227             goodFit =
 228                 (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin)
 229              && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin))
 230              && maxCoeffi <= 1;
 231         }
 232         else
 233         {
 234             // Upwind: weight on face is 1.
 235             goodFit =
 236                 (mag(coeffsi[0] - 1.0) < linearLimitFactor_*1.0)
 237              && maxCoeffi <= 1;
 238         }
 239
 240         // if (goodFit && iIt > 0)
 241         // {
 242             // Info<< "FitData<Polynomial>::calcFit"
 243             //     << "(const List<point>& C, const label facei" << nl
 244             //     << "Can now fit face " << facei << " iteration " << iIt
 245             //     << " with sum of weights " << sum(coeffsi) << nl
 246             //     << "    Weights " << coeffsi << nl
 247             //     << "    Linear weights " << wLin << " " << 1 - wLin << nl
 248             //     << "    sing vals " << svd.S() << endl;
 249         // }
 250
 251         if (!goodFit) // (not good fit so increase weight in the centre and
 252                       //  weight for constant and linear terms)
 253         {
 254             // if (iIt == 7)
 255             // {
 256             //     WarningIn
 257             //     (
 258             //         "FitData<Polynomial>::calcFit"
 259             //         "(const List<point>& C, const label facei"
 260             //     )   << "Cannot fit face " << facei << " iteration " << iIt
 261             //         << " with sum of weights " << sum(coeffsi) << nl
 262             //         << "    Weights " << coeffsi << nl
 263             //         << "    Linear weights " << wLin << " " << 1 - wLin << nl
 264             //         << "    sing vals " << svd.S() << endl;
 265             // }
 266
 267             wts[0] *= 10;
 268             if (linearCorrection_)
 269             {
 270                 wts[1] *= 10;
 271             }
 272
 273             for (label j = 0; j < B.m(); j++)
 274             {
 275                 B[0][j] *= 10;
 276                 B[1][j] *= 10;
 277             }
 278
 279             for (label i = 0; i < B.n(); i++)
 280             {
 281                 B[i][0] *= 10;
 282                 B[i][1] *= 10;
 283             }
 284         }
 285     }
 286
 287     if (goodFit)
 288     {
 289         if (linearCorrection_)
 290         {
 291             // Remove the uncorrected linear coefficients
 292             coeffsi[0] -= wLin;
 293             coeffsi[1] -= 1 - wLin;
 294         }
 295         else
 296         {
 297             // Remove the uncorrected upwind coefficients
 298             coeffsi[0] -= 1.0;
 299         }
 300     }
 301     else
 302     {
 303         // if (debug)
 304         // {
 305             WarningIn
 306             (
 307                 "FitData<Polynomial>::calcFit(..)"
 308             )   << "Could not fit face " << facei
 309                 << "    Weights = " << coeffsi
 310                 << ", reverting to linear." << nl
 311                 << "    Linear weights " << wLin << " " << 1 - wLin << endl;
 312         // }
 313
 314         coeffsi = 0;
 315     }
 316 }
 317
 318
 319 template<class FitDataType, class ExtendedStencil, class Polynomial>
 320 bool Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::movePoints()
 321 {
 322     calcFit();
 323     return true;
 324 }
 325
 326 // ************************************************************************* //