FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

FitData.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "FitData.H"
00027 #include <finiteVolume/surfaceFields.H>
00028 #include <finiteVolume/volFields.H>
00029 #include <OpenFOAM/SVD.H>
00030 
00031 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
00032 
00033 template<class Form, class ExtendedStencil, class Polynomial>
00034 Foam::FitData<Form, ExtendedStencil, Polynomial>::FitData
00035 (
00036     const fvMesh& mesh,
00037     const ExtendedStencil& stencil,
00038     const bool linearCorrection,
00039     const scalar linearLimitFactor,
00040     const scalar centralWeight
00041 )
00042 :
00043     MeshObject<fvMesh, Form>(mesh),
00044     stencil_(stencil),
00045     linearCorrection_(linearCorrection),
00046     linearLimitFactor_(linearLimitFactor),
00047     centralWeight_(centralWeight),
00048 #   ifdef SPHERICAL_GEOMETRY
00049     dim_(2),
00050 #   else
00051     dim_(mesh.nGeometricD()),
00052 #   endif
00053     minSize_(Polynomial::nTerms(dim_))
00054 {
00055     // Check input
00056     if (linearLimitFactor <= SMALL || linearLimitFactor > 3)
00057     {
00058         FatalErrorIn("FitData<Polynomial>::FitData(..)")
00059             << "linearLimitFactor requested = " << linearLimitFactor
00060             << " should be between zero and 3"
00061             << exit(FatalError);
00062     }
00063 }
00064 
00065 
00066 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00067 
00068 template<class FitDataType, class ExtendedStencil, class Polynomial>
00069 void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs
00070 (
00071     vector& idir,        // value changed in return
00072     vector& jdir,        // value changed in return
00073     vector& kdir,        // value changed in return
00074     const label facei
00075 )
00076 {
00077     const fvMesh& mesh = this->mesh();
00078 
00079     idir = mesh.faceAreas()[facei];
00080     idir /= mag(idir);
00081 
00082 #   ifndef SPHERICAL_GEOMETRY
00083     if (mesh.nGeometricD() <= 2) // find the normal direction
00084     {
00085         if (mesh.geometricD()[0] == -1)
00086         {
00087             kdir = vector(1, 0, 0);
00088         }
00089         else if (mesh.geometricD()[1] == -1)
00090         {
00091             kdir = vector(0, 1, 0);
00092         }
00093         else
00094         {
00095             kdir = vector(0, 0, 1);
00096         }
00097     }
00098     else // 3D so find a direction in the plane of the face
00099     {
00100         const face& f = mesh.faces()[facei];
00101         kdir = mesh.points()[f[0]] - mesh.faceCentres()[facei];
00102     }
00103 #   else
00104     // Spherical geometry so kdir is the radial direction
00105     kdir = mesh.faceCentres()[facei];
00106 #   endif
00107 
00108     if (mesh.nGeometricD() == 3)
00109     {
00110         // Remove the idir component from kdir and normalise
00111         kdir -= (idir & kdir)*idir;
00112 
00113         scalar magk = mag(kdir);
00114 
00115         if (magk < SMALL)
00116         {
00117             FatalErrorIn("findFaceDirs(..)") << " calculated kdir = zero"
00118                 << exit(FatalError);
00119         }
00120         else
00121         {
00122             kdir /= magk;
00123         }
00124     }
00125 
00126     jdir = kdir ^ idir;
00127 }
00128 
00129 
00130 template<class FitDataType, class ExtendedStencil, class Polynomial>
00131 void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
00132 (
00133     scalarList& coeffsi,
00134     const List<point>& C,
00135     const scalar wLin,
00136     const label facei
00137 )
00138 {
00139     vector idir(1,0,0);
00140     vector jdir(0,1,0);
00141     vector kdir(0,0,1);
00142     findFaceDirs(idir, jdir, kdir, facei);
00143 
00144     // Setup the point weights
00145     scalarList wts(C.size(), scalar(1));
00146     wts[0] = centralWeight_;
00147     if (linearCorrection_)
00148     {
00149         wts[1] = centralWeight_;
00150     }
00151 
00152     // Reference point
00153     point p0 = this->mesh().faceCentres()[facei];
00154 
00155     // Info << "Face " << facei << " at " << p0 << " stencil points at:\n"
00156     //     << C - p0 << endl;
00157 
00158     // p0 -> p vector in the face-local coordinate system
00159     vector d;
00160 
00161     // Local coordinate scaling
00162     scalar scale = 1;
00163 
00164     // Matrix of the polynomial components
00165     scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
00166 
00167     for(label ip = 0; ip < C.size(); ip++)
00168     {
00169         const point& p = C[ip];
00170 
00171         d.x() = (p - p0)&idir;
00172         d.y() = (p - p0)&jdir;
00173 #       ifndef SPHERICAL_GEOMETRY
00174         d.z() = (p - p0)&kdir;
00175 #       else
00176         d.z() = mag(p) - mag(p0);
00177 #       endif
00178 
00179         if (ip == 0)
00180         {
00181             scale = cmptMax(cmptMag((d)));
00182         }
00183 
00184         // Scale the radius vector
00185         d /= scale;
00186 
00187         Polynomial::addCoeffs
00188         (
00189             B[ip],
00190             d,
00191             wts[ip],
00192             dim_
00193         );
00194     }
00195 
00196     // Additional weighting for constant and linear terms
00197     for(label i = 0; i < B.n(); i++)
00198     {
00199         B[i][0] *= wts[0];
00200         B[i][1] *= wts[0];
00201     }
00202 
00203     // Set the fit
00204     label stencilSize = C.size();
00205     coeffsi.setSize(stencilSize);
00206 
00207     bool goodFit = false;
00208     for(int iIt = 0; iIt < 8 && !goodFit; iIt++)
00209     {
00210         SVD svd(B, SMALL);
00211 
00212         scalar maxCoeff = 0;
00213         label maxCoeffi = 0;
00214 
00215         for(label i=0; i<stencilSize; i++)
00216         {
00217             coeffsi[i] = wts[0]*wts[i]*svd.VSinvUt()[0][i];
00218             if (mag(coeffsi[i]) > maxCoeff)
00219             {
00220                 maxCoeff = mag(coeffsi[i]);
00221                 maxCoeffi = i;
00222             }
00223         }
00224 
00225         if (linearCorrection_)
00226         {
00227             goodFit =
00228                 (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin)
00229              && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin))
00230              && maxCoeffi <= 1;
00231         }
00232         else
00233         {
00234             // Upwind: weight on face is 1.
00235             goodFit =
00236                 (mag(coeffsi[0] - 1.0) < linearLimitFactor_*1.0)
00237              && maxCoeffi <= 1;
00238         }
00239 
00240         // if (goodFit && iIt > 0)
00241         // {
00242             // Info << "FitData<Polynomial>::calcFit"
00243             //     << "(const List<point>& C, const label facei" << nl
00244             //     << "Can now fit face " << facei << " iteration " << iIt
00245             //     << " with sum of weights " << sum(coeffsi) << nl
00246             //     << "    Weights " << coeffsi << nl
00247             //     << "    Linear weights " << wLin << " " << 1 - wLin << nl
00248             //     << "    sing vals " << svd.S() << endl;
00249         // }
00250 
00251         if (!goodFit) // (not good fit so increase weight in the centre and weight
00252                       //  for constant and linear terms)
00253         {
00254             // if (iIt == 7)
00255             // {
00256             //     WarningIn
00257             //     (
00258             //         "FitData<Polynomial>::calcFit"
00259             //         "(const List<point>& C, const label facei"
00260             //     )   << "Cannot fit face " << facei << " iteration " << iIt
00261             //         << " with sum of weights " << sum(coeffsi) << nl
00262             //         << "    Weights " << coeffsi << nl
00263             //         << "    Linear weights " << wLin << " " << 1 - wLin << nl
00264             //         << "    sing vals " << svd.S() << endl;
00265             // }
00266 
00267             wts[0] *= 10;
00268             if (linearCorrection_)
00269             {
00270                 wts[1] *= 10;
00271             }
00272 
00273             for(label j = 0; j < B.m(); j++)
00274             {
00275                 B[0][j] *= 10;
00276                 B[1][j] *= 10;
00277             }
00278 
00279             for(label i = 0; i < B.n(); i++)
00280             {
00281                 B[i][0] *= 10;
00282                 B[i][1] *= 10;
00283             }
00284         }
00285     }
00286 
00287     if (goodFit)
00288     {
00289         if (linearCorrection_)
00290         {
00291             // Remove the uncorrected linear coefficients
00292             coeffsi[0] -= wLin;
00293             coeffsi[1] -= 1 - wLin;
00294         }
00295         else
00296         {
00297             // Remove the uncorrected upwind coefficients
00298             coeffsi[0] -= 1.0;
00299         }
00300     }
00301     else
00302     {
00303         // if (debug)
00304         // {
00305             WarningIn
00306             (
00307                 "FitData<Polynomial>::calcFit(..)"
00308             )   << "Could not fit face " << facei
00309                 << "    Weights = " << coeffsi
00310                 << ", reverting to linear." << nl
00311                 << "    Linear weights " << wLin << " " << 1 - wLin << endl;
00312         // }
00313 
00314         coeffsi = 0;
00315     }
00316 }
00317 
00318 
00319 template<class FitDataType, class ExtendedStencil, class Polynomial>
00320 bool Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::movePoints()
00321 {
00322     calcFit();
00323     return true;
00324 }
00325 
00326 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines