src/finiteArea/interpolation/edgeInterpolation/edgeInterpolationScheme/edgeInterpolationScheme.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Abstract base class for edge interpolation schemes.
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "edgeInterpolationScheme.H"
  30 #include "areaFields.H"
  31 #include "edgeFields.H"
  32 #include "faPatchFields.H"
  33 #include "coupledFaPatchField.H"
  34 #include "transform.H"
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 namespace Foam
  39 {
  40
  41 // * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
  42
  43 // Return weighting factors for scheme given by name in dictionary
  44 template<class Type>
  45 tmp<edgeInterpolationScheme<Type> > edgeInterpolationScheme<Type>::New
  46 (
  47     const faMesh& mesh,
  48     Istream& schemeData
  49 )
  50 {
  51     if (edgeInterpolation::debug)
  52     {
  53         Info<< "edgeInterpolationScheme<Type>::New(const faMesh&, Istream&)"
  54                " : constructing edgeInterpolationScheme<Type>"
  55             << endl;
  56     }
  57
  58     if (schemeData.eof())
  59     {
  60         FatalIOErrorIn
  61         (
  62             "edgeInterpolationScheme<Type>::New(const faMesh&, Istream&)",
  63             schemeData
  64         )   << "Discretisation scheme not specified"
  65             << endl << endl
  66             << "Valid schemes are :" << endl
  67             << MeshConstructorTablePtr_->sortedToc()
  68             << exit(FatalIOError);
  69     }
  70
  71     word schemeName(schemeData);
  72
  73     typename MeshConstructorTable::iterator constructorIter =
  74         MeshConstructorTablePtr_->find(schemeName);
  75
  76     if (constructorIter == MeshConstructorTablePtr_->end())
  77     {
  78         FatalIOErrorIn
  79         (
  80             "edgeInterpolationScheme<Type>::New(const faMesh&, Istream&)",
  81             schemeData
  82         )   << "Unknown discretisation scheme " << schemeName
  83             << endl << endl
  84             << "Valid schemes are :" << endl
  85             << MeshConstructorTablePtr_->sortedToc()
  86             << exit(FatalIOError);
  87     }
  88
  89     return constructorIter()(mesh, schemeData);
  90 }
  91
  92
  93 // Return weighting factors for scheme given by name in dictionary
  94 template<class Type>
  95 tmp<edgeInterpolationScheme<Type> > edgeInterpolationScheme<Type>::New
  96 (
  97     const faMesh& mesh,
  98     const edgeScalarField& faceFlux,
  99     Istream& schemeData
 100 )
 101 {
 102     if (edgeInterpolation::debug)
 103     {
 104         Info<< "edgeInterpolationScheme<Type>::New"
 105                "(const faMesh&, const edgeScalarField&, Istream&) : "
 106                "constructing edgeInterpolationScheme<Type>"
 107             << endl;
 108     }
 109
 110     if (schemeData.eof())
 111     {
 112         FatalIOErrorIn
 113         (
 114             "edgeInterpolationScheme<Type>::New"
 115             "(const faMesh&, const edgeScalarField&, Istream&)",
 116             schemeData
 117         )   << "Discretisation scheme not specified"
 118             << endl << endl
 119             << "Valid schemes are :" << endl
 120             << MeshConstructorTablePtr_->sortedToc()
 121             << exit(FatalIOError);
 122     }
 123
 124     word schemeName(schemeData);
 125
 126     typename MeshFluxConstructorTable::iterator constructorIter =
 127         MeshFluxConstructorTablePtr_->find(schemeName);
 128
 129     if (constructorIter == MeshFluxConstructorTablePtr_->end())
 130     {
 131         FatalIOErrorIn
 132         (
 133             "edgeInterpolationScheme<Type>::New"
 134             "(const faMesh&, const edgeScalarField&, Istream&)",
 135             schemeData
 136         )   << "Unknown discretisation scheme " << schemeName
 137             << endl << endl
 138             << "Valid schemes are :" << endl
 139             << MeshFluxConstructorTablePtr_->sortedToc()
 140             << exit(FatalIOError);
 141     }
 142
 143     return constructorIter()(mesh, faceFlux, schemeData);
 144 }
 145
 146
 147 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 148
 149 template<class Type>
 150 edgeInterpolationScheme<Type>::~edgeInterpolationScheme()
 151 {}
 152
 153
 154 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 155
 156 //- Return the face-interpolate of the given cell field
 157 //  with the given owner and neighbour weighting factors
 158 template<class Type>
 159 tmp<GeometricField<Type, faePatchField, edgeMesh> >
 160 edgeInterpolationScheme<Type>::interpolate
 161 (
 162     const GeometricField<Type, faPatchField, areaMesh>& vf,
 163     const tmp<edgeScalarField>& tlambdas,
 164     const tmp<edgeScalarField>& tys
 165 )
 166 {
 167     if (edgeInterpolation::debug)
 168     {
 169         Info<< "edgeInterpolationScheme<Type>::uncorrectedInterpolate"
 170                "(const GeometricField<Type, faPatchField, areaMesh>&, "
 171                "const tmp<edgeScalarField>&, "
 172                "const tmp<edgeScalarField>&) : "
 173                "interpolating areaTypeField from cells to faces "
 174                "without explicit correction"
 175             << endl;
 176     }
 177
 178     const edgeScalarField& lambdas = tlambdas();
 179     const edgeScalarField& ys = tys();
 180
 181     const Field<Type>& vfi = vf.internalField();
 182     const scalarField& lambda = lambdas.internalField();
 183     const scalarField& y = ys.internalField();
 184
 185     const faMesh& mesh = vf.mesh();
 186     const unallocLabelList& P = mesh.owner();
 187     const unallocLabelList& N = mesh.neighbour();
 188
 189     tmp<GeometricField<Type, faePatchField, edgeMesh> > tsf
 190     (
 191         new GeometricField<Type, faePatchField, edgeMesh>
 192         (
 193             IOobject
 194             (
 195                 "interpolate("+vf.name()+')',
 196                 vf.instance(),
 197                 vf.db()
 198             ),
 199             mesh,
 200             vf.dimensions()
 201         )
 202     );
 203     GeometricField<Type, faePatchField, edgeMesh>& sf = tsf();
 204
 205     Field<Type>& sfi = sf.internalField();
 206
 207     for (register label fi=0; fi<P.size(); fi++)
 208     {
 209         // ZT, 22/Apr/2003
 210         const tensorField& curT = mesh.edgeTransformTensors()[fi];
 211
 212         const tensor& Te = curT[0];
 213         const tensor& TP = curT[1];
 214         const tensor& TN = curT[2];
 215
 216         sfi[fi] =
 217             transform
 218             (
 219                 Te.T(),
 220                 lambda[fi]*transform(TP, vfi[P[fi]])
 221               + y[fi]*transform(TN, vfi[N[fi]])
 222             );
 223     }
 224
 225
 226     // Interpolate across coupled patches using given lambdas and ys
 227
 228     forAll (lambdas.boundaryField(), pi)
 229     {
 230         const faePatchScalarField& pLambda = lambdas.boundaryField()[pi];
 231         const faePatchScalarField& pY = ys.boundaryField()[pi];
 232
 233         if (vf.boundaryField()[pi].coupled())
 234         {
 235             label size = vf.boundaryField()[pi].patch().size();
 236             label start = vf.boundaryField()[pi].patch().start();
 237
 238             Field<Type> pOwnVf = vf.boundaryField()[pi].patchInternalField();
 239             Field<Type> pNgbVf = vf.boundaryField()[pi].patchNeighbourField();
 240
 241             Field<Type>& pSf = sf.boundaryField()[pi];
 242
 243             for (label i=0; i<size; i++)
 244             {
 245                 const tensorField& curT =
 246                     mesh.edgeTransformTensors()[start + i];
 247
 248                 const tensor& Te = curT[0];
 249                 const tensor& TP = curT[1];
 250                 const tensor& TN = curT[2];
 251
 252                 pSf[i] =
 253                     transform
 254                     (
 255                         Te.T(),
 256                         pLambda[i]*transform(TP, pOwnVf[i])
 257                       + pY[i]*transform(TN, pNgbVf[i])
 258                     );
 259             }
 260
 261 //             sf.boundaryField()[pi] =
 262 //                 pLambda*vf.boundaryField()[pi].patchInternalField()
 263 //               + pY*vf.boundaryField()[pi].patchNeighbourField();
 264         }
 265         else
 266         {
 267             sf.boundaryField()[pi] = vf.boundaryField()[pi];
 268         }
 269     }
 270
 271     tlambdas.clear();
 272     tys.clear();
 273
 274     return tsf;
 275 }
 276
 277
 278 //- Return the face-interpolate of the given cell field
 279 //  with the given weigting factors
 280 template<class Type>
 281 tmp<GeometricField<Type, faePatchField, edgeMesh> >
 282 edgeInterpolationScheme<Type>::interpolate
 283 (
 284     const GeometricField<Type, faPatchField, areaMesh>& vf,
 285     const tmp<edgeScalarField>& tlambdas
 286 )
 287 {
 288     if (edgeInterpolation::debug)
 289     {
 290         Info<< "edgeInterpolationScheme<Type>::interpolate"
 291                "(const GeometricField<Type, faPatchField, areaMesh>&, "
 292                "const tmp<edgeScalarField>&) : "
 293                "interpolating areaTypeField from cells to faces "
 294                "without explicit correction"
 295             << endl;
 296     }
 297
 298     const edgeScalarField& lambdas = tlambdas();
 299
 300     const Field<Type>& vfi = vf.internalField();
 301     const scalarField& lambda = lambdas.internalField();
 302
 303     const faMesh& mesh = vf.mesh();
 304     const unallocLabelList& P = mesh.owner();
 305     const unallocLabelList& N = mesh.neighbour();
 306
 307     tmp<GeometricField<Type, faePatchField, edgeMesh> > tsf
 308     (
 309         new GeometricField<Type, faePatchField, edgeMesh>
 310         (
 311             IOobject
 312             (
 313                 "interpolate("+vf.name()+')',
 314                 vf.instance(),
 315                 vf.db()
 316             ),
 317             mesh,
 318             vf.dimensions()
 319         )
 320     );
 321     GeometricField<Type, faePatchField, edgeMesh>& sf = tsf();
 322
 323     Field<Type>& sfi = sf.internalField();
 324
 325     for (register label eI = 0; eI < P.size(); eI++)
 326     {
 327         // ZT, 22/Apr/2003
 328         const tensorField& curT = mesh.edgeTransformTensors()[eI];
 329
 330         const tensor& Te = curT[0];
 331         const tensor& TP = curT[1];
 332         const tensor& TN = curT[2];
 333
 334         sfi[eI] =
 335             transform
 336             (
 337                 Te.T(),
 338                 lambda[eI]*transform(TP, vfi[P[eI]])
 339               + (1 - lambda[eI])*transform(TN, vfi[N[eI]])
 340             );
 341     }
 342
 343
 344     // Interpolate across coupled patches using given lambdas
 345
 346     forAll (lambdas.boundaryField(), pi)
 347     {
 348         const faePatchScalarField& pLambda = lambdas.boundaryField()[pi];
 349
 350         if (vf.boundaryField()[pi].coupled())
 351         {
 352             label size = vf.boundaryField()[pi].patch().size();
 353             label start = vf.boundaryField()[pi].patch().start();
 354
 355             Field<Type> pOwnVf = vf.boundaryField()[pi].patchInternalField();
 356             Field<Type> pNgbVf = vf.boundaryField()[pi].patchNeighbourField();
 357
 358             Field<Type>& pSf = sf.boundaryField()[pi];
 359
 360             for (label i=0; i<size; i++)
 361             {
 362                 const tensorField& curT =
 363                     mesh.edgeTransformTensors()[start + i];
 364
 365                 const tensor& Te = curT[0];
 366                 const tensor& TP = curT[1];
 367                 const tensor& TN = curT[2];
 368
 369                 pSf[i] =
 370                     transform
 371                     (
 372                         Te.T(),
 373                         pLambda[i]*transform(TP, pOwnVf[i])
 374                       + (1 - pLambda[i])*transform(TN, pNgbVf[i])
 375                     );
 376             }
 377
 378 //             tsf().boundaryField()[pi] =
 379 //                 pLambda*vf.boundaryField()[pi].patchInternalField()
 380 //              + (1 - pLambda)*vf.boundaryField()[pi].patchNeighbourField();
 381         }
 382         else
 383         {
 384             sf.boundaryField()[pi] = vf.boundaryField()[pi];
 385         }
 386     }
 387
 388     tlambdas.clear();
 389
 390     return tsf;
 391 }
 392
 393
 394 //- Return the euclidian edge-interpolate of the given area field
 395 //  with the given weigting factors
 396 template<class Type>
 397 tmp<GeometricField<Type, faePatchField, edgeMesh> >
 398 edgeInterpolationScheme<Type>::euclidianInterpolate
 399 (
 400     const GeometricField<Type, faPatchField, areaMesh>& vf,
 401     const tmp<edgeScalarField>& tlambdas
 402 )
 403 {
 404     if (edgeInterpolation::debug)
 405     {
 406         Info<< "edgeInterpolationScheme<Type>::euclidianInterpolate"
 407                "(const GeometricField<Type, faPatchField, areaMesh>&, "
 408                "const tmp<edgeScalarField>&) : "
 409                "interpolating areaTypeField from cells to faces "
 410                "without explicit correction"
 411             << endl;
 412     }
 413
 414     const edgeScalarField& lambdas = tlambdas();
 415
 416     const Field<Type>& vfi = vf.internalField();
 417     const scalarField& lambda = lambdas.internalField();
 418
 419     const faMesh& mesh = vf.mesh();
 420     const unallocLabelList& P = mesh.owner();
 421     const unallocLabelList& N = mesh.neighbour();
 422
 423     tmp<GeometricField<Type, faePatchField, edgeMesh> > tsf
 424     (
 425         new GeometricField<Type, faePatchField, edgeMesh>
 426         (
 427             IOobject
 428             (
 429                 "interpolate("+vf.name()+')',
 430                 vf.instance(),
 431                 vf.db()
 432             ),
 433             mesh,
 434             vf.dimensions()
 435         )
 436     );
 437     GeometricField<Type, faePatchField, edgeMesh>& sf = tsf();
 438
 439     Field<Type>& sfi = sf.internalField();
 440
 441     for (register label eI = 0; eI < P.size(); eI++)
 442     {
 443         sfi[eI] = lambda[eI]*vfi[P[eI]] + (1 - lambda[eI])*vfi[N[eI]];
 444     }
 445
 446
 447     // Interpolate across coupled patches using given lambdas
 448
 449     forAll (lambdas.boundaryField(), pi)
 450     {
 451         const faePatchScalarField& pLambda = lambdas.boundaryField()[pi];
 452
 453         if (vf.boundaryField()[pi].coupled())
 454         {
 455             tsf().boundaryField()[pi] =
 456                 pLambda*vf.boundaryField()[pi].patchInternalField()
 457              + (1.0 - pLambda)*vf.boundaryField()[pi].patchNeighbourField();
 458         }
 459         else
 460         {
 461             sf.boundaryField()[pi] = vf.boundaryField()[pi];
 462         }
 463     }
 464
 465     tlambdas.clear();
 466
 467     return tsf;
 468 }
 469
 470
 471 //- Return the face-interpolate of the given cell field
 472 //  with explicit correction
 473 template<class Type>
 474 tmp<GeometricField<Type, faePatchField, edgeMesh> >
 475 edgeInterpolationScheme<Type>::interpolate
 476 (
 477     const GeometricField<Type, faPatchField, areaMesh>& vf
 478 ) const
 479 {
 480     if (edgeInterpolation::debug)
 481     {
 482         Info<< "edgeInterpolationScheme<Type>::interpolate"
 483                "(const GeometricField<Type, faPatchField, areaMesh>&) : "
 484             << "interpolating areaTypeField from cells to faces"
 485             << endl;
 486     }
 487
 488     tmp<GeometricField<Type, faePatchField, edgeMesh> > tsf
 489         = interpolate(vf, weights(vf));
 490
 491     if (corrected())
 492     {
 493         tsf() += correction(vf);
 494     }
 495
 496     return tsf;
 497 }
 498
 499 //- Return the euclidian edge-interpolate of the given area field
 500 //  without explicit correction
 501 template<class Type>
 502 tmp<GeometricField<Type, faePatchField, edgeMesh> >
 503 edgeInterpolationScheme<Type>::euclidianInterpolate
 504 (
 505     const GeometricField<Type, faPatchField, areaMesh>& vf
 506 ) const
 507 {
 508     if (edgeInterpolation::debug)
 509     {
 510         Info<< "edgeInterpolationScheme<Type>::interpolate"
 511                "(const GeometricField<Type, faPatchField, areaMesh>&) : "
 512             << "interpolating areaTypeField from cells to faces"
 513             << endl;
 514     }
 515
 516     tmp<GeometricField<Type, faePatchField, edgeMesh> > tsf
 517         = euclidianInterpolate(vf, weights(vf));
 518
 519     return tsf;
 520 }
 521
 522 //- Return the face-interpolate of the given cell field
 523 //  with explicit correction
 524 template<class Type>
 525 tmp<GeometricField<Type, faePatchField, edgeMesh> >
 526 edgeInterpolationScheme<Type>::interpolate
 527 (
 528     const tmp<GeometricField<Type, faPatchField, areaMesh> >& tvf
 529 ) const
 530 {
 531     tmp<GeometricField<Type, faePatchField, edgeMesh> > tinterpVf
 532         = interpolate(tvf());
 533     tvf.clear();
 534     return tinterpVf;
 535 }
 536
 537
 538 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 539
 540 } // End namespace Foam
 541
 542 // ************************************************************************* //