src/decompositionMethods/decompositionMethods/hierarchGeomDecomp/hierarchGeomDecomp.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "hierarchGeomDecomp.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "PstreamReduceOps.H"
  30 #include "SortableList.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36     defineTypeNameAndDebug(hierarchGeomDecomp, 0);
  37
  38     addToRunTimeSelectionTable
  39     (
  40         decompositionMethod,
  41         hierarchGeomDecomp,
  42         dictionary
  43     );
  44
  45     addToRunTimeSelectionTable
  46     (
  47         decompositionMethod,
  48         hierarchGeomDecomp,
  49         dictionaryMesh
  50     );
  51 }
  52
  53 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  54
  55 void Foam::hierarchGeomDecomp::setDecompOrder()
  56 {
  57     word order(geomDecomDict_.lookup("order"));
  58
  59     if (order.size() != 3)
  60     {
  61         FatalIOErrorIn
  62         (
  63             "hierarchGeomDecomp::hierarchGeomDecomp"
  64             "(const dictionary& decompositionDict)",
  65             decompositionDict_
  66         )   << "number of characters in order (" << order << ") != 3"
  67             << exit(FatalIOError);
  68     }
  69
  70     for (label i = 0; i < 3; i++)
  71     {
  72         if (order[i] == 'x')
  73         {
  74             decompOrder_[i] = 0;
  75         }
  76         else if (order[i] == 'y')
  77         {
  78             decompOrder_[i] = 1;
  79         }
  80         else if (order[i] == 'z')
  81         {
  82             decompOrder_[i] = 2;
  83         }
  84         else
  85         {
  86             FatalIOErrorIn
  87             (
  88                 "hierarchGeomDecomp::hierarchGeomDecomp"
  89                 "(const dictionary& decompositionDict)",
  90                 decompositionDict_
  91             )   << "Illegal decomposition order " << order << endl
  92                 << "It should only contain x, y or z" << exit(FatalError);
  93         }
  94     }
  95 }
  96
  97
  98 Foam::label Foam::hierarchGeomDecomp::findLower
  99 (
 100     const List<scalar>& l,
 101     const scalar t,
 102     const label initLow,
 103     const label initHigh
 104 )
 105 {
 106     if (initHigh <= initLow)
 107     {
 108         return initLow;
 109     }
 110
 111     label low = initLow;
 112     label high = initHigh;
 113
 114     while ((high - low) > 1)
 115     {
 116         label mid = (low + high)/2;
 117
 118         if (l[mid] < t)
 119         {
 120             low = mid;
 121         }
 122         else
 123         {
 124             high = mid;
 125         }
 126     }
 127
 128     // high and low can still differ by one. Choose best.
 129
 130     label tIndex = -1;
 131
 132     if (l[high-1] < t)
 133     {
 134         tIndex = high;
 135     }
 136     else
 137     {
 138         tIndex = low;
 139     }
 140
 141     return tIndex;
 142 }
 143
 144
 145 // Find position in values so between minIndex and this position there
 146 // are wantedSize elements.
 147 void Foam::hierarchGeomDecomp::findBinary
 148 (
 149     const label sizeTol,
 150     const List<scalar>& values,
 151     const label minIndex,       // index of previous value
 152     const scalar minValue,      // value at minIndex
 153     const scalar maxValue,      // global max of values
 154     const scalar wantedSize,    // wanted size
 155
 156     label& mid,                 // index where size of bin is
 157                                 // wantedSize (to within sizeTol)
 158     scalar& midValue            // value at mid
 159 )
 160 {
 161     label low = minIndex;
 162     scalar lowValue = minValue;
 163
 164     scalar highValue = maxValue;
 165     // (one beyond) index of highValue
 166     label high = values.size();
 167
 168     // Safeguards to avoid infinite loop.
 169     label lowPrev = -1;
 170     label midPrev = -1;
 171     label highPrev = -1;
 172
 173     //while (low <= high)
 174     while (true)
 175     {
 176         label size = returnReduce(mid-minIndex, sumOp<label>());
 177
 178         if (debug)
 179         {
 180             Pout<< "low:" << low << " lowValue:" << lowValue
 181                 << " high:" << high << " highValue:" << highValue
 182                 << " mid:" << mid << " midValue:" << midValue << nl
 183                 << "globalSize:" << size << " wantedSize:" << wantedSize
 184                 << " sizeTol:" << sizeTol << endl;
 185         }
 186
 187         if (wantedSize < size - sizeTol)
 188         {
 189             high = mid;
 190             highValue = midValue;
 191         }
 192         else if (wantedSize > size + sizeTol)
 193         {
 194             low = mid;
 195             lowValue = midValue;
 196         }
 197         else
 198         {
 199             break;
 200         }
 201
 202         // Update mid, midValue
 203         midValue = 0.5*(lowValue+highValue);
 204         mid = findLower(values, midValue, low, high);
 205
 206         // Safeguard if same as previous.
 207         bool hasNotChanged = (mid == midPrev) && (low == lowPrev) && (high == highPrev);
 208         if (returnReduce(hasNotChanged, andOp<bool>()))
 209         {
 210             WarningIn("hierarchGeomDecomp::findBinary(..)")
 211                 << "unable to find desired deomposition split, making do!"
 212                 << endl;
 213             break;
 214         }
 215
 216         midPrev = mid;
 217         lowPrev = low;
 218         highPrev = high;
 219     }
 220 }
 221
 222
 223 // Sort points into bins according to one component. Recurses to next component.
 224 void Foam::hierarchGeomDecomp::sortComponent
 225 (
 226     const label sizeTol,
 227     const pointField& points,
 228     const labelList& current,       // slice of points to decompose
 229     const direction componentIndex, // index in decompOrder_
 230     const label mult,               // multiplication factor for finalDecomp
 231     labelList& finalDecomp
 232 )
 233 {
 234     // Current component
 235     label compI = decompOrder_[componentIndex];
 236
 237     if (debug)
 238     {
 239         Pout<< "sortComponent : Sorting slice of size " << current.size()
 240             << " in component " << compI << endl;
 241     }
 242
 243     // Storage for sorted component compI
 244     SortableList<scalar> sortedCoord(current.size());
 245
 246     forAll(current, i)
 247     {
 248         label pointI = current[i];
 249
 250         sortedCoord[i] = points[pointI][compI];
 251     }
 252     sortedCoord.sort();
 253
 254     label globalCurrentSize = returnReduce(current.size(), sumOp<label>());
 255
 256     scalar minCoord = returnReduce
 257     (
 258         (
 259             sortedCoord.size()
 260           ? sortedCoord[0]
 261           : GREAT
 262         ),
 263         minOp<scalar>()
 264     );
 265
 266     scalar maxCoord = returnReduce
 267     (
 268         (
 269             sortedCoord.size()
 270           ? sortedCoord[sortedCoord.size()-1]
 271           : -GREAT
 272         ),
 273         maxOp<scalar>()
 274     );
 275
 276     if (debug)
 277     {
 278         Pout<< "sortComponent : minCoord:" << minCoord
 279             << " maxCoord:" << maxCoord << endl;
 280     }
 281
 282     // starting index (in sortedCoord) of bin (= local)
 283     label leftIndex = 0;
 284     // starting value of bin (= global since coordinate)
 285     scalar leftCoord = minCoord;
 286
 287     // Sort bins of size n
 288     for (label bin = 0; bin < n_[compI]; bin++)
 289     {
 290         // Now we need to determine the size of the bin (dx). This is
 291         // determined by the 'pivot' values - everything to the left of this
 292         // value goes in the current bin, everything to the right into the next
 293         // bins.
 294
 295         // Local number of elements
 296         label localSize = -1;     // offset from leftOffset
 297
 298         // Value at right of bin (leftIndex+localSize-1)
 299         scalar rightCoord = -GREAT;
 300
 301         if (bin == n_[compI]-1)
 302         {
 303             // Last bin. Copy all.
 304             localSize = current.size()-leftIndex;
 305             rightCoord = maxCoord;                  // note: not used anymore
 306         }
 307         else if (Pstream::nProcs() == 1)
 308         {
 309             // No need for binary searching of bin size
 310             localSize = label(current.size()/n_[compI]);
 311             rightCoord = sortedCoord[leftIndex+localSize];
 312         }
 313         else
 314         {
 315             // For the current bin (starting at leftCoord) we want a rightCoord
 316             // such that the sum of all sizes are globalCurrentSize/n_[compI].
 317             // We have to iterate to obtain this.
 318
 319             label rightIndex = current.size();
 320             rightCoord = maxCoord;
 321
 322             // Calculate rightIndex/rightCoord to have wanted size
 323             findBinary
 324             (
 325                 sizeTol,
 326                 sortedCoord,
 327                 leftIndex,
 328                 leftCoord,
 329                 maxCoord,
 330                 globalCurrentSize/n_[compI],  // wanted size
 331
 332                 rightIndex,
 333                 rightCoord
 334             );
 335             localSize = rightIndex - leftIndex;
 336         }
 337
 338         if (debug)
 339         {
 340             Pout<< "For component " << compI << ", bin " << bin
 341                 << " copying" << nl
 342                 << "from " << leftCoord << " at local index "
 343                 << leftIndex << nl
 344                 << "to " << rightCoord << " localSize:"
 345                 << localSize << nl
 346                 << endl;
 347         }
 348
 349
 350         // Copy localSize elements starting from leftIndex.
 351         labelList slice(localSize);
 352
 353         forAll(slice, i)
 354         {
 355             label pointI = current[sortedCoord.indices()[leftIndex+i]];
 356
 357             // Mark point into correct bin
 358             finalDecomp[pointI] += bin*mult;
 359
 360             // And collect for next sorting action
 361             slice[i] = pointI;
 362         }
 363
 364         // Sort slice in next component
 365         if (componentIndex < 2)
 366         {
 367             string oldPrefix;
 368             if (debug)
 369             {
 370                 oldPrefix = Pout.prefix();
 371                 Pout.prefix() = "  " + oldPrefix;
 372             }
 373
 374             sortComponent
 375             (
 376                 sizeTol,
 377                 points,
 378                 slice,
 379                 componentIndex+1,
 380                 mult*n_[compI],     // Multiplier to apply to decomposition.
 381                 finalDecomp
 382             );
 383
 384             if (debug)
 385             {
 386                 Pout.prefix() = oldPrefix;
 387             }
 388         }
 389
 390         // Step to next bin.
 391         leftIndex += localSize;
 392         leftCoord = rightCoord;
 393     }
 394 }
 395
 396
 397 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 398
 399 Foam::hierarchGeomDecomp::hierarchGeomDecomp
 400 (
 401     const dictionary& decompositionDict
 402 )
 403 :
 404     geomDecomp(decompositionDict, typeName),
 405     decompOrder_()
 406 {
 407     setDecompOrder();
 408 }
 409
 410
 411 Foam::hierarchGeomDecomp::hierarchGeomDecomp
 412 (
 413     const dictionary& decompositionDict,
 414     const polyMesh&
 415 )
 416 :
 417     geomDecomp(decompositionDict, hierarchGeomDecomp::typeName),
 418     decompOrder_()
 419 {
 420     setDecompOrder();
 421 }
 422
 423
 424 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 425
 426 Foam::labelList Foam::hierarchGeomDecomp::decompose
 427 (
 428     const pointField& points
 429 )
 430 {
 431     // construct a list for the final result
 432     labelList finalDecomp(points.size(), 0);
 433
 434     // Start off with every point sorted onto itself.
 435     labelList slice(points.size());
 436     forAll(slice, i)
 437     {
 438         slice[i] = i;
 439     }
 440
 441     pointField rotatedPoints = rotDelta_ & points;
 442
 443
 444     // Calculate tolerance of cell distribution. For large cases finding
 445     // distibution to the cell exact would cause too many iterations so allow
 446     // some slack.
 447     label allSize = points.size();
 448     reduce(allSize, sumOp<label>());
 449
 450     const label sizeTol = max(1, label(1E-3*allSize/nProcessors_));
 451
 452     // Sort recursive
 453     sortComponent
 454     (
 455         sizeTol,
 456         rotatedPoints,
 457         slice,
 458         0,              // Sort first component in decompOrder.
 459         1,              // Offset for different x bins.
 460         finalDecomp
 461     );
 462
 463     return finalDecomp;
 464 }
 465
 466
 467 // ************************************************************************* //