src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshAddressing.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 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 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  45
  46 void meshToMesh::calcAddressing()
  47 {
  48     if (debug)
  49     {
  50         Info<< "meshToMesh::calculateAddressing() : "
  51             << "calculating mesh-to-mesh cell addressing" << endl;
  52     }
  53
  54     // set reference to cells
  55     const cellList& fromCells = fromMesh_.cells();
  56     const pointField& fromPoints = fromMesh_.points();
  57
  58     // In an attempt to preserve the efficiency of linear search and prevent
  59     // failure, a RESCUE mechanism will be set up. Here, we shall mark all
  60     // cells next to the solid boundaries. If such a cell is found as the
  61     // closest, the relationship between the origin and cell will be examined.
  62     // If the origin is outside the cell, a global n-squared search is
  63     // triggered.
  64
  65     // SETTING UP RESCUE
  66
  67     // visit all boundaries and mark the cell next to the boundary.
  68
  69     if (debug)
  70     {
  71         Info<< "meshToMesh::calculateAddressing() : "
  72             << "Setting up rescue" << endl;
  73     }
  74
  75     List<bool> boundaryCell(fromCells.size(), false);
  76
  77     // set reference to boundary
  78     const polyPatchList& patchesFrom = fromMesh_.boundaryMesh();
  79
  80     forAll (patchesFrom, patchI)
  81     {
  82         // get reference to cells next to the boundary
  83         const unallocLabelList& bCells = patchesFrom[patchI].faceCells();
  84
  85         forAll (bCells, faceI)
  86         {
  87             boundaryCell[bCells[faceI]] = true;
  88         }
  89     }
  90
  91     treeBoundBox meshBb(fromPoints);
  92
  93     scalar typDim = meshBb.avgDim()/(2.0*cbrt(scalar(fromCells.size())));
  94
  95     treeBoundBox shiftedBb
  96     (
  97         meshBb.min(),
  98         meshBb.max() + vector(typDim, typDim, typDim)
  99     );
 100
 101     if (debug)
 102     {
 103         Info<< "\nMesh" << endl;
 104         Info<< "   bounding box           : " << meshBb << endl;
 105         Info<< "   bounding box (shifted) : " << shiftedBb << endl;
 106         Info<< "   typical dimension      :" << shiftedBb.typDim() << endl;
 107     }
 108
 109     // Wrap indices and mesh information into helper object
 110     octreeDataCell shapes(fromMesh_);
 111
 112     octree<octreeDataCell> oc
 113     (
 114         shiftedBb,  // overall bounding box
 115         shapes,     // all information needed to do checks on cells
 116         1,          // min. levels
 117         10.0,       // max. size of leaves
 118         10.0        // maximum ratio of cubes v.s. cells
 119     );
 120
 121     if (debug)
 122     {
 123         oc.printStats(Info);
 124     }
 125
 126     cellAddresses
 127     (
 128         cellAddressing_,
 129         toMesh_.cellCentres(),
 130         fromMesh_,
 131         boundaryCell,
 132         oc
 133     );
 134
 135     forAll (toMesh_.boundaryMesh(), patchi)
 136     {
 137         const polyPatch& toPatch = toMesh_.boundaryMesh()[patchi];
 138
 139         if (cuttingPatches_.found(toPatch.name()))
 140         {
 141             boundaryAddressing_[patchi].setSize(toPatch.size());
 142
 143             cellAddresses
 144             (
 145                 boundaryAddressing_[patchi],
 146                 toPatch.faceCentres(),
 147                 fromMesh_,
 148                 boundaryCell,
 149                 oc
 150             );
 151         }
 152         else if
 153         (
 154             patchMap_.found(toPatch.name())
 155          && fromMeshPatches_.found(patchMap_.find(toPatch.name())())
 156         )
 157         {
 158             const polyPatch& fromPatch = fromMesh_.boundaryMesh()
 159             [
 160                 fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
 161             ];
 162
 163             if (fromPatch.empty())
 164             {
 165                 WarningIn("meshToMesh::calcAddressing()")
 166                     << "Source patch " << fromPatch.name()
 167                     << " has no faces. Not performing mapping for it."
 168                     << endl;
 169                 boundaryAddressing_[patchi] = -1;
 170             }
 171             else
 172             {
 173                 treeBoundBox wallBb(fromPatch.localPoints());
 174                 scalar typDim =
 175                     wallBb.avgDim()/(2.0*sqrt(scalar(fromPatch.size())));
 176
 177                 treeBoundBox shiftedBb
 178                 (
 179                     wallBb.min(),
 180                     wallBb.max() + vector(typDim, typDim, typDim)
 181                 );
 182
 183                 // Wrap data for octree into container
 184                 octreeDataFace shapes(fromPatch);
 185
 186                 octree<octreeDataFace> oc
 187                 (
 188                     shiftedBb,  // overall search domain
 189                     shapes,     // all information needed to do checks on cells
 190                     1,          // min levels
 191                     20.0,       // maximum ratio of cubes v.s. cells
 192                     2.0
 193                 );
 194
 195
 196                 const vectorField::subField centresToBoundary =
 197                     toPatch.faceCentres();
 198
 199                 boundaryAddressing_[patchi].setSize(toPatch.size());
 200
 201                 treeBoundBox tightest;
 202                 scalar tightestDist;
 203
 204                 forAll(toPatch, toi)
 205                 {
 206                     tightest = wallBb;                  // starting search bb
 207                     tightestDist = Foam::GREAT;        // starting max distance
 208
 209                     boundaryAddressing_[patchi][toi] = oc.findNearest
 210                     (
 211                         centresToBoundary[toi],
 212                         tightest,
 213                         tightestDist
 214                     );
 215                 }
 216             }
 217         }
 218     }
 219
 220     if (debug)
 221     {
 222         Info<< "meshToMesh::calculateAddressing() : "
 223             << "finished calculating mesh-to-mesh acell ddressing" << endl;
 224     }
 225 }
 226
 227
 228 void meshToMesh::cellAddresses
 229 (
 230     labelList& cellAddressing_,
 231     const pointField& points,
 232     const fvMesh& fromMesh,
 233     const List<bool>& boundaryCell,
 234     const octree<octreeDataCell>& oc
 235 ) const
 236 {
 237     // the implemented search method is a simple neighbour array search.
 238     // It starts from a cell zero, searches its neighbours and finds one
 239     // which is nearer to the target point than the current position.
 240     // The location of the "current position" is reset to that cell and
 241     // search through the neighbours continues. The search is finished
 242     // when all the neighbours of the cell are farther from the target
 243     // point than the current cell
 244
 245     // set curCell label to zero (start)
 246     register label curCell = 0;
 247
 248     // set reference to cell to cell addressing
 249     const vectorField& centresFrom = fromMesh.cellCentres();
 250     const labelListList& cc = fromMesh.cellCells();
 251
 252     forAll (points, toI)
 253     {
 254         // pick up target position
 255         const vector& p = points[toI];
 256
 257         // set the sqr-distance
 258         scalar distSqr = magSqr(p - centresFrom[curCell]);
 259
 260         bool closer;
 261
 262         do
 263         {
 264             closer = false;
 265
 266             // set the current list of neighbouring cells
 267             const labelList& neighbours = cc[curCell];
 268
 269             forAll (neighbours, nI)
 270             {
 271                 scalar curDistSqr =
 272                     magSqr(p - centresFrom[neighbours[nI]]);
 273
 274                 // search through all the neighbours.
 275                 // If the cell is closer, reset current cell and distance
 276                 if (curDistSqr < (1 - SMALL)*distSqr)
 277                 {
 278                     curCell = neighbours[nI];
 279                     distSqr = curDistSqr;
 280                     closer = true;    // a closer neighbour has been found
 281                 }
 282             }
 283         } while (closer);
 284
 285         cellAddressing_[toI] = -1;
 286
 287         // Check point is actually in the nearest cell
 288         if (fromMesh.pointInCell(p, curCell))
 289         {
 290             cellAddressing_[toI] = curCell;
 291         }
 292         else
 293         {
 294             // If curCell is a boundary cell then the point maybe either outside
 295             // the domain or in an other region of the doamin, either way use
 296             // the octree search to find it.
 297             if (boundaryCell[curCell])
 298             {
 299                 cellAddressing_[toI] = oc.find(p);
 300             }
 301             else
 302             {
 303                 // If not on the boundary search the neighbours
 304                 bool found = false;
 305
 306                 // set the current list of neighbouring cells
 307                 const labelList& neighbours = cc[curCell];
 308
 309                 forAll (neighbours, nI)
 310                 {
 311                     // search through all the neighbours.
 312                     // If point is in neighbour reset current cell
 313                     if (fromMesh.pointInCell(p, neighbours[nI]))
 314                     {
 315                         cellAddressing_[toI] = neighbours[nI];
 316                         found = true;
 317                         break;
 318                     }
 319                 }
 320
 321                 if (!found)
 322                 {
 323                     // If still not found search the neighbour-neighbours
 324
 325                     // set the current list of neighbouring cells
 326                     const labelList& neighbours = cc[curCell];
 327
 328                     forAll (neighbours, nI)
 329                     {
 330                         // set the current list of neighbour-neighbouring cells
 331                         const labelList& nn = cc[neighbours[nI]];
 332
 333                         forAll (nn, nI)
 334                         {
 335                             // search through all the neighbours.
 336                             // If point is in neighbour reset current cell
 337                             if (fromMesh.pointInCell(p, nn[nI]))
 338                             {
 339                                 cellAddressing_[toI] = nn[nI];
 340                                 found = true;
 341                                 break;
 342                             }
 343                         }
 344                         if (found) break;
 345                     }
 346                 }
 347
 348                 if (!found)
 349                 {
 350                     // Still not found so us the octree
 351                     cellAddressing_[toI] = oc.find(p);
 352                 }
 353             }
 354         }
 355     }
 356 }
 357
 358
 359 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 360
 361 } // End namespace Foam
 362
 363 // ************************************************************************* //