src/dynamicMesh/meshCut/directions/directions.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 "directions.H"
  28 #include "polyMesh.H"
  29 #include "twoDPointCorrector.H"
  30 #include "directionInfo.H"
  31 #include "MeshWave.H"
  32 #include "OFstream.H"
  33 #include "meshTools.H"
  34 #include "hexMatcher.H"
  35 #include "Switch.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 template<>
  40 const char* Foam::NamedEnum<Foam::directions::directionType, 3>::names[] =
  41 {
  42     "tan1",
  43     "tan2",
  44     "normal"
  45 };
  46
  47 const Foam::NamedEnum<Foam::directions::directionType, 3>
  48     Foam::directions::directionTypeNames_;
  49
  50
  51
  52 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  53
  54 // For debugging
  55 void Foam::directions::writeOBJ(Ostream& os, const point& pt)
  56 {
  57     os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
  58 }
  59
  60 // For debugging
  61 void Foam::directions::writeOBJ
  62 (
  63     Ostream& os,
  64     const point& pt0,
  65     const point& pt1,
  66     label& vertI
  67 )
  68 {
  69     writeOBJ(os, pt0);
  70     writeOBJ(os, pt1);
  71
  72     os << "l " << vertI + 1 << ' ' << vertI + 2 << endl;
  73
  74     vertI += 2;
  75 }
  76
  77
  78 // Dump to file.
  79 void Foam::directions::writeOBJ
  80 (
  81     const fileName& fName,
  82     const primitiveMesh& mesh,
  83     const vectorField& dirs
  84 )
  85 {
  86     Pout<< "Writing cell info to " << fName << " as vectors at the cellCentres"
  87         << endl << endl;
  88
  89     OFstream xDirStream(fName);
  90
  91     label vertI = 0;
  92
  93     forAll(dirs, cellI)
  94     {
  95         const point& ctr = mesh.cellCentres()[cellI];
  96
  97         // Calculate local length scale
  98         scalar minDist = GREAT;
  99
 100         const labelList& nbrs = mesh.cellCells()[cellI];
 101
 102         forAll(nbrs, nbrI)
 103         {
 104             minDist = min(minDist, mag(mesh.cellCentres()[nbrs[nbrI]] - ctr));
 105         }
 106
 107         scalar scale = 0.5*minDist;
 108
 109         writeOBJ(xDirStream, ctr, ctr + scale*dirs[cellI], vertI);
 110     }
 111 }
 112
 113
 114 void Foam::directions::check2D
 115 (
 116     const twoDPointCorrector* correct2DPtr,
 117     const vector& vec
 118 )
 119 {
 120     if (correct2DPtr)
 121     {
 122         if (mag(correct2DPtr->planeNormal() & vec) > 1E-6)
 123         {
 124             FatalErrorIn("check2D") << "Specified vector " << vec
 125                 << "is not normal to plane defined in dynamicMeshDict."
 126                 << endl
 127                 << "Either make case 3D or adjust vector."
 128                 << exit(FatalError);
 129         }
 130     }
 131 }
 132
 133
 134 // Get direction on all cells
 135 Foam::vectorField Foam::directions::propagateDirection
 136 (
 137     const polyMesh& mesh,
 138     const bool useTopo,
 139     const polyPatch& pp,
 140     const vectorField& ppField,
 141     const vector& defaultDir
 142 )
 143 {
 144     // Seed all faces on patch
 145     labelList changedFaces(pp.size());
 146     List<directionInfo> changedFacesInfo(pp.size());
 147
 148     if (useTopo)
 149     {
 150         forAll(pp, patchFaceI)
 151         {
 152             label meshFaceI = pp.start() + patchFaceI;
 153
 154             label cellI = mesh.faceOwner()[meshFaceI];
 155
 156             if (!hexMatcher().isA(mesh, cellI))
 157             {
 158                 FatalErrorIn("propagateDirection")
 159                     << "useHexTopology specified but cell " << cellI
 160                     << " on face " << patchFaceI << " of patch " << pp.name()
 161                     << " is not a hex" << exit(FatalError);
 162             }
 163
 164             const vector& cutDir = ppField[patchFaceI];
 165
 166             // Get edge(bundle) on cell most in direction of cutdir
 167             label edgeI = meshTools::cutDirToEdge(mesh, cellI, cutDir);
 168
 169             // Convert edge into index on face
 170             label faceIndex =
 171                 directionInfo::edgeToFaceIndex
 172                 (
 173                     mesh,
 174                     cellI,
 175                     meshFaceI,
 176                     edgeI
 177                 );
 178
 179             // Set initial face and direction
 180             changedFaces[patchFaceI] = meshFaceI;
 181             changedFacesInfo[patchFaceI] =
 182                 directionInfo
 183                 (
 184                     faceIndex,
 185                     cutDir
 186                 );
 187         }
 188     }
 189     else
 190     {
 191         forAll(pp, patchFaceI)
 192         {
 193             changedFaces[patchFaceI] = pp.start() + patchFaceI;
 194             changedFacesInfo[patchFaceI] =
 195                 directionInfo
 196                 (
 197                     -2,         // Geometric information only
 198                     ppField[patchFaceI]
 199                 );
 200         }
 201     }
 202
 203     MeshWave<directionInfo> directionCalc
 204     (
 205         mesh,
 206         changedFaces,
 207         changedFacesInfo,
 208         mesh.nCells()
 209     );
 210
 211     const List<directionInfo>& cellInfo = directionCalc.allCellInfo();
 212
 213     vectorField dirField(cellInfo.size());
 214
 215     label nUnset = 0;
 216     label nGeom = 0;
 217     label nTopo = 0;
 218
 219     forAll(cellInfo, cellI)
 220     {
 221         label index = cellInfo[cellI].index();
 222
 223         if (index == -3)
 224         {
 225             // Never visited
 226             WarningIn("propagateDirection")
 227                 << "Cell " << cellI << " never visited to determine "
 228                 << "local coordinate system" << endl
 229                 << "Using direction " << defaultDir << " instead" << endl;
 230
 231             dirField[cellI] = defaultDir;
 232
 233             nUnset++;
 234         }
 235         else if (index == -2)
 236         {
 237             // Geometric direction
 238             dirField[cellI] = cellInfo[cellI].n();
 239
 240             nGeom++;
 241         }
 242         else if (index == -1)
 243         {
 244             FatalErrorIn("propagateDirection")
 245                 << "Illegal index " << index << endl
 246                 << "Value is only allowed on faces" << abort(FatalError);
 247         }
 248         else
 249         {
 250             // Topological edge cut. Convert into average cut direction.
 251             dirField[cellI] = meshTools::edgeToCutDir(mesh, cellI, index);
 252
 253             nTopo++;
 254         }
 255     }
 256
 257     Pout<< "Calculated local coords for " << defaultDir
 258         << endl
 259         << "    Geometric cut cells   : " << nGeom << endl
 260         << "    Topological cut cells : " << nTopo << endl
 261         << "    Unset cells           : " << nUnset << endl
 262         << endl;
 263
 264     return dirField;
 265 }
 266
 267
 268 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 269
 270 // Construct from dictionary
 271 Foam::directions::directions
 272 (
 273     const polyMesh& mesh,
 274     const dictionary& dict,
 275     const twoDPointCorrector* correct2DPtr
 276 )
 277 :
 278     List<vectorField>(wordList(dict.lookup("directions")).size())
 279 {
 280     const wordList wantedDirs(dict.lookup("directions"));
 281
 282     bool wantNormal = false;
 283     bool wantTan1 = false;
 284     bool wantTan2 = false;
 285
 286     forAll(wantedDirs, i)
 287     {
 288         directionType wantedDir = directionTypeNames_[wantedDirs[i]];
 289
 290         if (wantedDir == NORMAL)
 291         {
 292             wantNormal = true;
 293         }
 294         else if (wantedDir == TAN1)
 295         {
 296             wantTan1 = true;
 297         }
 298         else if (wantedDir == TAN2)
 299         {
 300             wantTan2 = true;
 301         }
 302     }
 303
 304
 305     label nDirs = 0;
 306
 307     word coordSystem(dict.lookup("coordinateSystem"));
 308
 309     if (coordSystem == "global")
 310     {
 311         const dictionary& globalDict = dict.subDict("globalCoeffs");
 312
 313         vector tan1(globalDict.lookup("tan1"));
 314         check2D(correct2DPtr, tan1);
 315
 316         vector tan2(globalDict.lookup("tan2"));
 317         check2D(correct2DPtr, tan2);
 318
 319         vector normal = tan1 ^ tan2;
 320         normal /= mag(normal);
 321
 322         Pout<< "Global Coordinate system:" << endl
 323             << "     normal : " << normal << endl
 324             << "     tan1   : " << tan1 << endl
 325             << "     tan2   : " << tan2
 326             << endl << endl;
 327
 328         if (wantNormal)
 329         {
 330             operator[](nDirs++) = vectorField(1, normal);
 331         }
 332         if (wantTan1)
 333         {
 334             operator[](nDirs++) = vectorField(1, tan1);
 335         }
 336         if (wantTan2)
 337         {
 338             operator[](nDirs++) = vectorField(1, tan2);
 339         }
 340     }
 341     else if (coordSystem == "patchLocal")
 342     {
 343         const dictionary& patchDict = dict.subDict("patchLocalCoeffs");
 344
 345         word patchName(patchDict.lookup("patch"));
 346
 347         label patchI = mesh.boundaryMesh().findPatchID(patchName);
 348
 349         if (patchI == -1)
 350         {
 351             FatalErrorIn
 352             (
 353                 "directions::directions(const polyMesh&, const dictionary&,"
 354                 "const twoDPointCorrector*"
 355             )   << "Cannot find patch "
 356                 << patchName << exit(FatalError);
 357         }
 358
 359         // Take zeroth face on patch
 360         const polyPatch& pp = mesh.boundaryMesh()[patchI];
 361
 362         vector tan1(patchDict.lookup("tan1"));
 363
 364         const vector& n0 = pp.faceNormals()[0];
 365
 366         if (correct2DPtr)
 367         {
 368             tan1 = correct2DPtr->planeNormal() ^ n0;
 369
 370             WarningIn
 371             (
 372                 "directions::directions(const polyMesh&, const dictionary&,"
 373                 "const twoDPointCorrector*"
 374             )   << "Discarding user specified tan1 since 2D case." << endl
 375                 << "Recalculated tan1 from face normal and planeNormal as "
 376                 << tan1 << endl << endl;
 377         }
 378
 379         Switch useTopo(dict.lookup("useHexTopology"));
 380
 381         vectorField normalDirs;
 382         vectorField tan1Dirs;
 383
 384         if (wantNormal || wantTan2)
 385         {
 386             normalDirs =
 387                 propagateDirection
 388                 (
 389                     mesh,
 390                     useTopo,
 391                     pp,
 392                     pp.faceNormals(),
 393                     n0
 394                 );
 395
 396             if (wantNormal)
 397             {
 398                 operator[](nDirs++) = normalDirs;
 399
 400                 //// Dump to file.
 401                 //writeOBJ("normal.obj", mesh, normalDirs);
 402             }
 403         }
 404
 405         if (wantTan1 || wantTan2)
 406         {
 407             tan1Dirs =
 408                 propagateDirection
 409                 (
 410                     mesh,
 411                     useTopo,
 412                     pp,
 413                     vectorField(pp.size(), tan1),
 414                     tan1
 415                 );
 416
 417
 418             if (wantTan1)
 419             {
 420                 operator[](nDirs++) = tan1Dirs;
 421
 422                 //// Dump to file.
 423                 //writeOBJ("tan1.obj", mesh, tan1Dirs);
 424             }
 425         }
 426         if (wantTan2)
 427         {
 428             vectorField tan2Dirs = normalDirs ^ tan1Dirs;
 429
 430             operator[](nDirs++) = tan2Dirs;
 431
 432             //// Dump to file.
 433             //writeOBJ("tan2.obj", mesh, tan2Dirs);
 434         }
 435     }
 436     else
 437     {
 438         FatalErrorIn
 439         (
 440             "directions::directions(const polyMesh&, const dictionary&,"
 441             "const twoDPointCorrector*"
 442         )   << "Unknown coordinate system "
 443             << coordSystem << endl
 444             << "Known types are global and patchLocal" << exit(FatalError);
 445     }
 446 }
 447
 448
 449 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 450
 451
 452 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 453
 454
 455 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 456
 457
 458 // * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
 459
 460
 461 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 462
 463
 464 // ************************************************************************* //