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