src/dynamicMesh/dynamicTopoFvMesh/fieldMapping/topoCellMapper.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 Class
  25     topoCellMapper
  26
  27 Description
  28     Implementation of the topoCellMapper class
  29
  30 Author
  31     Sandeep Menon
  32     University of Massachusetts Amherst
  33     All rights reserved
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "topoMapper.H"
  38 #include "mapPolyMesh.H"
  39 #include "topoCellMapper.H"
  40 #include "demandDrivenData.H"
  41
  42 namespace Foam
  43 {
  44
  45 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  46
  47 //- Clear out local storage
  48 void topoCellMapper::clearOut()
  49 {
  50     deleteDemandDrivenData(directAddrPtr_);
  51     deleteDemandDrivenData(interpolationAddrPtr_);
  52     deleteDemandDrivenData(weightsPtr_);
  53     deleteDemandDrivenData(insertedCellLabelsPtr_);
  54     deleteDemandDrivenData(volumesPtr_);
  55     deleteDemandDrivenData(centresPtr_);
  56 }
  57
  58
  59 //- Calculate addressing for interpolative mapping
  60 void topoCellMapper::calcAddressing() const
  61 {
  62     if (directAddrPtr_ || interpolationAddrPtr_)
  63     {
  64         FatalErrorIn("void topoCellMapper::calcAddressing() const")
  65             << "Addressing already calculated."
  66             << abort(FatalError);
  67     }
  68
  69     // Allocate for inserted cell labels
  70     label nInsertedCells = 0;
  71
  72     insertedCellLabelsPtr_ = new labelList(mesh_.nCells(), -1);
  73     labelList& insertedCells = *insertedCellLabelsPtr_;
  74
  75     if (direct())
  76     {
  77         // Direct addressing, no weights
  78         directAddrPtr_ = new labelList(mpm_.cellMap());
  79     }
  80     else
  81     {
  82         // Interpolative addressing
  83         interpolationAddrPtr_ = new labelListList(mesh_.nCells());
  84         labelListList& addr = *interpolationAddrPtr_;
  85
  86         const List<objectMap>& cfc = mpm_.cellsFromCellsMap();
  87
  88         forAll (cfc, cfcI)
  89         {
  90             // Get addressing
  91             const labelList& mo = cfc[cfcI].masterObjects();
  92
  93             label cellI = cfc[cfcI].index();
  94
  95             if (addr[cellI].size() > 0)
  96             {
  97                 FatalErrorIn("void topoCellMapper::calcAddressing() const")
  98                     << "Master cell " << cellI
  99                     << " mapped from cells " << mo
 100                     << " is already destination for mapping."
 101                     << abort(FatalError);
 102             }
 103
 104             // Set master objects
 105             addr[cellI] = mo;
 106         }
 107
 108         // Do mapped cells. Note that this can already be set by cellsFromCells
 109         // so check if addressing size still zero.
 110         const labelList& cm = mpm_.cellMap();
 111
 112         forAll (cm, cellI)
 113         {
 114             // Mapped from a single cell
 115             if (cm[cellI] > -1 && addr[cellI].empty())
 116             {
 117                 addr[cellI] = labelList(1, cm[cellI]);
 118             }
 119
 120             // Check for inserted cells without any addressing
 121             if (cm[cellI] < 0 && addr[cellI].empty())
 122             {
 123                 insertedCells[nInsertedCells++] = cellI;
 124             }
 125         }
 126     }
 127
 128     // Shorten inserted cells to actual size
 129     insertedCells.setSize(nInsertedCells);
 130
 131     if (nInsertedCells)
 132     {
 133         FatalErrorIn("void topoCellMapper::calcAddressing() const")
 134             << " Found " << nInsertedCells << " which are"
 135             << " not mapped from any parent cells." << nl
 136             << " List: " << nl
 137             << insertedCells
 138             << abort(FatalError);
 139     }
 140 }
 141
 142
 143 //- Calculate inverse-distance weights for interpolative mapping
 144 void topoCellMapper::calcInverseDistanceWeights() const
 145 {
 146     if (weightsPtr_)
 147     {
 148         FatalErrorIn
 149         (
 150             "void topoCellMapper::calcInverseDistanceWeights() const"
 151         )
 152             << "Weights already calculated."
 153             << abort(FatalError);
 154     }
 155
 156     // Fetch interpolative addressing
 157     const labelListList& addr = addressing();
 158
 159     // Allocate memory
 160     weightsPtr_ = new scalarListList(size());
 161     scalarListList& w = *weightsPtr_;
 162
 163     // Obtain cell-centre information from old/new meshes
 164     const vectorField& oldCentres = tMapper_.internalCentres();
 165     const vectorField& newCentres = mesh_.cellCentres();
 166
 167     forAll(addr, cellI)
 168     {
 169         const labelList& mo = addr[cellI];
 170
 171         // Do mapped cells
 172         if (mo.size() == 1)
 173         {
 174             w[cellI] = scalarList(1, 1.0);
 175         }
 176         else
 177         {
 178             // Map from masters, inverse-distance weights
 179             scalar totalWeight = 0.0;
 180             w[cellI] = scalarList(mo.size(), 0.0);
 181
 182             forAll (mo, oldCellI)
 183             {
 184                 w[cellI][oldCellI] =
 185                 (
 186                     1.0/stabilise
 187                     (
 188                         magSqr
 189                         (
 190                             newCentres[cellI]
 191                           - oldCentres[mo[oldCellI]]
 192                         ),
 193                         VSMALL
 194                     )
 195                 );
 196
 197                 totalWeight += w[cellI][oldCellI];
 198             }
 199
 200             // Normalize weights
 201             scalar normFactor = (1.0/totalWeight);
 202
 203             forAll (mo, oldCellI)
 204             {
 205                 w[cellI][oldCellI] *= normFactor;
 206             }
 207         }
 208     }
 209 }
 210
 211
 212 //- Calculate intersection weights for conservative mapping
 213 void topoCellMapper::calcIntersectionWeightsAndCentres() const
 214 {
 215     if (volumesPtr_ || centresPtr_)
 216     {
 217         FatalErrorIn
 218         (
 219             "void topoCellMapper::"
 220             "calcIntersectionWeightsAndCentres() const"
 221         )
 222             << "Weights already calculated."
 223             << abort(FatalError);
 224     }
 225
 226     // Fetch interpolative addressing
 227     const labelListList& addr = addressing();
 228
 229     // Allocate memory
 230     volumesPtr_ = new List<scalarField>(size(), scalarField(0));
 231     List<scalarField>& v = *volumesPtr_;
 232
 233     centresPtr_ = new List<vectorField>(size(), vectorField(0));
 234     List<vectorField>& x = *centresPtr_;
 235
 236     // Obtain stored cell-centres
 237     const vectorField& cellCentres = tMapper_.internalCentres();
 238
 239     // Fetch maps
 240     const List<objectMap>& cfc = mpm_.cellsFromCellsMap();
 241     const List<vectorField>& mapCellCentres = tMapper_.cellCentres();
 242     const List<scalarField>& mapCellWeights = tMapper_.cellWeights();
 243
 244     // Fill in maps first
 245     forAll(cfc, indexI)
 246     {
 247         x[cfc[indexI].index()] = mapCellCentres[indexI];
 248         v[cfc[indexI].index()] = mapCellWeights[indexI];
 249     }
 250
 251     // Now do mapped cells
 252     forAll(addr, cellI)
 253     {
 254         const labelList& mo = addr[cellI];
 255
 256         // Check if this is indeed a mapped cell
 257         if (mo.size() == 1 && x[cellI].empty() && v[cellI].empty())
 258         {
 259             x[cellI] = vectorField(1, cellCentres[mo[0]]);
 260             v[cellI] = scalarField(1, 1.0);
 261         }
 262     }
 263 }
 264
 265
 266 const List<scalarField>& topoCellMapper::intersectionWeights() const
 267 {
 268     if (direct())
 269     {
 270         FatalErrorIn
 271         (
 272             "const List<scalarField>& "
 273             "topoCellMapper::intersectionWeights() const"
 274         )   << "Requested interpolative weights for a direct mapper."
 275             << abort(FatalError);
 276     }
 277
 278     if (!volumesPtr_)
 279     {
 280         calcIntersectionWeightsAndCentres();
 281     }
 282
 283     return *volumesPtr_;
 284 }
 285
 286
 287 const List<vectorField>& topoCellMapper::intersectionCentres() const
 288 {
 289     if (direct())
 290     {
 291         FatalErrorIn
 292         (
 293             "const List<vectorField>& "
 294             "topoCellMapper::intersectionCentres() const"
 295         )   << "Requested interpolative weights for a direct mapper."
 296             << abort(FatalError);
 297     }
 298
 299     if (!centresPtr_)
 300     {
 301         calcIntersectionWeightsAndCentres();
 302     }
 303
 304     return *centresPtr_;
 305 }
 306
 307
 308 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 309
 310 // Construct from components
 311 topoCellMapper::topoCellMapper
 312 (
 313     const mapPolyMesh& mpm,
 314     const topoMapper& mapper
 315 )
 316 :
 317     mesh_(mpm.mesh()),
 318     mpm_(mpm),
 319     tMapper_(mapper),
 320     direct_(false),
 321     sizeBeforeMapping_(mpm.nOldCells()),
 322     directAddrPtr_(NULL),
 323     interpolationAddrPtr_(NULL),
 324     weightsPtr_(NULL),
 325     insertedCellLabelsPtr_(NULL),
 326     volumesPtr_(NULL),
 327     centresPtr_(NULL)
 328 {
 329     // Fetch offset sizes from topoMapper
 330     const labelList& sizes = tMapper_.cellSizes();
 331
 332     // Add offset sizes
 333     if (sizes.size())
 334     {
 335         forAll(sizes, pI)
 336         {
 337             sizeBeforeMapping_ += sizes[pI];
 338         }
 339     }
 340
 341     // Check for possibility of direct mapping
 342     if
 343     (
 344         (min(mpm_.cellMap()) > -1)
 345      && mpm_.cellsFromPointsMap().empty()
 346      && mpm_.cellsFromEdgesMap().empty()
 347      && mpm_.cellsFromFacesMap().empty()
 348      && mpm_.cellsFromCellsMap().empty()
 349     )
 350     {
 351         direct_ = true;
 352     }
 353     else
 354     {
 355         direct_ = false;
 356     }
 357 }
 358
 359 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 360
 361 topoCellMapper::~topoCellMapper()
 362 {
 363     clearOut();
 364 }
 365
 366 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 367
 368 //- Return size
 369 label topoCellMapper::size() const
 370 {
 371     return mesh_.nCells();
 372 }
 373
 374
 375 //- Return size before mapping
 376 label topoCellMapper::sizeBeforeMapping() const
 377 {
 378     return sizeBeforeMapping_;
 379 }
 380
 381
 382 //- Is the mapping direct
 383 bool topoCellMapper::direct() const
 384 {
 385     return direct_;
 386 }
 387
 388
 389 //- Return direct addressing
 390 const unallocLabelList& topoCellMapper::directAddressing() const
 391 {
 392     if (!direct())
 393     {
 394         FatalErrorIn
 395         (
 396             "const unallocLabelList& "
 397             "topoCellMapper::directAddressing() const"
 398         )   << "Requested direct addressing for an interpolative mapper."
 399             << abort(FatalError);
 400     }
 401
 402     if (!directAddrPtr_)
 403     {
 404         calcAddressing();
 405     }
 406
 407     return *directAddrPtr_;
 408 }
 409
 410
 411 //- Return interpolation addressing
 412 const labelListList& topoCellMapper::addressing() const
 413 {
 414     if (direct())
 415     {
 416         FatalErrorIn
 417         (
 418             "const labelListList& "
 419             "topoCellMapper::addressing() const"
 420         )   << "Requested interpolative addressing for a direct mapper."
 421             << abort(FatalError);
 422     }
 423
 424     if (!interpolationAddrPtr_)
 425     {
 426         calcAddressing();
 427     }
 428
 429     return *interpolationAddrPtr_;
 430 }
 431
 432
 433 //- Return weights
 434 const scalarListList& topoCellMapper::weights() const
 435 {
 436     if (direct())
 437     {
 438         FatalErrorIn
 439         (
 440             "const scalarListList& "
 441             "topoCellMapper::weights() const"
 442         )   << "Requested interpolative weights for a direct mapper."
 443             << abort(FatalError);
 444     }
 445
 446     if (!weightsPtr_)
 447     {
 448         calcInverseDistanceWeights();
 449     }
 450
 451     return *weightsPtr_;
 452 }
 453
 454
 455 //- Are there any inserted cells
 456 bool topoCellMapper::insertedObjects() const
 457 {
 458     return insertedObjectLabels().size();
 459 }
 460
 461
 462 //- Return list of inserted cells
 463 const labelList& topoCellMapper::insertedObjectLabels() const
 464 {
 465     if (!insertedCellLabelsPtr_)
 466     {
 467         calcAddressing();
 468     }
 469
 470     return *insertedCellLabelsPtr_;
 471 }
 472
 473
 474 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 475
 476 void topoCellMapper::operator=(const topoCellMapper& rhs)
 477 {
 478     // Check for assignment to self
 479     if (this == &rhs)
 480     {
 481         FatalErrorIn("void topoCellMapper::operator=")
 482             << "Attempted assignment to self"
 483             << abort(FatalError);
 484     }
 485 }
 486
 487 } // End namespace Foam
 488
 489 // ************************************************************************* //