src/OpenFOAM/primitives/globalIndexAndTransform/globalIndexAndTransform.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2010-2011 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 \*---------------------------------------------------------------------------*/
  25
  26 #include "globalIndexAndTransform.H"
  27 #include "cyclicPolyPatch.H"
  28
  29 // * * * * * * * * * * * * Private Static Data Members * * * * * * * * * * * //
  30
  31 defineTypeNameAndDebug(Foam::globalIndexAndTransform, 0);
  32
  33 const Foam::label Foam::globalIndexAndTransform::base_ = 32;
  34
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 Foam::label Foam::globalIndexAndTransform::matchTransform
  39 (
  40     const List<vectorTensorTransform>& refTransforms,
  41     label& matchedRefTransformI,
  42     const vectorTensorTransform& testTransform,
  43     scalar tolerance,
  44     bool checkBothSigns
  45 ) const
  46 {
  47     matchedRefTransformI = -1;
  48
  49     forAll(refTransforms, i)
  50     {
  51         const vectorTensorTransform& refTransform = refTransforms[i];
  52
  53         scalar maxVectorMag = sqrt
  54         (
  55             max(magSqr(testTransform.t()), magSqr(refTransform.t()))
  56         );
  57
  58         // Test the difference between vector parts to see if it is
  59         // less than tolerance times the larger vector part magnitude.
  60
  61         scalar vectorDiff =
  62             mag(refTransform.t() - testTransform.t())
  63            /(maxVectorMag + VSMALL)
  64            /tolerance;
  65
  66         // Test the difference between tensor parts to see if it is
  67         // less than the tolerance.  sqrt(3.0) factor used to scale
  68         // differnces as this is magnitude of a rotation tensor.  If
  69         // neither transform has a rotation, then the test is not
  70         // necessary.
  71
  72         scalar tensorDiff = 0;
  73
  74         if (refTransform.hasR() || testTransform.hasR())
  75         {
  76             tensorDiff =
  77                 mag(refTransform.R() - testTransform.R())
  78                /sqrt(3.0)
  79                /tolerance;
  80         }
  81
  82         // ...Diff result is < 1 if the test part matches the ref part
  83         // within tolerance
  84
  85         if (vectorDiff < 1 && tensorDiff < 1)
  86         {
  87             matchedRefTransformI = i;
  88
  89             return +1;
  90         }
  91
  92         if (checkBothSigns)
  93         {
  94             // Test the inverse transform differences too
  95
  96             vectorDiff =
  97                 mag(refTransform.t() + testTransform.t())
  98                /(maxVectorMag + VSMALL)
  99                /tolerance;
 100
 101             tensorDiff = 0;
 102
 103             if (refTransform.hasR() || testTransform.hasR())
 104             {
 105                 tensorDiff =
 106                     mag(refTransform.R() - testTransform.R().T())
 107                    /sqrt(3.0)
 108                    /tolerance;
 109             }
 110
 111             if (vectorDiff < 1 && tensorDiff < 1)
 112             {
 113                 matchedRefTransformI = i;
 114
 115                 return -1;
 116             }
 117         }
 118     }
 119
 120     return 0;
 121 }
 122
 123
 124 void Foam::globalIndexAndTransform::determineTransforms()
 125 {
 126     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 127
 128     transforms_ = List<vectorTensorTransform>(6);
 129     scalarField maxTol(6);
 130
 131     label nextTrans = 0;
 132
 133     label dummyMatch = -1;
 134
 135     forAll(patches, patchI)
 136     {
 137         const polyPatch& pp = patches[patchI];
 138
 139         // Note: special check for unordered cyclics. These are in fact
 140         // transform bcs and should probably be split off.
 141         if
 142         (
 143             isA<coupledPolyPatch>(pp)
 144         && !(
 145                 isA<cyclicPolyPatch>(pp)
 146              && (
 147                     refCast<const cyclicPolyPatch>(pp).transform()
 148                  == cyclicPolyPatch::NOORDERING
 149                 )
 150             )
 151         )
 152         {
 153             const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
 154
 155             if (cpp.separated())
 156             {
 157                 const vectorField& sepVecs = cpp.separation();
 158
 159                 forAll(sepVecs, sVI)
 160                 {
 161                     const vector& sepVec = sepVecs[sVI];
 162
 163                     if (mag(sepVec) > SMALL)
 164                     {
 165                         vectorTensorTransform transform(sepVec);
 166
 167                         if
 168                         (
 169                             matchTransform
 170                             (
 171                                 transforms_,
 172                                 dummyMatch,
 173                                 transform,
 174                                 cpp.matchTolerance(),
 175                                 false
 176                             ) == 0
 177                         )
 178                         {
 179                             if (nextTrans == 6)
 180                             {
 181                                 FatalErrorIn
 182                                 (
 183                                      "void Foam::globalIndexAndTransform::"
 184                                      "determineTransforms()"
 185                                 )   << "More than six unsigned transforms"
 186                                     << " detected:" << nl << transforms_
 187                                     << exit(FatalError);
 188                             }
 189                             transforms_[nextTrans] = transform;
 190                             maxTol[nextTrans++] = cpp.matchTolerance();
 191                         }
 192                     }
 193                 }
 194             }
 195             else if (!cpp.parallel())
 196             {
 197                 const tensorField& transTensors = cpp.reverseT();
 198
 199                 forAll(transTensors, tTI)
 200                 {
 201                     const tensor& transT = transTensors[tTI];
 202
 203                     if (mag(transT - I) > SMALL)
 204                     {
 205                         vectorTensorTransform transform(transT);
 206
 207                         if
 208                         (
 209                             matchTransform
 210                             (
 211                                 transforms_,
 212                                 dummyMatch,
 213                                 transform,
 214                                 cpp.matchTolerance(),
 215                                 false
 216                             ) == 0
 217                         )
 218                         {
 219                             if (nextTrans == 6)
 220                             {
 221                                 FatalErrorIn
 222                                 (
 223                                     "void Foam::globalIndexAndTransform::"
 224                                     "determineTransforms()"
 225                                 )   << "More than six unsigned transforms"
 226                                     << " detected:" << nl << transforms_
 227                                     << exit(FatalError);
 228                             }
 229                             transforms_[nextTrans] = transform;
 230                             maxTol[nextTrans++] = cpp.matchTolerance();
 231                         }
 232                     }
 233                 }
 234             }
 235         }
 236     }
 237
 238
 239     // Collect transforms on master
 240
 241     List<List<vectorTensorTransform> > allTransforms(Pstream::nProcs());
 242     allTransforms[Pstream::myProcNo()] = transforms_;
 243     Pstream::gatherList(allTransforms);
 244
 245     // Collect matching tolerance on master
 246     List<scalarField> allTols(Pstream::nProcs());
 247     allTols[Pstream::myProcNo()] = maxTol;
 248     Pstream::gatherList(allTols);
 249
 250     if (Pstream::master())
 251     {
 252         transforms_ = List<vectorTensorTransform>(3);
 253
 254         label nextTrans = 0;
 255
 256         forAll(allTransforms, procI)
 257         {
 258             const List<vectorTensorTransform>& procTransVecs =
 259                 allTransforms[procI];
 260
 261             forAll(procTransVecs, pSVI)
 262             {
 263                 const vectorTensorTransform& transform = procTransVecs[pSVI];
 264
 265                 if (mag(transform.t()) > SMALL || transform.hasR())
 266                 {
 267                     if
 268                     (
 269                         matchTransform
 270                         (
 271                             transforms_,
 272                             dummyMatch,
 273                             transform,
 274                             allTols[procI][pSVI],
 275                             true
 276                         ) ==  0
 277                     )
 278                     {
 279                         transforms_[nextTrans++] = transform;
 280                     }
 281
 282                     if (nextTrans > 3)
 283                     {
 284                         FatalErrorIn
 285                         (
 286                             "void Foam::globalIndexAndTransform::"
 287                             "determineTransforms()"
 288                         )
 289                             << "More than three independent basic "
 290                             << "transforms detected:" << nl
 291                             << allTransforms
 292                             << transforms_
 293                             << exit(FatalError);
 294                     }
 295                 }
 296             }
 297         }
 298
 299         transforms_.setSize(nextTrans);
 300     }
 301
 302     Pstream::scatter(transforms_);
 303
 304     if (transforms_.size() > 3)
 305     {
 306         WarningIn
 307         (
 308             "void globalIndexAndTransform::determineTransforms()"
 309         )   << "More than three independent basic "
 310             << "transforms detected:" << nl
 311             << transforms_ << nl
 312             << "This is not a space filling tiling and will probably"
 313             << " give problems for e.g. lagrangian tracking or interpolation"
 314             << endl;
 315     }
 316 }
 317
 318
 319 void Foam::globalIndexAndTransform::determineTransformPermutations()
 320 {
 321     label nTransformPermutations = pow(3, transforms_.size());
 322
 323     transformPermutations_.setSize(nTransformPermutations);
 324
 325     forAll(transformPermutations_, tPI)
 326     {
 327         vectorTensorTransform transform;
 328
 329         label transformIndex = tPI;
 330
 331         // Invert the ternary index encoding using repeated division by
 332         // three
 333
 334         forAll(transforms_, b)
 335         {
 336             const label w = (transformIndex % 3) - 1;
 337
 338             transformIndex /= 3;
 339
 340             if (w > 0)
 341             {
 342                 transform &= transforms_[b];
 343             }
 344             else if (w < 0)
 345             {
 346                 transform &= inv(transforms_[b]);
 347             }
 348         }
 349
 350         transformPermutations_[tPI] = transform;
 351     }
 352
 353
 354     // Encode index for 0 sign
 355     labelList permutationIndices(nIndependentTransforms(), 0);
 356     nullTransformIndex_ = encodeTransformIndex(permutationIndices);
 357 }
 358
 359
 360 void Foam::globalIndexAndTransform::determinePatchTransformSign()
 361 {
 362     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 363
 364     patchTransformSign_.setSize(patches.size(), Pair<label>(-1, 0));
 365
 366     label matchTransI = -1;
 367
 368     forAll(patches, patchI)
 369     {
 370         const polyPatch& pp = patches[patchI];
 371
 372         // Pout<< nl << patchI << " " << pp.name() << endl;
 373
 374         // Note: special check for unordered cyclics. These are in fact
 375         // transform bcs and should probably be split off.
 376         if
 377         (
 378             isA<coupledPolyPatch>(pp)
 379         && !(
 380                 isA<cyclicPolyPatch>(pp)
 381              && (
 382                     refCast<const cyclicPolyPatch>(pp).transform()
 383                  == cyclicPolyPatch::NOORDERING
 384                 )
 385             )
 386         )
 387         {
 388             const coupledPolyPatch& cpp =
 389             refCast<const coupledPolyPatch>(pp);
 390
 391             if (cpp.separated())
 392             {
 393                 const vectorField& sepVecs = cpp.separation();
 394
 395                 // Pout<< "sepVecs " << sepVecs << endl;
 396
 397                 // This loop is implicitly expecting only a single
 398                 // value for separation()
 399                 forAll(sepVecs, sVI)
 400                 {
 401                     const vector& sepVec = sepVecs[sVI];
 402
 403                     if (mag(sepVec) > SMALL)
 404                     {
 405                         vectorTensorTransform t(sepVec);
 406
 407                         label sign = matchTransform
 408                         (
 409                             transforms_,
 410                             matchTransI,
 411                             t,
 412                             cpp.matchTolerance(),
 413                             true
 414                         );
 415
 416                         // Pout<< sign << " " << matchTransI << endl;
 417
 418                         // List<label> permutation(transforms_.size(), 0);
 419
 420                         // permutation[matchTransI] = sign;
 421
 422                         // Pout<< encodeTransformIndex(permutation) << nl
 423                         //     << transformPermutations_
 424                         //        [
 425                         //            encodeTransformIndex(permutation)
 426                         //        ]
 427                         //     << endl;
 428
 429                         patchTransformSign_[patchI] =
 430                             Pair<label>(matchTransI, sign);
 431                     }
 432                 }
 433
 434             }
 435             else if (!cpp.parallel())
 436             {
 437                 const tensorField& transTensors = cpp.reverseT();
 438
 439                 // Pout<< "transTensors " << transTensors << endl;
 440
 441                 // This loop is implicitly expecting only a single
 442                 // value for reverseT()
 443                 forAll(transTensors, tTI)
 444                 {
 445                     const tensor& transT = transTensors[tTI];
 446
 447                     if (mag(transT - I) > SMALL)
 448                     {
 449                         vectorTensorTransform t(transT);
 450
 451                         label sign = matchTransform
 452                         (
 453                             transforms_,
 454                             matchTransI,
 455                             t,
 456                             cpp.matchTolerance(),
 457                             true
 458                         );
 459
 460                         // Pout<< sign << " " << matchTransI << endl;
 461
 462                         // List<label> permutation(transforms_.size(), 0);
 463
 464                         // permutation[matchTransI] = sign;
 465
 466                         // Pout<< encodeTransformIndex(permutation) << nl
 467                         //     << transformPermutations_
 468                         //        [
 469                         //            encodeTransformIndex(permutation)
 470                         //        ]
 471                         //     << endl;
 472
 473                         patchTransformSign_[patchI] =
 474                             Pair<label>(matchTransI, sign);
 475                     }
 476                 }
 477             }
 478         }
 479     }
 480
 481     // Pout<< patchTransformSign_ << endl;
 482 }
 483
 484
 485 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 486
 487 Foam::globalIndexAndTransform::globalIndexAndTransform
 488 (
 489     const polyMesh& mesh
 490 )
 491 :
 492     mesh_(mesh),
 493     transforms_(),
 494     transformPermutations_(),
 495     patchTransformSign_()
 496 {
 497     determineTransforms();
 498
 499     determineTransformPermutations();
 500
 501     determinePatchTransformSign();
 502
 503     if (debug && transforms_.size() > 0)
 504     {
 505         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 506
 507         Info<< "Determined global transforms :" << endl;
 508         Info<< "\t\ttranslation\trotation" << endl;
 509         forAll(transforms_, i)
 510         {
 511             Info<< '\t' << i << '\t';
 512             const vectorTensorTransform& trafo = transforms_[i];
 513             if (trafo.hasR())
 514             {
 515                  Info<< trafo.t() << '\t' << trafo.R();
 516             }
 517             else
 518             {
 519                  Info<< trafo.t() << '\t' << "---";
 520             }
 521             Info<< endl;
 522         }
 523         Info<< endl;
 524
 525
 526         Info<< "\tpatch\ttransform\tsign" << endl;
 527         forAll(patchTransformSign_, patchI)
 528         {
 529             if (patchTransformSign_[patchI].first() != -1)
 530             {
 531                 Info<< '\t' << patches[patchI].name()
 532                     << '\t' << patchTransformSign_[patchI].first()
 533                     << '\t' << patchTransformSign_[patchI].second()
 534                     << endl;
 535             }
 536         }
 537         Info<< endl;
 538
 539
 540         Info<< "Permutations of transformations:" << endl
 541             << "\t\ttranslation\trotation" << endl;
 542         forAll(transformPermutations_, i)
 543         {
 544             Info<< '\t' << i << '\t';
 545             const vectorTensorTransform& trafo = transformPermutations_[i];
 546             if (trafo.hasR())
 547             {
 548                  Info<< trafo.t() << '\t' << trafo.R();
 549             }
 550             else
 551             {
 552                  Info<< trafo.t() << '\t' << "---";
 553             }
 554             Info<< endl;
 555         }
 556         Info<< "nullTransformIndex:" << nullTransformIndex() << endl
 557             << endl;
 558     }
 559 }
 560
 561
 562 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 563
 564 Foam::globalIndexAndTransform::~globalIndexAndTransform()
 565 {}
 566
 567
 568 // ************************************************************************* //