src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.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 Description
  26     Cell to face interpolation scheme. Included in fvMesh.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "fvMesh.H"
  31 #include "volFields.H"
  32 #include "surfaceFields.H"
  33 #include "demandDrivenData.H"
  34 #include "coupledFvPatch.H"
  35 #include "mathematicalConstants.H"
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 defineTypeNameAndDebug(surfaceInterpolation, 0);
  45
  46
  47 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
  48
  49 void surfaceInterpolation::clearOut()
  50 {
  51     deleteDemandDrivenData(weightingFactors_);
  52     deleteDemandDrivenData(differenceFactors_);
  53     deleteDemandDrivenData(correctionVectors_);
  54 }
  55
  56
  57 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
  58
  59 surfaceInterpolation::surfaceInterpolation(const fvMesh& fvm)
  60 :
  61     fvSchemes(static_cast<const objectRegistry&>(fvm)),
  62     fvSolution(static_cast<const objectRegistry&>(fvm)),
  63     mesh_(fvm),
  64     weightingFactors_(NULL),
  65     differenceFactors_(NULL),
  66     orthogonal_(false),
  67     correctionVectors_(NULL)
  68 {}
  69
  70
  71 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
  72
  73 surfaceInterpolation::~surfaceInterpolation()
  74 {
  75     clearOut();
  76 }
  77
  78
  79 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  80
  81 const surfaceScalarField& surfaceInterpolation::weights() const
  82 {
  83     if (!weightingFactors_)
  84     {
  85         makeWeights();
  86     }
  87
  88     return (*weightingFactors_);
  89 }
  90
  91
  92 const surfaceScalarField& surfaceInterpolation::deltaCoeffs() const
  93 {
  94     if (!differenceFactors_)
  95     {
  96         makeDeltaCoeffs();
  97     }
  98
  99     return (*differenceFactors_);
 100 }
 101
 102
 103 bool surfaceInterpolation::orthogonal() const
 104 {
 105     if (orthogonal_ == false && !correctionVectors_)
 106     {
 107         makeCorrectionVectors();
 108     }
 109
 110     return orthogonal_;
 111 }
 112
 113
 114 const surfaceVectorField& surfaceInterpolation::correctionVectors() const
 115 {
 116     if (orthogonal())
 117     {
 118         FatalErrorIn("surfaceInterpolation::correctionVectors()")
 119             << "cannot return correctionVectors; mesh is orthogonal"
 120             << abort(FatalError);
 121     }
 122
 123     return (*correctionVectors_);
 124 }
 125
 126
 127 // Do what is neccessary if the mesh has moved
 128 bool surfaceInterpolation::movePoints()
 129 {
 130     deleteDemandDrivenData(weightingFactors_);
 131     deleteDemandDrivenData(differenceFactors_);
 132
 133     orthogonal_ = false;
 134     deleteDemandDrivenData(correctionVectors_);
 135
 136     return true;
 137 }
 138
 139
 140 void surfaceInterpolation::makeWeights() const
 141 {
 142     if (debug)
 143     {
 144         Info<< "surfaceInterpolation::makeWeights() : "
 145             << "Constructing weighting factors for face interpolation"
 146             << endl;
 147     }
 148
 149
 150     weightingFactors_ = new surfaceScalarField
 151     (
 152         IOobject
 153         (
 154             "weightingFactors",
 155             mesh_.pointsInstance(),
 156             mesh_
 157         ),
 158         mesh_,
 159         dimless
 160     );
 161     surfaceScalarField& weightingFactors = *weightingFactors_;
 162
 163
 164     // Set local references to mesh data
 165     // (note that we should not use fvMesh sliced fields at this point yet
 166     //  since this causes a loop when generating weighting factors in
 167     //  coupledFvPatchField evaluation phase)
 168     const unallocLabelList& owner = mesh_.owner();
 169     const unallocLabelList& neighbour = mesh_.neighbour();
 170
 171     const vectorField& Cf = mesh_.faceCentres();
 172     const vectorField& C = mesh_.cellCentres();
 173     const vectorField& Sf = mesh_.faceAreas();
 174
 175     // ... and reference to the internal field of the weighting factors
 176     scalarField& w = weightingFactors.internalField();
 177
 178     forAll(owner, facei)
 179     {
 180         // Note: mag in the dot-product.
 181         // For all valid meshes, the non-orthogonality will be less that
 182         // 90 deg and the dot-product will be positive.  For invalid
 183         // meshes (d & s <= 0), this will stabilise the calculation
 184         // but the result will be poor.
 185         scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]]));
 186         scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei]));
 187         w[facei] = SfdNei/(SfdOwn + SfdNei);
 188     }
 189
 190     forAll(mesh_.boundary(), patchi)
 191     {
 192         mesh_.boundary()[patchi].makeWeights
 193         (
 194             weightingFactors.boundaryField()[patchi]
 195         );
 196     }
 197
 198
 199     if (debug)
 200     {
 201         Info<< "surfaceInterpolation::makeWeights() : "
 202             << "Finished constructing weighting factors for face interpolation"
 203             << endl;
 204     }
 205 }
 206
 207
 208 void surfaceInterpolation::makeDeltaCoeffs() const
 209 {
 210     if (debug)
 211     {
 212         Info<< "surfaceInterpolation::makeDeltaCoeffs() : "
 213             << "Constructing differencing factors array for face gradient"
 214             << endl;
 215     }
 216
 217     // Force the construction of the weighting factors
 218     // needed to make sure deltaCoeffs are calculated for parallel runs.
 219     weights();
 220
 221     differenceFactors_ = new surfaceScalarField
 222     (
 223         IOobject
 224         (
 225             "differenceFactors_",
 226             mesh_.pointsInstance(),
 227             mesh_
 228         ),
 229         mesh_,
 230         dimless/dimLength
 231     );
 232     surfaceScalarField& DeltaCoeffs = *differenceFactors_;
 233
 234
 235     // Set local references to mesh data
 236     const volVectorField& C = mesh_.C();
 237     const unallocLabelList& owner = mesh_.owner();
 238     const unallocLabelList& neighbour = mesh_.neighbour();
 239     const surfaceVectorField& Sf = mesh_.Sf();
 240     const surfaceScalarField& magSf = mesh_.magSf();
 241
 242     forAll(owner, facei)
 243     {
 244         vector delta = C[neighbour[facei]] - C[owner[facei]];
 245         vector unitArea = Sf[facei]/magSf[facei];
 246
 247         // Standard cell-centre distance form
 248         //DeltaCoeffs[facei] = (unitArea & delta)/magSqr(delta);
 249
 250         // Slightly under-relaxed form
 251         //DeltaCoeffs[facei] = 1.0/mag(delta);
 252
 253         // More under-relaxed form
 254         //DeltaCoeffs[facei] = 1.0/(mag(unitArea & delta) + VSMALL);
 255
 256         // Stabilised form for bad meshes
 257         DeltaCoeffs[facei] = 1.0/max(unitArea & delta, 0.05*mag(delta));
 258     }
 259
 260     forAll(DeltaCoeffs.boundaryField(), patchi)
 261     {
 262         mesh_.boundary()[patchi].makeDeltaCoeffs
 263         (
 264             DeltaCoeffs.boundaryField()[patchi]
 265         );
 266     }
 267 }
 268
 269
 270 void surfaceInterpolation::makeCorrectionVectors() const
 271 {
 272     if (debug)
 273     {
 274         Info<< "surfaceInterpolation::makeCorrectionVectors() : "
 275             << "Constructing non-orthogonal correction vectors"
 276             << endl;
 277     }
 278
 279     correctionVectors_ = new surfaceVectorField
 280     (
 281         IOobject
 282         (
 283             "correctionVectors",
 284             mesh_.pointsInstance(),
 285             mesh_
 286         ),
 287         mesh_,
 288         dimless
 289     );
 290     surfaceVectorField& corrVecs = *correctionVectors_;
 291
 292     // Set local references to mesh data
 293     const volVectorField& C = mesh_.C();
 294     const unallocLabelList& owner = mesh_.owner();
 295     const unallocLabelList& neighbour = mesh_.neighbour();
 296     const surfaceVectorField& Sf = mesh_.Sf();
 297     const surfaceScalarField& magSf = mesh_.magSf();
 298     const surfaceScalarField& DeltaCoeffs = deltaCoeffs();
 299
 300     forAll(owner, facei)
 301     {
 302         vector unitArea = Sf[facei]/magSf[facei];
 303         vector delta = C[neighbour[facei]] - C[owner[facei]];
 304
 305         corrVecs[facei] = unitArea - delta*DeltaCoeffs[facei];
 306     }
 307
 308     // Boundary correction vectors set to zero for boundary patches
 309     // and calculated consistently with internal corrections for
 310     // coupled patches
 311
 312     forAll(corrVecs.boundaryField(), patchi)
 313     {
 314         fvsPatchVectorField& patchcorrVecs = corrVecs.boundaryField()[patchi];
 315
 316         if (!patchcorrVecs.coupled())
 317         {
 318             patchcorrVecs = vector::zero;
 319         }
 320         else
 321         {
 322             const fvsPatchScalarField& patchDeltaCoeffs
 323                 = DeltaCoeffs.boundaryField()[patchi];
 324
 325             const fvPatch& p = patchcorrVecs.patch();
 326
 327             vectorField patchDeltas = mesh_.boundary()[patchi].delta();
 328
 329             forAll(p, patchFacei)
 330             {
 331                 vector unitArea =
 332                     Sf.boundaryField()[patchi][patchFacei]
 333                    /magSf.boundaryField()[patchi][patchFacei];
 334
 335                 const vector& delta = patchDeltas[patchFacei];
 336
 337                 patchcorrVecs[patchFacei] =
 338                     unitArea - delta*patchDeltaCoeffs[patchFacei];
 339             }
 340         }
 341     }
 342
 343     scalar NonOrthogCoeff = 0.0;
 344
 345     // Calculate the non-orthogonality for meshes with 1 face or more
 346     if (returnReduce(magSf.size(), sumOp<label>()) > 0)
 347     {
 348         NonOrthogCoeff =
 349             asin
 350             (
 351                 min
 352                 (
 353                     (sum(magSf*mag(corrVecs))/sum(magSf)).value(),
 354                     1.0
 355                 )
 356             )*180.0/mathematicalConstant::pi;
 357     }
 358
 359     if (debug)
 360     {
 361         Info<< "surfaceInterpolation::makeCorrectionVectors() : "
 362             << "non-orthogonality coefficient = " << NonOrthogCoeff << " deg."
 363             << endl;
 364     }
 365
 366     //NonOrthogCoeff = 0.0;
 367
 368     if (NonOrthogCoeff < 0.1)
 369     {
 370         orthogonal_ = true;
 371         deleteDemandDrivenData(correctionVectors_);
 372     }
 373     else
 374     {
 375         orthogonal_ = false;
 376     }
 377
 378     if (debug)
 379     {
 380         Info<< "surfaceInterpolation::makeCorrectionVectors() : "
 381             << "Finished constructing non-orthogonal correction vectors"
 382             << endl;
 383     }
 384 }
 385
 386
 387 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 388
 389 } // End namespace Foam
 390
 391 // ************************************************************************* //