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