src/dynamicMesh/dynamicTopoFvMesh/dynamicTopoFvMeshMapping.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     dynamicTopoFvMesh
  26
  27 Description
  28     Functions specific to conservative mapping
  29
  30 Author
  31     Sandeep Menon
  32     University of Massachusetts Amherst
  33     All rights reserved
  34
  35 \*---------------------------------------------------------------------------*/
  36
  37 #include "dynamicTopoFvMesh.H"
  38
  39 #include "foamTime.H"
  40 #include "meshOps.H"
  41 #include "IOmanip.H"
  42 #include "triFace.H"
  43 #include "objectMap.H"
  44 #include "faceSetAlgorithm.H"
  45 #include "cellSetAlgorithm.H"
  46
  47 namespace Foam
  48 {
  49
  50 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  51
  52 defineTemplateTypeNameAndDebugWithName(IOList<objectMap>, "objectMapIOList", 0);
  53
  54 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  55
  56 // Compute mapping weights for modified entities
  57 void dynamicTopoFvMesh::computeMapping
  58 (
  59     const scalar matchTol,
  60     const bool skipMapping,
  61     const bool mappingOutput,
  62     const label faceStart,
  63     const label faceSize,
  64     const label cellStart,
  65     const label cellSize
  66 )
  67 {
  68     // Convex-set algorithm for cells
  69     cellSetAlgorithm cellAlgorithm
  70     (
  71         (*this),
  72         oldPoints_,
  73         edges_,
  74         faces_,
  75         cells_,
  76         owner_,
  77         neighbour_
  78     );
  79
  80     label nInconsistencies = 0;
  81     scalar maxFaceError = 0.0, maxCellError = 0.0;
  82     DynamicList<scalar> cellErrors(10), faceErrors(10);
  83     DynamicList<objectMap> failedCells(10), failedFaces(10);
  84
  85     // Compute cell mapping
  86     for (label cellI = cellStart; cellI < (cellStart + cellSize); cellI++)
  87     {
  88         label cIndex = cellsFromCells_[cellI].index();
  89         labelList& masterObjects = cellsFromCells_[cellI].masterObjects();
  90
  91         if (skipMapping)
  92         {
  93             // Dummy map from cell[0]
  94             masterObjects = labelList(1, 0);
  95             cellWeights_[cellI].setSize(1, 1.0);
  96             cellCentres_[cellI].setSize(1, vector::zero);
  97         }
  98         else
  99         {
 100             // Obtain weighting factors for this cell.
 101             cellAlgorithm.computeWeights
 102             (
 103                 cIndex,
 104                 0,
 105                 cellParents_[cIndex],
 106                 polyMesh::cellCells(),
 107                 masterObjects,
 108                 cellWeights_[cellI],
 109                 cellCentres_[cellI]
 110             );
 111
 112             // Add contributions from subMeshes, if any.
 113             computeCoupledWeights
 114             (
 115                 cIndex,
 116                 cellAlgorithm.dimension(),
 117                 masterObjects,
 118                 cellWeights_[cellI],
 119                 cellCentres_[cellI]
 120             );
 121
 122             // Compute error
 123             scalar error = mag(1.0 - sum(cellWeights_[cellI]));
 124
 125             if (error > matchTol)
 126             {
 127                 bool consistent = false;
 128
 129                 // Check whether any edges lie on boundary patches.
 130                 // These cells can have relaxed weights to account
 131                 // for mild convexity.
 132                 const cell& cellToCheck = cells_[cIndex];
 133
 134                 if (is2D())
 135                 {
 136                     const labelList& parents = cellParents_[cIndex];
 137
 138                     forAll(parents, cI)
 139                     {
 140                         const cell& pCell = polyMesh::cells()[parents[cI]];
 141
 142                         forAll(pCell, fI)
 143                         {
 144                             const face& pFace = polyMesh::faces()[pCell[fI]];
 145
 146                             if (pFace.size() == 3)
 147                             {
 148                                 continue;
 149                             }
 150
 151                             label fP = boundaryMesh().whichPatch(pCell[fI]);
 152
 153                             // Disregard processor patches
 154                             if (getNeighbourProcessor(fP) > -1)
 155                             {
 156                                 continue;
 157                             }
 158
 159                             if (fP > -1)
 160                             {
 161                                 consistent = true;
 162                                 break;
 163                             }
 164                         }
 165
 166                         if (consistent)
 167                         {
 168                             break;
 169                         }
 170                     }
 171                 }
 172                 else
 173                 {
 174                     forAll(cellToCheck, fI)
 175                     {
 176                         const labelList& fE = faceEdges_[cellToCheck[fI]];
 177
 178                         forAll(fE, eI)
 179                         {
 180                             label eP = whichEdgePatch(fE[eI]);
 181
 182                             // Disregard processor patches
 183                             if (getNeighbourProcessor(eP) > -1)
 184                             {
 185                                 continue;
 186                             }
 187
 188                             if (eP > -1)
 189                             {
 190                                 consistent = true;
 191                                 break;
 192                             }
 193                         }
 194
 195                         if (consistent)
 196                         {
 197                             break;
 198                         }
 199                     }
 200                 }
 201
 202                 if (!consistent)
 203                 {
 204                     nInconsistencies++;
 205
 206                     // Add to list
 207                     cellErrors.append(error);
 208
 209                     // Accumulate error stats
 210                     maxCellError = Foam::max(maxCellError, error);
 211
 212                     failedCells.append
 213                     (
 214                         objectMap
 215                         (
 216                             cIndex,
 217                             cellParents_[cIndex]
 218                         )
 219                     );
 220                 }
 221             }
 222         }
 223     }
 224
 225     // Convex-set algorithm for faces
 226     faceSetAlgorithm faceAlgorithm
 227     (
 228         (*this),
 229         oldPoints_,
 230         edges_,
 231         faces_,
 232         cells_,
 233         owner_,
 234         neighbour_
 235     );
 236
 237     // Compute face mapping
 238     for (label faceI = faceStart; faceI < (faceStart + faceSize); faceI++)
 239     {
 240         label fIndex = facesFromFaces_[faceI].index();
 241         labelList& masterObjects = facesFromFaces_[faceI].masterObjects();
 242
 243         label patchIndex = whichPatch(fIndex);
 244         label neiProc = getNeighbourProcessor(patchIndex);
 245
 246         // Skip mapping for internal / processor faces.
 247         if (patchIndex == -1 || neiProc > -1)
 248         {
 249             // Set dummy masters, so that the conventional
 250             // faceMapper doesn't crash-and-burn
 251             masterObjects = labelList(1, 0);
 252
 253             continue;
 254         }
 255
 256         if (skipMapping)
 257         {
 258             // Dummy map from patch[0]
 259             masterObjects = labelList(1, 0);
 260             faceWeights_[faceI].setSize(1, 1.0);
 261             faceCentres_[faceI].setSize(1, vector::zero);
 262         }
 263         else
 264         {
 265             // Obtain weighting factors for this face.
 266             faceAlgorithm.computeWeights
 267             (
 268                 fIndex,
 269                 boundaryMesh()[patchIndex].start(),
 270                 faceParents_[fIndex],
 271                 boundaryMesh()[patchIndex].faceFaces(),
 272                 masterObjects,
 273                 faceWeights_[faceI],
 274                 faceCentres_[faceI]
 275             );
 276
 277             // Add contributions from subMeshes, if any.
 278             computeCoupledWeights
 279             (
 280                 fIndex,
 281                 faceAlgorithm.dimension(),
 282                 masterObjects,
 283                 faceWeights_[faceI],
 284                 faceCentres_[faceI]
 285             );
 286
 287             // Compute error
 288             scalar error = mag(1.0 - sum(faceWeights_[faceI]));
 289
 290             if (error > matchTol)
 291             {
 292                 bool consistent = false;
 293
 294                 // Check whether any edges lie on bounding curves.
 295                 // These faces can have relaxed weights to account
 296                 // for addressing into patches on the other side
 297                 // of the curve.
 298                 const labelList& fEdges = faceEdges_[fIndex];
 299
 300                 forAll(fEdges, eI)
 301                 {
 302                     if (checkBoundingCurve(fEdges[eI]))
 303                     {
 304                         consistent = true;
 305                     }
 306                 }
 307
 308                 if (!consistent)
 309                 {
 310                     nInconsistencies++;
 311
 312                     // Add to list
 313                     faceErrors.append(error);
 314
 315                     // Accumulate error stats
 316                     maxFaceError = Foam::max(maxFaceError, error);
 317
 318                     failedFaces.append
 319                     (
 320                         objectMap
 321                         (
 322                             fIndex,
 323                             faceParents_[fIndex]
 324                         )
 325                     );
 326                 }
 327             }
 328         }
 329     }
 330
 331     if (nInconsistencies)
 332     {
 333         Pout<< " Mapping errors: "
 334             << " max cell error: " << maxCellError
 335             << " max face error: " << maxFaceError
 336             << endl;
 337
 338         if (debug || mappingOutput)
 339         {
 340             if (failedCells.size())
 341             {
 342                 Pout<< " failedCells: " << endl;
 343
 344                 forAll(failedCells, cellI)
 345                 {
 346                     label index = failedCells[cellI].index();
 347
 348                     Pout<< "  Cell: " << index
 349                         << "  Error: " << cellErrors[cellI]
 350                         << endl;
 351                 }
 352             }
 353
 354             if (failedFaces.size())
 355             {
 356                 Pout<< " failedFaces: " << endl;
 357
 358                 forAll(failedFaces, faceI)
 359                 {
 360                     label index = failedFaces[faceI].index();
 361                     label patchIndex = whichPatch(index);
 362
 363                     const polyBoundaryMesh& boundary = boundaryMesh();
 364
 365                     Pout<< "  Face: " << index
 366                         << "  Patch: " << boundary[patchIndex].name()
 367                         << "  Error: " << faceErrors[faceI]
 368                         << endl;
 369                 }
 370             }
 371
 372             if (debug > 3 || mappingOutput)
 373             {
 374                 // Prepare lists
 375                 labelList objects;
 376                 scalarField weights;
 377                 vectorField centres;
 378
 379                 if (failedCells.size())
 380                 {
 381                     forAll(failedCells, cellI)
 382                     {
 383                         label cIndex = failedCells[cellI].index();
 384
 385                         cellAlgorithm.computeWeights
 386                         (
 387                             cIndex,
 388                             0,
 389                             cellParents_[cIndex],
 390                             polyMesh::cellCells(),
 391                             objects,
 392                             weights,
 393                             centres,
 394                             true
 395                         );
 396
 397                         computeCoupledWeights
 398                         (
 399                             cIndex,
 400                             cellAlgorithm.dimension(),
 401                             objects,
 402                             weights,
 403                             centres,
 404                             true
 405                         );
 406
 407                         // Clear lists
 408                         objects.clear();
 409                         weights.clear();
 410                         centres.clear();
 411                     }
 412                 }
 413
 414                 if (failedFaces.size())
 415                 {
 416                     forAll(failedFaces, faceI)
 417                     {
 418                         label fIndex = failedFaces[faceI].index();
 419                         label patchIndex = whichPatch(fIndex);
 420
 421                         const polyBoundaryMesh& boundary = boundaryMesh();
 422
 423                         faceAlgorithm.computeWeights
 424                         (
 425                             fIndex,
 426                             boundary[patchIndex].start(),
 427                             faceParents_[fIndex],
 428                             boundary[patchIndex].faceFaces(),
 429                             objects,
 430                             weights,
 431                             centres,
 432                             true
 433                         );
 434
 435                         computeCoupledWeights
 436                         (
 437                             fIndex,
 438                             faceAlgorithm.dimension(),
 439                             objects,
 440                             weights,
 441                             centres,
 442                             true
 443                         );
 444
 445                         // Clear lists
 446                         objects.clear();
 447                         weights.clear();
 448                         centres.clear();
 449                     }
 450                 }
 451             }
 452         }
 453     }
 454 }
 455
 456
 457 // Static equivalent for multiThreading
 458 void dynamicTopoFvMesh::computeMappingThread(void *argument)
 459 {
 460     // Recast the argument
 461     meshHandler *thread = static_cast<meshHandler*>(argument);
 462
 463     if (thread->slave())
 464     {
 465         thread->sendSignal(meshHandler::START);
 466     }
 467
 468     dynamicTopoFvMesh& mesh = thread->reference();
 469
 470     // Recast the pointers for the argument
 471     scalar& matchTol  = *(static_cast<scalar*>(thread->operator()(0)));
 472     bool& skipMapping = *(static_cast<bool*>(thread->operator()(1)));
 473     bool& mappingOutput = *(static_cast<bool*>(thread->operator()(2)));
 474     label& faceStart = *(static_cast<label*>(thread->operator()(3)));
 475     label& faceSize = *(static_cast<label*>(thread->operator()(4)));
 476     label& cellStart = *(static_cast<label*>(thread->operator()(5)));
 477     label& cellSize = *(static_cast<label*>(thread->operator()(6)));
 478
 479     // Now calculate addressing
 480     mesh.computeMapping
 481     (
 482         matchTol,
 483         skipMapping,
 484         mappingOutput,
 485         faceStart, faceSize,
 486         cellStart, cellSize
 487     );
 488
 489     if (thread->slave())
 490     {
 491         thread->sendSignal(meshHandler::STOP);
 492     }
 493 }
 494
 495
 496 // Routine to invoke threaded mapping
 497 void dynamicTopoFvMesh::threadedMapping
 498 (
 499     scalar matchTol,
 500     bool skipMapping,
 501     bool mappingOutput
 502 )
 503 {
 504     label nThreads = threader_->getNumThreads();
 505
 506     // If mapping is being skipped, issue a warning.
 507     if (skipMapping)
 508     {
 509         Info<< " *** Mapping is being skipped *** " << endl;
 510     }
 511
 512     // Check if single-threaded
 513     if (nThreads == 1)
 514     {
 515         computeMapping
 516         (
 517             matchTol,
 518             skipMapping,
 519             mappingOutput,
 520             0, facesFromFaces_.size(),
 521             0, cellsFromCells_.size()
 522         );
 523
 524         return;
 525     }
 526
 527     // Set one handler per thread
 528     PtrList<meshHandler> hdl(nThreads);
 529
 530     forAll(hdl, i)
 531     {
 532         hdl.set(i, new meshHandler(*this, threader()));
 533     }
 534
 535     // Simple load-balancing scheme
 536     FixedList<label, 2> index(-1);
 537     FixedList<labelList, 2> tStarts(labelList(nThreads, 0));
 538     FixedList<labelList, 2> tSizes(labelList(nThreads, 0));
 539
 540     index[0] = facesFromFaces_.size();
 541     index[1] = cellsFromCells_.size();
 542
 543     if (debug > 2)
 544     {
 545         Pout<< " Mapping Faces: " << index[0] << nl
 546             << " Mapping Cells: " << index[1] << endl;
 547     }
 548
 549     forAll(index, indexI)
 550     {
 551         label j = 0, total = 0;
 552
 553         while (index[indexI]--)
 554         {
 555             tSizes[indexI][(j = tSizes[indexI].fcIndex(j))]++;
 556         }
 557
 558         for (label i = 1; i < tStarts[indexI].size(); i++)
 559         {
 560             tStarts[indexI][i] = tSizes[indexI][i-1] + total;
 561
 562             total += tSizes[indexI][i-1];
 563         }
 564
 565         if (debug > 2)
 566         {
 567             Pout<< " Load starts: " << tStarts[indexI] << nl
 568                 << " Load sizes: " << tSizes[indexI] << endl;
 569         }
 570     }
 571
 572     // Set the argument list for each thread
 573     forAll(hdl, i)
 574     {
 575         // Size up the argument list
 576         hdl[i].setSize(7);
 577
 578         // Set match tolerance
 579         hdl[i].set(0, &matchTol);
 580
 581         // Set the skipMapping flag
 582         hdl[i].set(1, &skipMapping);
 583
 584         // Set the mappingOutput flag
 585         hdl[i].set(2, &mappingOutput);
 586
 587         // Set the start/size indices
 588         hdl[i].set(3, &(tStarts[0][i]));
 589         hdl[i].set(4, &(tSizes[0][i]));
 590         hdl[i].set(5, &(tStarts[1][i]));
 591         hdl[i].set(6, &(tSizes[1][i]));
 592     }
 593
 594     // Prior to multi-threaded operation,
 595     // force calculation of demand-driven data.
 596     polyMesh::cells();
 597     primitiveMesh::cellCells();
 598
 599     const polyBoundaryMesh& boundary = boundaryMesh();
 600
 601     forAll(boundary, patchI)
 602     {
 603         boundary[patchI].faceFaces();
 604     }
 605
 606     // Force calculation of demand-driven data on subMeshes
 607     initCoupledWeights();
 608
 609     // Execute threads in linear sequence
 610     executeThreads(identity(nThreads), hdl, &computeMappingThread);
 611 }
 612
 613
 614 // Set fill-in mapping information for a particular cell
 615 void dynamicTopoFvMesh::setCellMapping
 616 (
 617     const label cIndex,
 618     const labelList& mapCells,
 619     bool addEntry
 620 )
 621 {
 622     if (addEntry)
 623     {
 624         if (debug > 3)
 625         {
 626             Pout<< "Inserting mapping cell: " << cIndex
 627                 << " mapCells: " << mapCells
 628                 << endl;
 629         }
 630
 631         // Insert index into the list, and overwrite if necessary
 632         label index = -1;
 633
 634         forAll(cellsFromCells_, indexI)
 635         {
 636             if (cellsFromCells_[indexI].index() == cIndex)
 637             {
 638                 index = indexI;
 639                 break;
 640             }
 641         }
 642
 643         if (index == -1)
 644         {
 645             meshOps::sizeUpList
 646             (
 647                 objectMap(cIndex, labelList(0)),
 648                 cellsFromCells_
 649             );
 650         }
 651         else
 652         {
 653             cellsFromCells_[index].masterObjects() = labelList(0);
 654         }
 655     }
 656
 657     // Update cell-parents information
 658     dynamicLabelList masterCells(5);
 659
 660     forAll(mapCells, cellI)
 661     {
 662         if (mapCells[cellI] < 0)
 663         {
 664             continue;
 665         }
 666
 667         if (mapCells[cellI] < nOldCells_)
 668         {
 669             if (findIndex(masterCells, mapCells[cellI]) == -1)
 670             {
 671                 masterCells.append(mapCells[cellI]);
 672             }
 673         }
 674
 675         if (cellParents_.found(mapCells[cellI]))
 676         {
 677             const labelList& nParents = cellParents_[mapCells[cellI]];
 678
 679             forAll(nParents, cI)
 680             {
 681                 if (findIndex(masterCells, nParents[cI]) == -1)
 682                 {
 683                     masterCells.append(nParents[cI]);
 684                 }
 685             }
 686         }
 687     }
 688
 689     cellParents_.set(cIndex, masterCells);
 690 }
 691
 692
 693 // Set fill-in mapping information for a particular face
 694 void dynamicTopoFvMesh::setFaceMapping
 695 (
 696     const label fIndex,
 697     const labelList& mapFaces
 698 )
 699 {
 700     label patch = whichPatch(fIndex);
 701     label neiProc = getNeighbourProcessor(patch);
 702
 703     if (debug > 3)
 704     {
 705         const polyBoundaryMesh& boundary = boundaryMesh();
 706
 707         word pName;
 708
 709         if (patch == -1)
 710         {
 711             pName = "Internal";
 712         }
 713         else
 714         if (patch < boundary.size())
 715         {
 716             pName = boundaryMesh()[patch].name();
 717         }
 718         else
 719         {
 720             pName = "New patch: " + Foam::name(patch);
 721         }
 722
 723         Pout<< "Inserting mapping face: " << fIndex
 724             << " patch: " << pName
 725             << " mapFaces: " << mapFaces
 726             << " neiProc: "  << neiProc
 727             << endl;
 728     }
 729
 730     bool foundError = false;
 731
 732     // Check to ensure that internal faces are not mapped
 733     // from any faces in the mesh
 734     if (patch == -1 && mapFaces.size())
 735     {
 736         foundError = true;
 737     }
 738
 739     // Check to ensure that mapFaces is non-empty for physical patches
 740     if (patch > -1 && mapFaces.empty() && neiProc == -1)
 741     {
 742         foundError = true;
 743     }
 744
 745     if (foundError)
 746     {
 747         writeVTK("mapFace_" + Foam::name(fIndex), fIndex, 2);
 748
 749         FatalErrorIn
 750         (
 751             "\n"
 752             "void dynamicTopoFvMesh::setFaceMapping\n"
 753             "(\n"
 754             "    const label fIndex,\n"
 755             "    const labelList& mapFaces\n"
 756             ")"
 757         )
 758             << nl << " Incompatible mapping. " << nl
 759             << "  Possible reasons: " << nl
 760             << "    1. No mapping specified for a boundary face; " << nl
 761             << "    2. Mapping specified for an internal face, " << nl
 762             << "       when none was expected." << nl << nl
 763             << " Face: " << fIndex << nl
 764             << " Patch: "
 765             << (patch > -1 ? boundaryMesh()[patch].name() : "Internal") << nl
 766             << " Owner: " << owner_[fIndex] << nl
 767             << " Neighbour: " << neighbour_[fIndex] << nl
 768             << " mapFaces: " << mapFaces << nl
 769             << abort(FatalError);
 770     }
 771
 772     // Insert addressing into the list, and overwrite if necessary
 773     label index = -1;
 774
 775     forAll(facesFromFaces_, indexI)
 776     {
 777         if (facesFromFaces_[indexI].index() == fIndex)
 778         {
 779             index = indexI;
 780             break;
 781         }
 782     }
 783
 784     if (index == -1)
 785     {
 786         meshOps::sizeUpList
 787         (
 788             objectMap(fIndex, labelList(0)),
 789             facesFromFaces_
 790         );
 791     }
 792     else
 793     {
 794         facesFromFaces_[index].masterObjects() = labelList(0);
 795     }
 796
 797     // For internal / processor faces, bail out
 798     if (patch == -1 || neiProc > -1)
 799     {
 800         return;
 801     }
 802
 803     // Update face-parents information
 804     dynamicLabelList masterFaces(5);
 805
 806     forAll(mapFaces, faceI)
 807     {
 808         if (mapFaces[faceI] < 0)
 809         {
 810             continue;
 811         }
 812
 813         if (mapFaces[faceI] < nOldFaces_)
 814         {
 815             if (findIndex(masterFaces, mapFaces[faceI]) == -1)
 816             {
 817                 masterFaces.append(mapFaces[faceI]);
 818             }
 819         }
 820         else
 821         if (faceParents_.found(mapFaces[faceI]))
 822         {
 823             const labelList& nParents = faceParents_[mapFaces[faceI]];
 824
 825             forAll(nParents, fI)
 826             {
 827                 if (findIndex(masterFaces, nParents[fI]) == -1)
 828                 {
 829                     masterFaces.append(nParents[fI]);
 830                 }
 831             }
 832         }
 833     }
 834
 835     faceParents_.set(fIndex, masterFaces);
 836 }
 837
 838
 839 } // End namespace Foam
 840
 841 // ************************************************************************* //