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