src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Description
  26     private member of meshToMesh.
  27     Calculates mesh to mesh addressing pattern (for each cell from one mesh
  28     find the closest cell centre in the other mesh).
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "meshToMesh.H"
  33 #include "SubField.H"
  34
  35 #include "octree.H"
  36 #include "octreeDataCell.H"
  37 #include "octreeDataFace.H"
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 const scalar meshToMesh::cellCentreDistanceTol
  45 (
  46     debug::tolerances("meshToMeshCellCentreDistanceTol", 1e-3)
  47 );
  48
  49
  50 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  51
  52 void meshToMesh::calcAddressing()
  53 {
  54     if (debug)
  55     {
  56         Info<< "meshToMesh::calculateAddressing() : "
  57             << "calculating mesh-to-mesh cell addressing" << endl;
  58     }
  59
  60     // set reference to cells
  61     const cellList& fromCells = fromMesh_.cells();
  62     const pointField& fromPoints = fromMesh_.points();
  63
  64     // In an attempt to preserve the efficiency of linear search and prevent
  65     // failure, a RESCUE mechanism will be set up. Here, we shall mark all
  66     // cells next to the solid boundaries. If such a cell is found as the
  67     // closest, the relationship between the origin and cell will be examined.
  68     // If the origin is outside the cell, a global n-squared search is
  69     // triggered.
  70
  71     // SETTING UP RESCUE
  72
  73     // visit all boundaries and mark the cell next to the boundary.
  74
  75     if (debug)
  76     {
  77         Info<< "meshToMesh::calculateAddressing() : "
  78             << "Setting up rescue" << endl;
  79     }
  80
  81     List<bool> boundaryCell(fromCells.size(), false);
  82
  83     // set reference to boundary
  84     const polyPatchList& patchesFrom = fromMesh_.boundaryMesh();
  85
  86     forAll (patchesFrom, patchI)
  87     {
  88         // get reference to cells next to the boundary
  89         const unallocLabelList& bCells = patchesFrom[patchI].faceCells();
  90
  91         forAll (bCells, faceI)
  92         {
  93             boundaryCell[bCells[faceI]] = true;
  94         }
  95     }
  96
  97     treeBoundBox meshBb(fromPoints);
  98
  99     scalar typDim = meshBb.avgDim()/(2.0*cbrt(scalar(fromCells.size())));
 100
 101     treeBoundBox shiftedBb
 102     (
 103         meshBb.min(),
 104         meshBb.max() + vector(typDim, typDim, typDim)
 105     );
 106
 107     if (debug)
 108     {
 109         Info<< "\nMesh" << endl;
 110         Info<< "   bounding box           : " << meshBb << endl;
 111         Info<< "   bounding box (shifted) : " << shiftedBb << endl;
 112         Info<< "   typical dimension      :" << shiftedBb.typDim() << endl;
 113     }
 114
 115     // Wrap indices and mesh information into helper object
 116     octreeDataCell shapes(fromMesh_);
 117
 118     octree<octreeDataCell> oc
 119     (
 120         shiftedBb,  // overall bounding box
 121         shapes,     // all information needed to do checks on cells
 122         1,          // min. levels
 123         10.0,       // max. size of leaves
 124         10.0        // maximum ratio of cubes v.s. cells
 125     );
 126
 127     if (debug)
 128     {
 129         oc.printStats(Info);
 130     }
 131
 132     cellAddresses
 133     (
 134         cellAddressing_,
 135         toMesh_.cellCentres(),
 136         fromMesh_,
 137         boundaryCell,
 138         oc,
 139         true
 140     );
 141
 142     forAll (toMesh_.boundaryMesh(), patchi)
 143     {
 144
 145         const polyPatch& toPatch = toMesh_.boundaryMesh()[patchi];
 146
 147         if (cuttingPatches_.found(toPatch.name()))
 148         {
 149
 150             boundaryAddressing_[patchi].setSize(toPatch.size());
 151
 152             cellAddresses
 153             (
 154                 boundaryAddressing_[patchi],
 155                 toPatch.faceCentres(),
 156                 fromMesh_,
 157                 boundaryCell,
 158                 oc,
 159                 false
 160             );
 161         }
 162         else if
 163         (
 164             patchMap_.found(toPatch.name())
 165          && fromMeshPatches_.found(patchMap_.find(toPatch.name())())
 166         )
 167         {
 168             const polyPatch& fromPatch = fromMesh_.boundaryMesh()
 169             [
 170                 fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
 171             ];
 172
 173             if (fromPatch.empty())
 174             {
 175                 WarningIn("meshToMesh::calcAddressing()")
 176                     << "Source patch " << fromPatch.name()
 177                     << " has no faces. Not performing mapping for it."
 178                     << endl;
 179                 boundaryAddressing_[patchi] = -1;
 180             }
 181             else
 182             {
 183                 treeBoundBox wallBb(fromPatch.localPoints());
 184                 scalar typDim =
 185                     wallBb.avgDim()/(2.0*sqrt(scalar(fromPatch.size())));
 186
 187                 treeBoundBox shiftedBb
 188                 (
 189                     wallBb.min(),
 190                     wallBb.max() + vector(typDim, typDim, typDim)
 191                 );
 192
 193                 // Wrap data for octree into container
 194                 octreeDataFace shapes(fromPatch);
 195
 196                 octree<octreeDataFace> oc
 197                 (
 198                     shiftedBb,  // overall search domain
 199                     shapes,     // all information needed to do checks on cells
 200                     1,          // min levels
 201                     20.0,       // maximum ratio of cubes v.s. cells
 202                     2.0
 203                 );
 204
 205
 206                 const vectorField::subField centresToBoundary =
 207                     toPatch.faceCentres();
 208
 209                 boundaryAddressing_[patchi].setSize(toPatch.size());
 210
 211                 treeBoundBox tightest;
 212                 scalar tightestDist;
 213
 214                 forAll(toPatch, toi)
 215                 {
 216                     tightest = wallBb;                  // starting search bb
 217                     tightestDist = Foam::GREAT;        // starting max distance
 218
 219                     boundaryAddressing_[patchi][toi] = oc.findNearest
 220                     (
 221                         centresToBoundary[toi],
 222                         tightest,
 223                         tightestDist
 224                     );
 225                 }
 226             }
 227         }
 228     }
 229
 230     if (debug)
 231     {
 232         Info<< "meshToMesh::calculateAddressing() : "
 233             << "finished calculating mesh-to-mesh acell ddressing" << endl;
 234     }
 235 }
 236
 237
 238 void meshToMesh::cellAddresses
 239 (
 240     labelList& cellAddressing_,
 241     const pointField& points,
 242     const fvMesh& fromMesh,
 243     const List<bool>& boundaryCell,
 244     const octree<octreeDataCell>& oc,
 245     bool forceFind
 246 ) const
 247 {
 248
 249     label nCellsOutsideAddressing = 0;
 250
 251     // the implemented search method is a simple neighbour array search.
 252     // It starts from a cell zero, searches its neighbours and finds one
 253     // which is nearer to the target point than the current position.
 254     // The location of the "current position" is reset to that cell and
 255     // search through the neighbours continues. The search is finished
 256     // when all the neighbours of the cell are farther from the target
 257     // point than the current cell
 258
 259     // set curCell label to zero (start)
 260     register label curCell = 0;
 261
 262     // set reference to cell to cell addressing
 263     const vectorField& centresFrom = fromMesh.cellCentres();
 264     const labelListList& cc = fromMesh.cellCells();
 265
 266     forAll (points, toI)
 267     {
 268
 269         scalar localTol = cellCentreDistanceTol;
 270
 271         bool isBoundary = false;
 272
 273         // pick up target position
 274         const vector& p = points[toI];
 275
 276         // set the sqr-distance
 277         scalar distSqr = magSqr(p - centresFrom[curCell]);
 278
 279         bool closer;
 280
 281         do
 282         {
 283             closer = false;
 284
 285             // set the current list of neighbouring cells
 286             const labelList& neighbours = cc[curCell];
 287
 288             forAll (neighbours, nI)
 289             {
 290                 scalar curDistSqr =
 291                     magSqr(p - centresFrom[neighbours[nI]]);
 292
 293                 // search through all the neighbours.
 294                 // If the cell is closer, reset current cell and distance
 295                 if (curDistSqr < (1 - SMALL)*distSqr)
 296                 {
 297                     curCell = neighbours[nI];
 298                     distSqr = curDistSqr;
 299                     closer = true;    // a closer neighbour has been found
 300                 }
 301             }
 302         } while (closer);
 303
 304         cellAddressing_[toI] = -1;
 305
 306         // Check point is actually in the nearest cell
 307         if (fromMesh.pointInCell(p, curCell))
 308         {
 309             cellAddressing_[toI] = curCell;
 310         }
 311         else
 312         {
 313             // If curCell is a boundary cell then the point maybe either
 314             // outside the domain or in an other region of the doamin,
 315             //  either way use the octree search to find it.
 316             if (boundaryCell[curCell])
 317             {
 318                 isBoundary = true;
 319                 cellAddressing_[toI] = oc.find(p);
 320             }
 321             else
 322             {
 323                 // If not on the boundary search the neighbours
 324                 bool found = false;
 325
 326                 // set the current list of neighbouring cells
 327                 const labelList& neighbours = cc[curCell];
 328
 329                 forAll (neighbours, nI)
 330                 {
 331                     // search through all the neighbours.
 332                     // If point is in neighbour reset current cell
 333                     if (fromMesh.pointInCell(p, neighbours[nI]))
 334                     {
 335                         cellAddressing_[toI] = neighbours[nI];
 336                         found = true;
 337                         break;
 338                     }
 339                 }
 340
 341                 if (!found)
 342                 {
 343                     // If still not found search the neighbour-neighbours
 344
 345                     // set the current list of neighbouring cells
 346                     const labelList& neighbours = cc[curCell];
 347
 348                     forAll (neighbours, nI)
 349                     {
 350                         // set the current list of neighbour-neighbouring cells
 351                         const labelList& nn = cc[neighbours[nI]];
 352
 353                         forAll (nn, nI)
 354                         {
 355                             // search through all the neighbours.
 356                             // If point is in neighbour reset current cell
 357                             if (fromMesh.pointInCell(p, nn[nI]))
 358                             {
 359                                 cellAddressing_[toI] = nn[nI];
 360                                 found = true;
 361                                 break;
 362                             }
 363                         }
 364                         if (found) break;
 365                     }
 366                 }
 367
 368                 if (!found)
 369                 {
 370                     // Still not found so use the octree
 371                     cellAddressing_[toI] = oc.find(p);
 372                 }
 373             }
 374
 375             if (cellAddressing_[toI] < 0)
 376             {
 377                 nCellsOutsideAddressing++;
 378
 379                 if (isBoundary && forceFind)
 380                 {
 381                     // If still not found, get the closest cell within the
 382                     // specified tolerance
 383
 384                     forAll(fromMesh.boundary(), patchi)
 385                     {
 386                         const fvPatch& patch = fromMesh.boundary()[patchi];
 387
 388                         word name = patch.name();
 389
 390                         label patchID =
 391                             toMesh_.boundaryMesh().findPatchID(name);
 392
 393                         label sizePatch = 0;
 394                         if (patchID > -1)
 395                         {
 396                             sizePatch = toMesh_.boundary()[patchID].size();
 397                         }
 398
 399                         if
 400                         (
 401                             sizePatch > 0
 402                         )
 403                         {
 404                             forAll(patch, facei)
 405                             {
 406                                 label celli = patch.faceCells()[facei];
 407
 408                                 const vector& centre = fromMesh.C()[celli];
 409                                 if (mag(points[toI] - centre) < localTol)
 410                                 {
 411                                     localTol = mag(points[toI] - centre);
 412                                     cellAddressing_[toI] = celli;
 413                                 }
 414
 415                             }
 416                         }
 417
 418                     }
 419                 }
 420             }
 421         }
 422     }
 423
 424     if (nCellsOutsideAddressing > 0)
 425     {
 426         Info<< "Found " << nCellsOutsideAddressing
 427             << " cells outside of the addressing" << nl
 428             << "Cell addressing size = " << cellAddressing_.size() << endl;
 429     }
 430 }
 431
 432
 433 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 434
 435 } // End namespace Foam
 436
 437 // ************************************************************************* //