src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.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 "coupledPolyPatch.H"
  27 #include "ListOps.H"
  28 #include "transform.H"
  29 #include "OFstream.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35     defineTypeNameAndDebug(coupledPolyPatch, 0);
  36
  37     const scalar coupledPolyPatch::defaultMatchTol_ = 1E-4;
  38
  39     template<>
  40     const char* NamedEnum<coupledPolyPatch::transformType, 4>::names[] =
  41     {
  42         "unknown",
  43         "rotational",
  44         "translational",
  45         "noOrdering"
  46     };
  47
  48     const NamedEnum<coupledPolyPatch::transformType, 4>
  49         coupledPolyPatch::transformTypeNames;
  50 }
  51
  52
  53 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  54
  55 void Foam::coupledPolyPatch::writeOBJ(Ostream& os, const point& pt)
  56 {
  57     os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
  58 }
  59
  60
  61 void Foam::coupledPolyPatch::writeOBJ
  62 (
  63     Ostream& os,
  64     const pointField& points,
  65     const labelList& pointLabels
  66 )
  67 {
  68     forAll(pointLabels, i)
  69     {
  70         writeOBJ(os, points[pointLabels[i]]);
  71     }
  72 }
  73
  74
  75 void Foam::coupledPolyPatch::writeOBJ
  76 (
  77     Ostream& os,
  78     const point& p0,
  79     const point& p1,
  80     label& vertI
  81 )
  82 {
  83     writeOBJ(os, p0);
  84     vertI++;
  85
  86     writeOBJ(os, p1);
  87     vertI++;
  88
  89     os  << "l " << vertI-1 << ' ' << vertI << nl;
  90 }
  91
  92
  93 void Foam::coupledPolyPatch::writeOBJ
  94 (
  95     const fileName& fName,
  96     const UList<face>& faces,
  97     const pointField& points
  98 )
  99 {
 100     OFstream os(fName);
 101
 102     Map<label> foamToObj(4*faces.size());
 103
 104     label vertI = 0;
 105
 106     forAll(faces, i)
 107     {
 108         const face& f = faces[i];
 109
 110         forAll(f, fp)
 111         {
 112             if (foamToObj.insert(f[fp], vertI))
 113             {
 114                 writeOBJ(os, points[f[fp]]);
 115                 vertI++;
 116             }
 117         }
 118
 119         os << 'l';
 120         forAll(f, fp)
 121         {
 122             os << ' ' << foamToObj[f[fp]]+1;
 123         }
 124         os << ' ' << foamToObj[f[0]]+1 << nl;
 125     }
 126 }
 127
 128
 129 Foam::pointField Foam::coupledPolyPatch::getAnchorPoints
 130 (
 131     const UList<face>& faces,
 132     const pointField& points
 133 )
 134 {
 135     pointField anchors(faces.size());
 136
 137     forAll(faces, faceI)
 138     {
 139         anchors[faceI] = points[faces[faceI][0]];
 140     }
 141
 142     return anchors;
 143 }
 144
 145
 146 Foam::scalarField Foam::coupledPolyPatch::calcFaceTol
 147 (
 148     const scalar matchTol,
 149     const UList<face>& faces,
 150     const pointField& points,
 151     const pointField& faceCentres
 152 )
 153 {
 154     // Calculate typical distance per face
 155     scalarField tols(faces.size());
 156
 157     forAll(faces, faceI)
 158     {
 159         const point& cc = faceCentres[faceI];
 160
 161         const face& f = faces[faceI];
 162
 163         // 1. calculate a typical size of the face. Use maximum distance
 164         //    to face centre
 165         scalar maxLenSqr = -GREAT;
 166         // 2. as measure of truncation error when comparing two coordinates
 167         //    use SMALL * maximum component
 168         scalar maxCmpt = -GREAT;
 169
 170         forAll(f, fp)
 171         {
 172             const point& pt = points[f[fp]];
 173             maxLenSqr = max(maxLenSqr, magSqr(pt - cc));
 174             maxCmpt = max(maxCmpt, cmptMax(cmptMag(pt)));
 175         }
 176         tols[faceI] = max(SMALL*maxCmpt, matchTol*Foam::sqrt(maxLenSqr));
 177     }
 178     return tols;
 179 }
 180
 181
 182 Foam::label Foam::coupledPolyPatch::getRotation
 183 (
 184     const pointField& points,
 185     const face& f,
 186     const point& anchor,
 187     const scalar tol
 188 )
 189 {
 190     label anchorFp = -1;
 191     scalar minDistSqr = GREAT;
 192
 193     forAll(f, fp)
 194     {
 195         scalar distSqr = magSqr(anchor - points[f[fp]]);
 196
 197         if (distSqr < minDistSqr)
 198         {
 199             minDistSqr = distSqr;
 200             anchorFp = fp;
 201         }
 202     }
 203
 204     if (anchorFp == -1 || mag(minDistSqr) > tol)
 205     {
 206         return -1;
 207     }
 208     else
 209     {
 210         // Positive rotation
 211         return (f.size() - anchorFp) % f.size();
 212     }
 213 }
 214
 215
 216 void Foam::coupledPolyPatch::calcTransformTensors
 217 (
 218     const vectorField& Cf,
 219     const vectorField& Cr,
 220     const vectorField& nf,
 221     const vectorField& nr,
 222     const scalarField& smallDist,
 223     const scalar absTol,
 224     const transformType transform
 225 ) const
 226 {
 227     if (debug)
 228     {
 229         Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
 230             << "    transform:" << transformTypeNames[transform] << nl
 231             << "    (half)size:" << Cf.size() << nl
 232             << "    absTol:" << absTol << nl
 233             << "    smallDist min:" << min(smallDist) << nl
 234             << "    smallDist max:" << max(smallDist) << nl
 235             << "    sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl;
 236     }
 237
 238     // Tolerance calculation.
 239     // - normal calculation: assume absTol is the absolute error in a
 240     // single normal/transformation calculation. Consists both of numerical
 241     // precision (on the order of SMALL and of writing precision
 242     // (from e.g. decomposition)
 243     // Then the overall error of summing the normals is sqrt(size())*absTol
 244     // - separation calculation: pass in from the outside an allowable error.
 245
 246     if (Cf.size() == 0)
 247     {
 248         // Dummy geometry.
 249         separation_.setSize(0);
 250         forwardT_ = I;
 251         reverseT_ = I;
 252         collocated_.setSize(0);
 253     }
 254     else
 255     {
 256         scalar error = absTol*Foam::sqrt(1.0*Cf.size());
 257
 258         if (debug)
 259         {
 260             Pout<< "    error:" << error << endl;
 261         }
 262
 263         if
 264         (
 265             transform == ROTATIONAL
 266          || (
 267                 transform != TRANSLATIONAL
 268              && (sum(mag(nf & nr)) < Cf.size()-error)
 269             )
 270         )
 271         {
 272             // Type is rotation or unknown and normals not aligned
 273
 274             // Assume per-face differing transformation, correct later
 275
 276             separation_.setSize(0);
 277
 278             forwardT_.setSize(Cf.size());
 279             reverseT_.setSize(Cf.size());
 280             collocated_.setSize(Cf.size());
 281             collocated_ = false;
 282
 283             forAll(forwardT_, facei)
 284             {
 285                 forwardT_[facei] = rotationTensor(-nr[facei], nf[facei]);
 286                 reverseT_[facei] = rotationTensor(nf[facei], -nr[facei]);
 287             }
 288
 289             if (debug)
 290             {
 291                 Pout<< "    sum(mag(forwardT_ - forwardT_[0])):"
 292                     << sum(mag(forwardT_ - forwardT_[0]))
 293                     << endl;
 294             }
 295
 296             if (sum(mag(forwardT_ - forwardT_[0])) < error)
 297             {
 298                 forwardT_.setSize(1);
 299                 reverseT_.setSize(1);
 300                 collocated_.setSize(1);
 301
 302                 if (debug)
 303                 {
 304                     Pout<< "    difference in rotation less than"
 305                         << " local tolerance "
 306                         << error << ". Assuming uniform rotation." << endl;
 307                 }
 308             }
 309         }
 310         else
 311         {
 312             // Translational or (unknown and normals aligned)
 313
 314             forwardT_.setSize(0);
 315             reverseT_.setSize(0);
 316
 317             separation_ = Cr - Cf;
 318
 319             collocated_.setSize(separation_.size());
 320
 321             // Three situations:
 322             // - separation is zero. No separation.
 323             // - separation is same. Single separation vector.
 324             // - separation differs per face. Separation vectorField.
 325
 326             // Check for different separation per face
 327             bool sameSeparation = true;
 328             bool doneWarning = false;
 329
 330             forAll(separation_, facei)
 331             {
 332                 scalar smallSqr = sqr(smallDist[facei]);
 333
 334                 collocated_[facei] = (magSqr(separation_[facei]) < smallSqr);
 335
 336                 // Check if separation differing w.r.t. face 0.
 337                 if (magSqr(separation_[facei] - separation_[0]) > smallSqr)
 338                 {
 339                     sameSeparation = false;
 340
 341                     if (!doneWarning && debug)
 342                     {
 343                         doneWarning = true;
 344
 345                         Pout<< "    separation " << separation_[facei]
 346                             << " at " << facei
 347                             << " differs from separation[0] " << separation_[0]
 348                             << " by more than local tolerance "
 349                             << smallDist[facei]
 350                             << ". Assuming non-uniform separation." << endl;
 351                     }
 352                 }
 353             }
 354
 355             if (sameSeparation)
 356             {
 357                 // Check for zero separation (at 0 so everywhere)
 358                 if (collocated_[0])
 359                 {
 360                     if (debug)
 361                     {
 362                         Pout<< "    separation " << mag(separation_[0])
 363                             << " less than local tolerance " << smallDist[0]
 364                             << ". Assuming zero separation." << endl;
 365                     }
 366
 367                     separation_.setSize(0);
 368                     collocated_ = boolList(1, true);
 369                 }
 370                 else
 371                 {
 372                     if (debug)
 373                     {
 374                         Pout<< "    separation " << mag(separation_[0])
 375                             << " more than local tolerance " << smallDist[0]
 376                             << ". Assuming uniform separation." << endl;
 377                     }
 378
 379                     separation_.setSize(1);
 380                     collocated_ = boolList(1, false);
 381                 }
 382             }
 383         }
 384     }
 385
 386     if (debug)
 387     {
 388         Pout<< "    separation_:" << separation_.size() << nl
 389             << "    forwardT size:" << forwardT_.size() << endl;
 390     }
 391 }
 392
 393
 394 // * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * * * * //
 395
 396 Foam::coupledPolyPatch::coupledPolyPatch
 397 (
 398     const word& name,
 399     const label size,
 400     const label start,
 401     const label index,
 402     const polyBoundaryMesh& bm
 403 )
 404 :
 405     polyPatch(name, size, start, index, bm),
 406     matchTolerance_(defaultMatchTol_)
 407 {}
 408
 409
 410 Foam::coupledPolyPatch::coupledPolyPatch
 411 (
 412     const word& name,
 413     const dictionary& dict,
 414     const label index,
 415     const polyBoundaryMesh& bm
 416 )
 417 :
 418     polyPatch(name, dict, index, bm),
 419     matchTolerance_(dict.lookupOrDefault("matchTolerance", defaultMatchTol_))
 420 {}
 421
 422
 423 Foam::coupledPolyPatch::coupledPolyPatch
 424 (
 425     const coupledPolyPatch& pp,
 426     const polyBoundaryMesh& bm
 427 )
 428 :
 429     polyPatch(pp, bm),
 430     matchTolerance_(pp.matchTolerance_)
 431 {}
 432
 433
 434 Foam::coupledPolyPatch::coupledPolyPatch
 435 (
 436     const coupledPolyPatch& pp,
 437     const polyBoundaryMesh& bm,
 438     const label index,
 439     const label newSize,
 440     const label newStart
 441 )
 442 :
 443     polyPatch(pp, bm, index, newSize, newStart),
 444     matchTolerance_(pp.matchTolerance_)
 445 {}
 446
 447
 448 Foam::coupledPolyPatch::coupledPolyPatch
 449 (
 450     const coupledPolyPatch& pp,
 451     const polyBoundaryMesh& bm,
 452     const label index,
 453     const labelUList& mapAddressing,
 454     const label newStart
 455 )
 456 :
 457     polyPatch(pp, bm, index, mapAddressing, newStart),
 458     matchTolerance_(pp.matchTolerance_)
 459 {}
 460
 461
 462 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 463
 464 Foam::coupledPolyPatch::~coupledPolyPatch()
 465 {}
 466
 467
 468 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 469
 470 void Foam::coupledPolyPatch::write(Ostream& os) const
 471 {
 472     polyPatch::write(os);
 473     //if (matchTolerance_ != defaultMatchTol_)
 474     {
 475         os.writeKeyword("matchTolerance") << matchTolerance_
 476             << token::END_STATEMENT << nl;
 477     }
 478 }
 479
 480
 481 // ************************************************************************* //