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