src/meshTools/searchableSurface/searchableSurfaceCollection.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 "searchableSurfaceCollection.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "SortableList.H"
  30 #include "Time.H"
  31 #include "ListOps.H"
  32
  33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 defineTypeNameAndDebug(searchableSurfaceCollection, 0);
  39 addToRunTimeSelectionTable
  40 (
  41     searchableSurface,
  42     searchableSurfaceCollection,
  43     dict
  44 );
  45
  46 }
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 void Foam::searchableSurfaceCollection::findNearest
  52 (
  53     const pointField& samples,
  54     scalarField& minDistSqr,
  55     List<pointIndexHit>& nearestInfo,
  56     labelList& nearestSurf
  57 ) const
  58 {
  59     // Initialise
  60     nearestInfo.setSize(samples.size());
  61     nearestInfo = pointIndexHit();
  62     nearestSurf.setSize(samples.size());
  63     nearestSurf = -1;
  64
  65     List<pointIndexHit> hitInfo(samples.size());
  66
  67     const scalarField localMinDistSqr(samples.size(), GREAT);
  68
  69     forAll(subGeom_, surfI)
  70     {
  71         subGeom_[surfI].findNearest
  72         (
  73             cmptDivide  // Transform then divide
  74             (
  75                 transform_[surfI].localPosition(samples),
  76                 scale_[surfI]
  77             ),
  78             localMinDistSqr,
  79             hitInfo
  80         );
  81
  82         forAll(hitInfo, pointI)
  83         {
  84             if (hitInfo[pointI].hit())
  85             {
  86                 // Rework back into global coordinate sys. Multiply then
  87                 // transform
  88                 point globalPt = transform_[surfI].globalPosition
  89                 (
  90                     cmptMultiply
  91                     (
  92                         hitInfo[pointI].rawPoint(),
  93                         scale_[surfI]
  94                     )
  95                 );
  96
  97                 scalar distSqr = magSqr(globalPt - samples[pointI]);
  98
  99                 if (distSqr < minDistSqr[pointI])
 100                 {
 101                     minDistSqr[pointI] = distSqr;
 102                     nearestInfo[pointI].setPoint(globalPt);
 103                     nearestInfo[pointI].setHit();
 104                     nearestInfo[pointI].setIndex
 105                     (
 106                         hitInfo[pointI].index()
 107                       + indexOffset_[surfI]
 108                     );
 109                     nearestSurf[pointI] = surfI;
 110                 }
 111             }
 112         }
 113     }
 114 }
 115
 116
 117 // Sort hits into per-surface bins. Misses are rejected. Maintains map back
 118 // to position
 119 void Foam::searchableSurfaceCollection::sortHits
 120 (
 121     const List<pointIndexHit>& info,
 122     List<List<pointIndexHit> >& surfInfo,
 123     labelListList& infoMap
 124 ) const
 125 {
 126     // Count hits per surface.
 127     labelList nHits(subGeom_.size(), 0);
 128
 129     forAll(info, pointI)
 130     {
 131         if (info[pointI].hit())
 132         {
 133             label index = info[pointI].index();
 134             label surfI = findLower(indexOffset_, index+1);
 135             nHits[surfI]++;
 136         }
 137     }
 138
 139     // Per surface the hit
 140     surfInfo.setSize(subGeom_.size());
 141     // Per surface the original position
 142     infoMap.setSize(subGeom_.size());
 143
 144     forAll(surfInfo, surfI)
 145     {
 146         surfInfo[surfI].setSize(nHits[surfI]);
 147         infoMap[surfI].setSize(nHits[surfI]);
 148     }
 149     nHits = 0;
 150
 151     forAll(info, pointI)
 152     {
 153         if (info[pointI].hit())
 154         {
 155             label index = info[pointI].index();
 156             label surfI = findLower(indexOffset_, index+1);
 157
 158             // Store for correct surface and adapt indices back to local
 159             // ones
 160             label localI = nHits[surfI]++;
 161             surfInfo[surfI][localI] = pointIndexHit
 162             (
 163                 info[pointI].hit(),
 164                 info[pointI].rawPoint(),
 165                 index-indexOffset_[surfI]
 166             );
 167             infoMap[surfI][localI] = pointI;
 168         }
 169     }
 170 }
 171
 172
 173 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 174
 175 Foam::searchableSurfaceCollection::searchableSurfaceCollection
 176 (
 177     const IOobject& io,
 178     const dictionary& dict
 179 )
 180 :
 181     searchableSurface(io),
 182     instance_(dict.size()),
 183     scale_(dict.size()),
 184     transform_(dict.size()),
 185     subGeom_(dict.size()),
 186     mergeSubRegions_(dict.lookup("mergeSubRegions")),
 187     indexOffset_(dict.size()+1)
 188 {
 189     Info<< "SearchableCollection : " << name() << endl;
 190
 191     label surfI = 0;
 192     label startIndex = 0;
 193     forAllConstIter(dictionary, dict, iter)
 194     {
 195         if (dict.isDict(iter().keyword()))
 196         {
 197             instance_[surfI] = iter().keyword();
 198
 199             const dictionary& subDict = dict.subDict(instance_[surfI]);
 200
 201             scale_[surfI] = subDict.lookup("scale");
 202             transform_.set
 203             (
 204                 surfI,
 205                 coordinateSystem::New
 206                 (
 207                     "",
 208                     subDict.subDict("transform")
 209                 )
 210             );
 211
 212             const word subGeomName(subDict.lookup("surface"));
 213             //Pout<< "Trying to find " << subGeomName << endl;
 214
 215             const searchableSurface& s =
 216                 io.db().lookupObject<searchableSurface>(subGeomName);
 217
 218             // I don't know yet how to handle the globalSize combined with
 219             // regionOffset. Would cause non-consecutive indices locally
 220             // if all indices offset by globalSize() of the local region...
 221             if (s.size() != s.globalSize())
 222             {
 223                 FatalErrorIn
 224                 (
 225                     "searchableSurfaceCollection::searchableSurfaceCollection"
 226                     "(const IOobject&, const dictionary&)"
 227                 )   << "Cannot use a distributed surface in a collection."
 228                     << exit(FatalError);
 229             }
 230
 231             subGeom_.set(surfI, &const_cast<searchableSurface&>(s));
 232
 233             indexOffset_[surfI] = startIndex;
 234             startIndex += subGeom_[surfI].size();
 235
 236             Info<< "    instance : " << instance_[surfI] << endl;
 237             Info<< "    surface  : " << s.name() << endl;
 238             Info<< "    scale    : " << scale_[surfI] << endl;
 239             Info<< "    coordsys : " << transform_[surfI] << endl;
 240
 241             surfI++;
 242         }
 243     }
 244     indexOffset_[surfI] = startIndex;
 245
 246     instance_.setSize(surfI);
 247     scale_.setSize(surfI);
 248     transform_.setSize(surfI);
 249     subGeom_.setSize(surfI);
 250     indexOffset_.setSize(surfI+1);
 251 }
 252
 253
 254 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 255
 256 Foam::searchableSurfaceCollection::~searchableSurfaceCollection()
 257 {}
 258
 259
 260 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 261
 262 const Foam::wordList& Foam::searchableSurfaceCollection::regions() const
 263 {
 264     if (regions_.size() == 0)
 265     {
 266         regionOffset_.setSize(subGeom_.size());
 267
 268         DynamicList<word> allRegions;
 269         forAll(subGeom_, surfI)
 270         {
 271             regionOffset_[surfI] = allRegions.size();
 272
 273             if (mergeSubRegions_)
 274             {
 275                 // Single name regardless how many regions subsurface has
 276                 allRegions.append(instance_[surfI] + "_" + Foam::name(surfI));
 277             }
 278             else
 279             {
 280                 const wordList& subRegions = subGeom_[surfI].regions();
 281
 282                 forAll(subRegions, i)
 283                 {
 284                     allRegions.append(instance_[surfI] + "_" + subRegions[i]);
 285                 }
 286             }
 287         }
 288         regions_.transfer(allRegions.shrink());
 289     }
 290     return regions_;
 291 }
 292
 293
 294 Foam::label Foam::searchableSurfaceCollection::size() const
 295 {
 296     return indexOffset_[indexOffset_.size()-1];
 297 }
 298
 299
 300 Foam::pointField Foam::searchableSurfaceCollection::coordinates() const
 301 {
 302     // Get overall size
 303     pointField coords(size());
 304
 305     // Append individual coordinates
 306     label coordI = 0;
 307
 308     forAll(subGeom_, surfI)
 309     {
 310         const pointField subCoords = subGeom_[surfI].coordinates();
 311
 312         forAll(subCoords, i)
 313         {
 314             coords[coordI++] = transform_[surfI].globalPosition
 315             (
 316                 cmptMultiply
 317                 (
 318                     subCoords[i],
 319                     scale_[surfI]
 320                 )
 321             );
 322         }
 323     }
 324
 325     return coords;
 326 }
 327
 328
 329 void Foam::searchableSurfaceCollection::findNearest
 330 (
 331     const pointField& samples,
 332     const scalarField& nearestDistSqr,
 333     List<pointIndexHit>& nearestInfo
 334 ) const
 335 {
 336     // How to scale distance?
 337     scalarField minDistSqr(nearestDistSqr);
 338
 339     labelList nearestSurf;
 340     findNearest
 341     (
 342         samples,
 343         minDistSqr,
 344         nearestInfo,
 345         nearestSurf
 346     );
 347 }
 348
 349
 350 void Foam::searchableSurfaceCollection::findLine
 351 (
 352     const pointField& start,
 353     const pointField& end,
 354     List<pointIndexHit>& info
 355 ) const
 356 {
 357     info.setSize(start.size());
 358     info = pointIndexHit();
 359
 360     // Current nearest (to start) intersection
 361     pointField nearest(end);
 362
 363     List<pointIndexHit> hitInfo(start.size());
 364
 365     forAll(subGeom_, surfI)
 366     {
 367         // Starting point
 368         tmp<pointField> e0 = cmptDivide
 369         (
 370             transform_[surfI].localPosition
 371             (
 372                 start
 373             ),
 374             scale_[surfI]
 375         );
 376
 377         // Current best end point
 378         tmp<pointField> e1 = cmptDivide
 379         (
 380             transform_[surfI].localPosition
 381             (
 382                 nearest
 383             ),
 384             scale_[surfI]
 385         );
 386
 387         subGeom_[surfI].findLine(e0, e1, hitInfo);
 388
 389         forAll(hitInfo, pointI)
 390         {
 391             if (hitInfo[pointI].hit())
 392             {
 393                 // Transform back to global coordinate sys.
 394                 nearest[pointI] = transform_[surfI].globalPosition
 395                 (
 396                     cmptMultiply
 397                     (
 398                         hitInfo[pointI].rawPoint(),
 399                         scale_[surfI]
 400                     )
 401                 );
 402                 info[pointI] = hitInfo[pointI];
 403                 info[pointI].rawPoint() = nearest[pointI];
 404                 info[pointI].setIndex
 405                 (
 406                     hitInfo[pointI].index()
 407                   + indexOffset_[surfI]
 408                 );
 409             }
 410         }
 411     }
 412
 413
 414     // Debug check
 415     if (false)
 416     {
 417         forAll(info, pointI)
 418         {
 419             if (info[pointI].hit())
 420             {
 421                 vector n(end[pointI] - start[pointI]);
 422                 scalar magN = mag(n);
 423
 424                 if (magN > SMALL)
 425                 {
 426                     n /= mag(n);
 427
 428                     scalar s = ((info[pointI].rawPoint()-start[pointI])&n);
 429
 430                     if (s < 0 || s > 1)
 431                     {
 432                         FatalErrorIn
 433                         (
 434                             "searchableSurfaceCollection::findLine(..)"
 435                         )   << "point:" << info[pointI]
 436                             << " s:" << s
 437                             << " outside vector "
 438                             << " start:" << start[pointI]
 439                             << " end:" << end[pointI]
 440                             << abort(FatalError);
 441                     }
 442                 }
 443             }
 444         }
 445     }
 446 }
 447
 448
 449 void Foam::searchableSurfaceCollection::findLineAny
 450 (
 451     const pointField& start,
 452     const pointField& end,
 453     List<pointIndexHit>& info
 454 ) const
 455 {
 456     // To be done ...
 457     findLine(start, end, info);
 458 }
 459
 460
 461 void Foam::searchableSurfaceCollection::findLineAll
 462 (
 463     const pointField& start,
 464     const pointField& end,
 465     List<List<pointIndexHit> >& info
 466 ) const
 467 {
 468     // To be done. Assume for now only one intersection.
 469     List<pointIndexHit> nearestInfo;
 470     findLine(start, end, nearestInfo);
 471
 472     info.setSize(start.size());
 473     forAll(info, pointI)
 474     {
 475         if (nearestInfo[pointI].hit())
 476         {
 477             info[pointI].setSize(1);
 478             info[pointI][0] = nearestInfo[pointI];
 479         }
 480         else
 481         {
 482             info[pointI].clear();
 483         }
 484     }
 485 }
 486
 487
 488 void Foam::searchableSurfaceCollection::getRegion
 489 (
 490     const List<pointIndexHit>& info,
 491     labelList& region
 492 ) const
 493 {
 494     if (subGeom_.size() == 0)
 495     {}
 496     else if (subGeom_.size() == 1)
 497     {
 498         if (mergeSubRegions_)
 499         {
 500             region.setSize(info.size());
 501             region = regionOffset_[0];
 502         }
 503         else
 504         {
 505             subGeom_[0].getRegion(info, region);
 506         }
 507     }
 508     else
 509     {
 510         // Multiple surfaces. Sort by surface.
 511
 512         // Per surface the hit
 513         List<List<pointIndexHit> > surfInfo;
 514         // Per surface the original position
 515         List<List<label> > infoMap;
 516         sortHits(info, surfInfo, infoMap);
 517
 518         region.setSize(info.size());
 519         region = -1;
 520
 521         // Do region tests
 522
 523         if (mergeSubRegions_)
 524         {
 525             // Actually no need for surfInfo. Just take region for surface.
 526             forAll(infoMap, surfI)
 527             {
 528                 const labelList& map = infoMap[surfI];
 529                 forAll(map, i)
 530                 {
 531                     region[map[i]] = regionOffset_[surfI];
 532                 }
 533             }
 534         }
 535         else
 536         {
 537             forAll(infoMap, surfI)
 538             {
 539                 labelList surfRegion;
 540                 subGeom_[surfI].getRegion(surfInfo[surfI], surfRegion);
 541
 542                 const labelList& map = infoMap[surfI];
 543                 forAll(map, i)
 544                 {
 545                     region[map[i]] = regionOffset_[surfI] + surfRegion[i];
 546                 }
 547             }
 548         }
 549     }
 550 }
 551
 552
 553 void Foam::searchableSurfaceCollection::getNormal
 554 (
 555     const List<pointIndexHit>& info,
 556     vectorField& normal
 557 ) const
 558 {
 559     if (subGeom_.size() == 0)
 560     {}
 561     else if (subGeom_.size() == 1)
 562     {
 563         subGeom_[0].getNormal(info, normal);
 564     }
 565     else
 566     {
 567         // Multiple surfaces. Sort by surface.
 568
 569         // Per surface the hit
 570         List<List<pointIndexHit> > surfInfo;
 571         // Per surface the original position
 572         List<List<label> > infoMap;
 573         sortHits(info, surfInfo, infoMap);
 574
 575         normal.setSize(info.size());
 576
 577         // Do region tests
 578         forAll(surfInfo, surfI)
 579         {
 580             vectorField surfNormal;
 581             subGeom_[surfI].getNormal(surfInfo[surfI], surfNormal);
 582
 583             const labelList& map = infoMap[surfI];
 584             forAll(map, i)
 585             {
 586                 normal[map[i]] = surfNormal[i];
 587             }
 588         }
 589     }
 590 }
 591
 592
 593 void Foam::searchableSurfaceCollection::getVolumeType
 594 (
 595     const pointField& points,
 596     List<volumeType>& volType
 597 ) const
 598 {
 599     FatalErrorIn
 600     (
 601         "searchableSurfaceCollection::getVolumeType(const pointField&"
 602         ", List<volumeType>&) const"
 603     )   << "Volume type not supported for collection."
 604         << exit(FatalError);
 605 }
 606
 607
 608 void Foam::searchableSurfaceCollection::distribute
 609 (
 610     const List<treeBoundBox>& bbs,
 611     const bool keepNonLocal,
 612     autoPtr<mapDistribute>& faceMap,
 613     autoPtr<mapDistribute>& pointMap
 614 )
 615 {
 616     forAll(subGeom_, surfI)
 617     {
 618         // Note:Tranform the bounding boxes? Something like
 619         // pointField bbPoints =
 620         // cmptDivide
 621         // (
 622         //     transform_[surfI].localPosition
 623         //     (
 624         //         bbs[i].points()
 625         //     ),
 626         //     scale_[surfI]
 627         // );
 628         // treeBoundBox newBb(bbPoints);
 629
 630         // Note: what to do with faceMap, pointMap from multiple surfaces?
 631         subGeom_[surfI].distribute
 632         (
 633             bbs,
 634             keepNonLocal,
 635             faceMap,
 636             pointMap
 637         );
 638     }
 639 }
 640
 641
 642 void Foam::searchableSurfaceCollection::setField(const labelList& values)
 643 {
 644     forAll(subGeom_, surfI)
 645     {
 646         subGeom_[surfI].setField
 647         (
 648             static_cast<const labelList&>
 649             (
 650                 SubList<label>
 651                 (
 652                     values,
 653                     subGeom_[surfI].size(),
 654                     indexOffset_[surfI]
 655                 )
 656             )
 657         );
 658     }
 659 }
 660
 661
 662 void Foam::searchableSurfaceCollection::getField
 663 (
 664     const List<pointIndexHit>& info,
 665     labelList& values
 666 ) const
 667 {
 668     if (subGeom_.size() == 0)
 669     {}
 670     else if (subGeom_.size() == 1)
 671     {
 672         subGeom_[0].getField(info, values);
 673     }
 674     else
 675     {
 676         // Multiple surfaces. Sort by surface.
 677
 678         // Per surface the hit
 679         List<List<pointIndexHit> > surfInfo;
 680         // Per surface the original position
 681         List<List<label> > infoMap;
 682         sortHits(info, surfInfo, infoMap);
 683
 684         // Do surface tests
 685         forAll(surfInfo, surfI)
 686         {
 687             labelList surfValues;
 688             subGeom_[surfI].getField(surfInfo[surfI], surfValues);
 689
 690             if (surfValues.size())
 691             {
 692                 // Size values only when we have a surface that supports it.
 693                 values.setSize(info.size());
 694
 695                 const labelList& map = infoMap[surfI];
 696                 forAll(map, i)
 697                 {
 698                     values[map[i]] = surfValues[i];
 699                 }
 700             }
 701         }
 702     }
 703 }
 704
 705
 706 // ************************************************************************* //