src/meshTools/searchableSurface/searchableSurfaceCollection.C

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