src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-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 "mapPolyMesh.H"
  27 #include "PstreamCombineReduceOps.H"
  28
  29 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  30
  31 //- combineReduce operator for lists. Used for counting.
  32 class listEq
  33 {
  34
  35 public:
  36
  37     template<class T>
  38     void operator()(T& x, const T& y) const
  39     {
  40         forAll(y, i)
  41         {
  42             if (y[i].size())
  43             {
  44                 x[i] = y[i];
  45             }
  46         }
  47     }
  48 };
  49
  50
  51 template <class Container, class T>
  52 void Foam::fvMeshDistribute::exchange
  53 (
  54     const List<Container >& sendBufs,
  55     List<Container >& recvBufs,
  56     labelListList& sizes
  57 )
  58 {
  59     if (Pstream::parRun())
  60     {
  61         if (!contiguous<T>())
  62         {
  63             FatalErrorIn
  64             (
  65                 "Pstream::exchange(..)"
  66             )   << "Continuous data only." << Foam::abort(FatalError);
  67         }
  68
  69         if (sendBufs.size() != Pstream::nProcs())
  70         {
  71             FatalErrorIn
  72             (
  73                 "Pstream::exchange(..)"
  74             )   << "Size of list:" << sendBufs.size()
  75                 << " does not equal the number of processors:"
  76                 << Pstream::nProcs()
  77                 << Foam::abort(FatalError);
  78         }
  79
  80         sizes.setSize(Pstream::nProcs());
  81         labelList& nsTransPs = sizes[Pstream::myProcNo()];
  82         nsTransPs.setSize(Pstream::nProcs());
  83
  84         forAll(sendBufs, procI)
  85         {
  86             nsTransPs[procI] = sendBufs[procI].size();
  87         }
  88
  89         Foam::combineReduce(sizes, listEq());
  90
  91
  92         // Set up receives
  93         // ~~~~~~~~~~~~~~~
  94
  95         recvBufs.setSize(sendBufs.size());
  96         forAll(sizes, procI)
  97         {
  98             label nRecv = sizes[procI][Pstream::myProcNo()];
  99
 100             if (procI != Pstream::myProcNo() && nRecv > 0)
 101             {
 102                 recvBufs[procI].setSize(nRecv);
 103                 IPstream::read
 104                 (
 105                     Pstream::nonBlocking,
 106                     procI,
 107                     reinterpret_cast<char*>(recvBufs[procI].begin()),
 108                     nRecv*sizeof(T)
 109                 );
 110             }
 111         }
 112
 113
 114         // Set up sends
 115         // ~~~~~~~~~~~~
 116
 117         forAll(sendBufs, procI)
 118         {
 119             if (procI != Pstream::myProcNo() && sendBufs[procI].size() > 0)
 120             {
 121                 if
 122                 (
 123                    !OPstream::write
 124                     (
 125                         Pstream::nonBlocking,
 126                         procI,
 127                         reinterpret_cast<const char*>(sendBufs[procI].begin()),
 128                         sendBufs[procI].size()*sizeof(T)
 129                     )
 130                 )
 131                 {
 132                     FatalErrorIn("Pstream::exchange(..)")
 133                         << "Cannot send outgoing message. "
 134                         << "to:" << procI << " nBytes:"
 135                         << label(sendBufs[procI].size()*sizeof(T))
 136                         << Foam::abort(FatalError);
 137                 }
 138             }
 139         }
 140
 141
 142         // Wait for all to finish
 143         // ~~~~~~~~~~~~~~~~~~~~~~
 144
 145         IPstream::waitRequests();
 146         OPstream::waitRequests();
 147     }
 148
 149     // Do myself
 150     recvBufs[Pstream::myProcNo()] = sendBufs[Pstream::myProcNo()];
 151 }
 152
 153
 154 template<class GeoField>
 155 void Foam::fvMeshDistribute::printFieldInfo(const fvMesh& mesh)
 156 {
 157     HashTable<const GeoField*> flds
 158     (
 159         mesh.objectRegistry::lookupClass<GeoField>()
 160     );
 161
 162     for
 163     (
 164         typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
 165         iter != flds.end();
 166         ++iter
 167     )
 168     {
 169         const GeoField& fld = *iter();
 170
 171         Pout<< "Field:" << iter.key() << " internalsize:" << fld.size()
 172             //<< " value:" << fld
 173             << endl;
 174
 175         forAll(fld.boundaryField(), patchI)
 176         {
 177             Pout<< "    " << patchI
 178                 << ' ' << fld.boundaryField()[patchI].patch().name()
 179                 << ' ' << fld.boundaryField()[patchI].type()
 180                 << ' ' << fld.boundaryField()[patchI].size()
 181                 << endl;
 182         }
 183     }
 184 }
 185
 186
 187 template<class GeoField>
 188 void Foam::fvMeshDistribute::addPatchFields(const word& patchFieldType)
 189 {
 190     HashTable<const GeoField*> flds
 191     (
 192         mesh_.objectRegistry::lookupClass<GeoField>()
 193     );
 194
 195     for
 196     (
 197         typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
 198         iter != flds.end();
 199         ++iter
 200     )
 201     {
 202         const GeoField& fld = *iter();
 203
 204         typename GeoField::GeometricBoundaryField& bfld =
 205             const_cast<typename GeoField::GeometricBoundaryField&>
 206             (
 207                 fld.boundaryField()
 208             );
 209
 210         label sz = bfld.size();
 211         bfld.setSize(sz + 1);
 212         bfld.set
 213         (
 214             sz,
 215             GeoField::PatchFieldType::New
 216             (
 217                 patchFieldType,
 218                 mesh_.boundary()[sz],
 219                 fld.dimensionedInternalField()
 220             )
 221         );
 222     }
 223 }
 224
 225
 226 // Delete trailing patch fields
 227 template<class GeoField>
 228 void Foam::fvMeshDistribute::deleteTrailingPatchFields()
 229 {
 230     HashTable<const GeoField*> flds
 231     (
 232         mesh_.objectRegistry::lookupClass<GeoField>()
 233     );
 234
 235     for
 236     (
 237         typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
 238         iter != flds.end();
 239         ++iter
 240     )
 241     {
 242         const GeoField& fld = *iter();
 243
 244         typename GeoField::GeometricBoundaryField& bfld =
 245             const_cast<typename GeoField::GeometricBoundaryField&>
 246             (
 247                 fld.boundaryField()
 248             );
 249
 250         // Shrink patchFields
 251         bfld.setSize(bfld.size() - 1);
 252     }
 253 }
 254
 255
 256 // Save whole boundary field
 257 template <class T, class Mesh>
 258 void Foam::fvMeshDistribute::saveBoundaryFields
 259 (
 260     PtrList<FieldField<fvsPatchField, T> >& bflds
 261 ) const
 262 {
 263     typedef GeometricField<T, fvsPatchField, Mesh> fldType;
 264
 265     HashTable<const fldType*> flds
 266     (
 267         mesh_.objectRegistry::lookupClass<fldType>()
 268     );
 269
 270     bflds.setSize(flds.size());
 271
 272     label i = 0;
 273
 274     for
 275     (
 276         typename HashTable<const fldType*>::const_iterator iter = flds.begin();
 277         iter != flds.end();
 278         ++iter
 279     )
 280     {
 281         const fldType& fld = *iter();
 282
 283         bflds.set(i, fld.boundaryField().clone().ptr());
 284
 285         i++;
 286     }
 287 }
 288
 289
 290 // Map boundary field
 291 template <class T, class Mesh>
 292 void Foam::fvMeshDistribute::mapBoundaryFields
 293 (
 294     const mapPolyMesh& map,
 295     const PtrList<FieldField<fvsPatchField, T> >& oldBflds
 296 )
 297 {
 298     const labelList& oldPatchStarts = map.oldPatchStarts();
 299     const labelList& faceMap = map.faceMap();
 300
 301     typedef GeometricField<T, fvsPatchField, Mesh> fldType;
 302
 303     HashTable<const fldType*> flds
 304     (
 305         mesh_.objectRegistry::lookupClass<fldType>()
 306     );
 307
 308     if (flds.size() != oldBflds.size())
 309     {
 310         FatalErrorIn("fvMeshDistribute::mapBoundaryFields(..)") << "problem"
 311             << abort(FatalError);
 312     }
 313
 314     label fieldI = 0;
 315
 316     for
 317     (
 318         typename HashTable<const fldType*>::const_iterator iter = flds.begin();
 319         iter != flds.end();
 320         ++iter
 321     )
 322     {
 323         const fldType& fld = *iter();
 324         typename fldType::GeometricBoundaryField& bfld =
 325             const_cast<typename fldType::GeometricBoundaryField&>
 326             (
 327                 fld.boundaryField()
 328             );
 329
 330
 331         const FieldField<fvsPatchField, T>& oldBfld = oldBflds[fieldI++];
 332
 333         // Pull from old boundary field into bfld.
 334
 335         forAll(bfld, patchI)
 336         {
 337             fvsPatchField<T>& patchFld = bfld[patchI];
 338             label faceI = patchFld.patch().patch().start();
 339
 340             forAll(patchFld, i)
 341             {
 342                 label oldFaceI = faceMap[faceI++];
 343
 344                 // Find patch and local patch face oldFaceI was in.
 345                 forAll(oldPatchStarts, oldPatchI)
 346                 {
 347                     label oldLocalI = oldFaceI - oldPatchStarts[oldPatchI];
 348
 349                     if (oldLocalI >= 0 && oldLocalI < oldBfld[oldPatchI].size())
 350                     {
 351                         patchFld[i] = oldBfld[oldPatchI][oldLocalI];
 352                     }
 353                 }
 354             }
 355         }
 356     }
 357 }
 358
 359
 360 // Init patch fields of certain type
 361 template<class GeoField, class PatchFieldType>
 362 void Foam::fvMeshDistribute::initPatchFields
 363 (
 364     const typename GeoField::value_type& initVal
 365 )
 366 {
 367     HashTable<const GeoField*> flds
 368     (
 369         mesh_.objectRegistry::lookupClass<GeoField>()
 370     );
 371
 372     for
 373     (
 374         typename HashTable<const GeoField*>::const_iterator iter = flds.begin();
 375         iter != flds.end();
 376         ++iter
 377     )
 378     {
 379         const GeoField& fld = *iter();
 380
 381         typename GeoField::GeometricBoundaryField& bfld =
 382             const_cast<typename GeoField::GeometricBoundaryField&>
 383             (
 384                 fld.boundaryField()
 385             );
 386
 387         forAll(bfld, patchI)
 388         {
 389             if (isA<PatchFieldType>(bfld[patchI]))
 390             {
 391                 bfld[patchI] == initVal;
 392             }
 393         }
 394     }
 395 }
 396
 397
 398 // Send fields. Note order supplied so we can receive in exactly the same order.
 399 // Note that field gets written as entry in dictionary so we
 400 // can construct from subdictionary.
 401 // (since otherwise the reading as-a-dictionary mixes up entries from
 402 // consecutive fields)
 403 // The dictionary constructed is:
 404 //  volScalarField
 405 //  {
 406 //      p {internalField ..; boundaryField ..;}
 407 //      k {internalField ..; boundaryField ..;}
 408 //  }
 409 //  volVectorField
 410 //  {
 411 //      U {internalField ...  }
 412 //  }
 413
 414 // volVectorField {U {internalField ..; boundaryField ..;}}
 415 //
 416 template<class GeoField>
 417 void Foam::fvMeshDistribute::sendFields
 418 (
 419     const label domain,
 420     const wordList& fieldNames,
 421     const fvMeshSubset& subsetter,
 422     OSstream& toNbr
 423 )
 424 {
 425     toNbr << GeoField::typeName << token::NL << token::BEGIN_BLOCK << token::NL;
 426     forAll(fieldNames, i)
 427     {
 428         if (debug)
 429         {
 430             Pout<< "Subsetting field " << fieldNames[i]
 431                 << " for domain:" << domain << endl;
 432         }
 433
 434         // Send all fieldNames. This has to be exactly the same set as is
 435         // being received!
 436         const GeoField& fld =
 437             subsetter.baseMesh().lookupObject<GeoField>(fieldNames[i]);
 438
 439         tmp<GeoField> tsubfld = subsetter.interpolate(fld);
 440
 441         toNbr
 442             << fieldNames[i] << token::NL << token::BEGIN_BLOCK
 443             << tsubfld
 444             << token::NL << token::END_BLOCK << token::NL;
 445     }
 446     toNbr << token::END_BLOCK << token::NL;
 447 }
 448
 449
 450 // Opposite of sendFields
 451 template<class GeoField>
 452 void Foam::fvMeshDistribute::receiveFields
 453 (
 454     const label domain,
 455     const wordList& fieldNames,
 456     fvMesh& mesh,
 457     PtrList<GeoField>& fields,
 458     const dictionary& fieldDicts
 459 )
 460 {
 461     if (debug)
 462     {
 463         Pout<< "Receiving fields " << fieldNames
 464             << " from domain:" << domain << endl;
 465     }
 466
 467     fields.setSize(fieldNames.size());
 468
 469     forAll(fieldNames, i)
 470     {
 471         if (debug)
 472         {
 473             Pout<< "Constructing field " << fieldNames[i]
 474                 << " from domain:" << domain << endl;
 475         }
 476
 477         fields.set
 478         (
 479             i,
 480             new GeoField
 481             (
 482                 IOobject
 483                 (
 484                     fieldNames[i],
 485                     mesh.time().timeName(),
 486                     mesh,
 487                     IOobject::NO_READ,
 488                     IOobject::AUTO_WRITE
 489                 ),
 490                 mesh,
 491                 fieldDicts.subDict(fieldNames[i])
 492             )
 493         );
 494     }
 495 }
 496
 497
 498 // ************************************************************************* //