src/dynamicMesh/dynamicTopoFvMesh/convexSetAlgorithm/convexSetAlgorithm.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     convexSetAlgorithm
  26
  27 Description
  28     Base class for convexSetAlgorithms
  29
  30 Author
  31     Sandeep Menon
  32     University of Massachusetts Amherst
  33     All rights reserved
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "objectRegistry.H"
  38 #include "foamTime.H"
  39 #include "IOMap.H"
  40 #include "meshOps.H"
  41 #include "polyMesh.H"
  42 #include "objectMap.H"
  43 #include "edgeIOList.H"
  44 #include "cellIOList.H"
  45
  46 #include "convexSetAlgorithm.H"
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 namespace Foam
  51 {
  52
  53 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  54
  55 // Construct from components
  56 convexSetAlgorithm::convexSetAlgorithm
  57 (
  58     const polyMesh& mesh,
  59     const pointField& newPoints,
  60     const UList<edge>& newEdges,
  61     const UList<face>& newFaces,
  62     const UList<cell>& newCells,
  63     const UList<label>& newOwner,
  64     const UList<label>& newNeighbour
  65 )
  66 :
  67     nOldPoints_(mesh.nPoints()),
  68     mesh_(mesh),
  69     newPoints_(newPoints),
  70     newEdges_(newEdges),
  71     newFaces_(newFaces),
  72     newCells_(newCells),
  73     newOwner_(newOwner),
  74     newNeighbour_(newNeighbour)
  75 {}
  76
  77
  78 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  79
  80 // Obtain map weighting factors
  81 void convexSetAlgorithm::computeWeights
  82 (
  83     const label index,
  84     const label offset,
  85     const labelList& mapCandidates,
  86     const labelListList& oldNeighbourList,
  87     labelList& parents,
  88     scalarField& weights,
  89     vectorField& centres,
  90     bool output
  91 )
  92 {
  93     if (parents.size() || weights.size() || centres.size())
  94     {
  95         FatalErrorIn
  96         (
  97             "\n\n"
  98             "void convexSetAlgorithm::computeWeights\n"
  99             "(\n"
 100             "    const label index,\n"
 101             "    const label offset,\n"
 102             "    const labelList& mapCandidates,\n"
 103             "    const labelListList& oldNeighbourList,\n"
 104             "    labelList& parents,\n"
 105             "    scalarField& weights,\n"
 106             "    vectorField& centres,\n"
 107             "    bool output\n"
 108             ")\n"
 109         )
 110             << " Addressing has already been calculated." << nl
 111             << " Index: " << index << nl
 112             << " Offset: " << offset << nl
 113             << " Type: " << (dimension() == 2 ? "Face" : "Cell") << nl
 114             << " mapCandidates: " << mapCandidates << nl
 115             << " Parents: " << parents << nl
 116             << " Weights: " << weights << nl
 117             << " Centres: " << centres << nl
 118             << abort(FatalError);
 119     }
 120
 121     // Do nothing for empty lists
 122     if (mapCandidates.empty() || oldNeighbourList.empty())
 123     {
 124         return;
 125     }
 126
 127     bool changed;
 128     label nAttempts = 0, nIntersects = 0;
 129
 130     // Calculate the algorithm normFactor
 131     computeNormFactor(index);
 132
 133     // Maintain a check-list
 134     labelHashSet checked, skipped;
 135
 136     // Loop and add intersections until nothing changes
 137     do
 138     {
 139         // Reset flag
 140         changed = false;
 141
 142         // Fetch the set of candidates
 143         labelList checkList;
 144
 145         if (nAttempts == 0)
 146         {
 147             checkList = mapCandidates;
 148         }
 149         else
 150         {
 151             checkList = checked.toc();
 152         }
 153
 154         forAll(checkList, indexI)
 155         {
 156             labelList checkEntities;
 157
 158             if (nAttempts == 0)
 159             {
 160                 checkEntities = labelList(1, checkList[indexI] - offset);
 161             }
 162             else
 163             {
 164                 checkEntities = oldNeighbourList[checkList[indexI]];
 165             }
 166
 167             forAll(checkEntities, entityI)
 168             {
 169                 label checkEntity = checkEntities[entityI];
 170
 171                 // Skip if this is already
 172                 // on the checked / skipped list
 173                 if
 174                 (
 175                     (checked.found(checkEntity)) ||
 176                     (skipped.found(checkEntity))
 177                 )
 178                 {
 179                     continue;
 180                 }
 181
 182                 bool intersect =
 183                 (
 184                     computeIntersection
 185                     (
 186                         index,
 187                         checkEntity + offset,
 188                         offset,
 189                         output
 190                     )
 191                 );
 192
 193                 if (intersect)
 194                 {
 195                     nIntersects++;
 196
 197                     if (!checked.found(checkEntity))
 198                     {
 199                         checked.insert(checkEntity);
 200                     }
 201
 202                     changed = true;
 203                 }
 204                 else
 205                 {
 206                     // Add to the skipped list
 207                     if (!skipped.found(checkEntity))
 208                     {
 209                         skipped.insert(checkEntity);
 210                     }
 211                 }
 212             }
 213         }
 214
 215         if (nAttempts == 0 && !changed)
 216         {
 217             // Need to setup a rescue mechanism.
 218             labelHashSet rescue;
 219
 220             forAll(mapCandidates, cI)
 221             {
 222                 if (!rescue.found(mapCandidates[cI] - offset))
 223                 {
 224                     rescue.insert(mapCandidates[cI] - offset);
 225                 }
 226             }
 227
 228             for (label level = 0; level < 10; level++)
 229             {
 230                 labelList initList = rescue.toc();
 231
 232                 forAll(initList, fI)
 233                 {
 234                     const labelList& ff = oldNeighbourList[initList[fI]];
 235
 236                     forAll(ff, entityI)
 237                     {
 238                         if (!rescue.found(ff[entityI]))
 239                         {
 240                             rescue.insert(ff[entityI]);
 241                         }
 242                     }
 243                 }
 244             }
 245
 246             labelList finalList = rescue.toc();
 247
 248             forAll(finalList, entityI)
 249             {
 250                 label checkEntity = finalList[entityI];
 251
 252                 bool intersect =
 253                 (
 254                     computeIntersection
 255                     (
 256                         index,
 257                         checkEntity + offset,
 258                         offset,
 259                         output
 260                     )
 261                 );
 262
 263                 if (intersect)
 264                 {
 265                     nIntersects++;
 266
 267                     if (!checked.found(checkEntity))
 268                     {
 269                         checked.insert(checkEntity);
 270                     }
 271
 272                     changed = true;
 273                     break;
 274                 }
 275             }
 276
 277             if (!changed)
 278             {
 279                 // No point in continuing further...
 280                 break;
 281             }
 282         }
 283
 284         nAttempts++;
 285
 286         // Break out if we're taking too long
 287         if (nAttempts > 20)
 288         {
 289             break;
 290         }
 291
 292     } while (changed);
 293
 294     // Normalize weights
 295     normalize(false);
 296
 297     // Populate lists
 298     populateLists(parents, centres, weights);
 299 }
 300
 301
 302 // Output an entity as a VTK file
 303 void convexSetAlgorithm::writeVTK
 304 (
 305     const word& name,
 306     const label entity,
 307     const label primitiveType,
 308     const bool useOldConnectivity
 309 ) const
 310 {
 311     writeVTK
 312     (
 313         name,
 314         labelList(1, entity),
 315         primitiveType,
 316         useOldConnectivity
 317     );
 318 }
 319
 320
 321 // Output a list of entities as a VTK file
 322 void convexSetAlgorithm::writeVTK
 323 (
 324     const word& name,
 325     const labelList& cList,
 326     const label primitiveType,
 327     const bool useOldConnectivity
 328 ) const
 329 {
 330     if (useOldConnectivity)
 331     {
 332         const polyMesh& mesh = this->mesh_;
 333
 334         meshOps::writeVTK
 335         (
 336             mesh,
 337             name,
 338             cList,
 339             primitiveType,
 340             mesh.points(),
 341             mesh.edges(),
 342             mesh.faces(),
 343             mesh.cells(),
 344             mesh.faceOwner()
 345         );
 346     }
 347     else
 348     {
 349         meshOps::writeVTK
 350         (
 351             this->mesh_,
 352             name,
 353             cList,
 354             primitiveType,
 355             newPoints_,
 356             newEdges_,
 357             newFaces_,
 358             newCells_,
 359             newOwner_
 360         );
 361     }
 362 }
 363
 364
 365 bool convexSetAlgorithm::consistent(const scalar tolerance) const
 366 {
 367     if (weights_.size())
 368     {
 369         scalar normError = mag(1.0 - (sum(weights_)/normFactor_));
 370
 371         if (normError < tolerance)
 372         {
 373             return true;
 374         }
 375         else
 376         {
 377             return false;
 378         }
 379     }
 380
 381     return false;
 382 }
 383
 384
 385 // Return the normFactor
 386 scalar convexSetAlgorithm::normFactor() const
 387 {
 388     return normFactor_;
 389 }
 390
 391
 392 // Normalize stored weights
 393 void convexSetAlgorithm::normalize(bool normSum) const
 394 {
 395     if (normSum)
 396     {
 397         if (weights_.size())
 398         {
 399             weights_ /= sum(weights_);
 400         }
 401     }
 402     else
 403     {
 404         if (weights_.size())
 405         {
 406             weights_ /= normFactor_;
 407         }
 408     }
 409 }
 410
 411
 412 // Extract weights and centres to lists
 413 void convexSetAlgorithm::populateLists
 414 (
 415     labelList& parents,
 416     vectorField& centres,
 417     scalarField& weights
 418 ) const
 419 {
 420     // Clear inputs
 421     parents.clear();
 422     centres.clear();
 423     weights.clear();
 424
 425     if (weights_.size())
 426     {
 427         parents = parents_;
 428         centres = centres_;
 429         weights = weights_;
 430     }
 431 }
 432
 433
 434 // Write out connectivity information to disk
 435 bool convexSetAlgorithm::write() const
 436 {
 437     pointIOField
 438     (
 439         IOobject
 440         (
 441             "newPoints",
 442             mesh_.time().timeName(),
 443             "convexSetAlgorithm",
 444             mesh_,
 445             IOobject::NO_READ,
 446             IOobject::NO_WRITE,
 447             false
 448         ),
 449         newPoints_
 450     ).writeObject
 451     (
 452         IOstream::BINARY,
 453         IOstream::currentVersion,
 454         IOstream::COMPRESSED
 455     );
 456
 457     edgeIOList
 458     (
 459         IOobject
 460         (
 461             "newEdges",
 462             mesh_.time().timeName(),
 463             "convexSetAlgorithm",
 464             mesh_,
 465             IOobject::NO_READ,
 466             IOobject::NO_WRITE,
 467             false
 468         ),
 469         newEdges_
 470     ).writeObject
 471     (
 472         IOstream::BINARY,
 473         IOstream::currentVersion,
 474         IOstream::COMPRESSED
 475     );
 476
 477     faceIOList
 478     (
 479         IOobject
 480         (
 481             "newFaces",
 482             mesh_.time().timeName(),
 483             "convexSetAlgorithm",
 484             mesh_,
 485             IOobject::NO_READ,
 486             IOobject::NO_WRITE,
 487             false
 488         ),
 489         newFaces_
 490     ).writeObject
 491     (
 492         IOstream::BINARY,
 493         IOstream::currentVersion,
 494         IOstream::COMPRESSED
 495     );
 496
 497     cellIOList
 498     (
 499         IOobject
 500         (
 501             "newCells",
 502             mesh_.time().timeName(),
 503             "convexSetAlgorithm",
 504             mesh_,
 505             IOobject::NO_READ,
 506             IOobject::NO_WRITE,
 507             false
 508         ),
 509         newCells_
 510     ).writeObject
 511     (
 512         IOstream::BINARY,
 513         IOstream::currentVersion,
 514         IOstream::COMPRESSED
 515     );
 516
 517     labelIOList
 518     (
 519         IOobject
 520         (
 521             "newOwner",
 522             mesh_.time().timeName(),
 523             "convexSetAlgorithm",
 524             mesh_,
 525             IOobject::NO_READ,
 526             IOobject::NO_WRITE,
 527             false
 528         ),
 529         newOwner_
 530     ).writeObject
 531     (
 532         IOstream::BINARY,
 533         IOstream::currentVersion,
 534         IOstream::COMPRESSED
 535     );
 536
 537     labelIOList
 538     (
 539         IOobject
 540         (
 541             "newNeighbour",
 542             mesh_.time().timeName(),
 543             "convexSetAlgorithm",
 544             mesh_,
 545             IOobject::NO_READ,
 546             IOobject::NO_WRITE,
 547             false
 548         ),
 549         newNeighbour_
 550     ).writeObject
 551     (
 552         IOstream::BINARY,
 553         IOstream::currentVersion,
 554         IOstream::COMPRESSED
 555     );
 556
 557     return true;
 558 }
 559
 560
 561 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 562
 563 } // End namespace Foam
 564
 565 // ************************************************************************* //