src/dynamicMesh/meshCut/cellLooper/hexCellLooper.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 "hexCellLooper.H"
  27 #include "cellFeatures.H"
  28 #include "polyMesh.H"
  29 #include "cellModeller.H"
  30 #include "plane.H"
  31 #include "ListOps.H"
  32 #include "meshTools.H"
  33 #include "OFstream.H"
  34
  35 #include "addToRunTimeSelectionTable.H"
  36
  37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41 defineTypeNameAndDebug(hexCellLooper, 0);
  42 addToRunTimeSelectionTable(cellLooper, hexCellLooper, word);
  43 }
  44
  45
  46 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  47
  48 // Starting from cut edge start walking.
  49 bool Foam::hexCellLooper::walkHex
  50 (
  51     const label cellI,
  52     const label startFaceI,
  53     const label startEdgeI,
  54
  55     labelList& loop,
  56     scalarField& loopWeights
  57 ) const
  58 {
  59     label faceI = startFaceI;
  60
  61     label edgeI = startEdgeI;
  62
  63     label cutI = 0;
  64
  65     do
  66     {
  67         if (debug & 2)
  68         {
  69             Pout<< "    walkHex : inserting cut onto edge:" << edgeI
  70                 << " vertices:" << mesh().edges()[edgeI] << endl;
  71         }
  72
  73         // Store cut through edge. For now cut edges halfway.
  74         loop[cutI] = edgeToEVert(edgeI);
  75         loopWeights[cutI] = 0.5;
  76         cutI++;
  77
  78         faceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);
  79
  80         const edge& e = mesh().edges()[edgeI];
  81
  82         // Walk two edges further
  83         edgeI = meshTools::walkFace(mesh(), faceI, edgeI, e.end(), 2);
  84
  85         if (edgeI == startEdgeI)
  86         {
  87             break;
  88         }
  89     }
  90     while (true);
  91
  92     // Checks.
  93     if (cutI > 4)
  94     {
  95         Pout<< "hexCellLooper::walkHex" << "Problem : cell:" << cellI
  96             << " collected loop:";
  97         writeCuts(Pout, loop, loopWeights);
  98         Pout<< "loopWeights:" << loopWeights << endl;
  99
 100         return false;
 101     }
 102     else
 103     {
 104         return true;
 105     }
 106 }
 107
 108
 109
 110 void Foam::hexCellLooper::makeFace
 111 (
 112     const labelList& loop,
 113     const scalarField& loopWeights,
 114
 115     labelList& faceVerts,
 116     pointField& facePoints
 117 ) const
 118 {
 119     facePoints.setSize(loop.size());
 120     faceVerts.setSize(loop.size());
 121
 122     forAll(loop, cutI)
 123     {
 124         label cut = loop[cutI];
 125
 126         if (isEdge(cut))
 127         {
 128             label edgeI = getEdge(cut);
 129
 130             const edge& e = mesh().edges()[edgeI];
 131
 132             const point& v0 = mesh().points()[e.start()];
 133             const point& v1 = mesh().points()[e.end()];
 134
 135             facePoints[cutI] =
 136                 loopWeights[cutI]*v1 + (1.0-loopWeights[cutI])*v0;
 137         }
 138         else
 139         {
 140             label vertI = getVertex(cut);
 141
 142             facePoints[cutI] = mesh().points()[vertI];
 143         }
 144
 145         faceVerts[cutI] = cutI;
 146     }
 147 }
 148
 149
 150 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 151
 152 // Construct from components
 153 Foam::hexCellLooper::hexCellLooper(const polyMesh& mesh)
 154 :
 155     geomCellLooper(mesh),
 156     hex_(*(cellModeller::lookup("hex")))
 157 {}
 158
 159
 160 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 161
 162 Foam::hexCellLooper::~hexCellLooper()
 163 {}
 164
 165
 166 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 167
 168 bool Foam::hexCellLooper::cut
 169 (
 170     const vector& refDir,
 171     const label cellI,
 172     const boolList& vertIsCut,
 173     const boolList& edgeIsCut,
 174     const scalarField& edgeWeight,
 175
 176     labelList& loop,
 177     scalarField& loopWeights
 178 ) const
 179 {
 180     bool success = false;
 181
 182     if (mesh().cellShapes()[cellI].model() == hex_)
 183     {
 184         // Get starting edge. Note: should be compatible with way refDir is
 185         // determined.
 186         label edgeI = meshTools::cutDirToEdge(mesh(), cellI, refDir);
 187
 188         // Get any face using edge
 189         label face0;
 190         label face1;
 191         meshTools::getEdgeFaces(mesh(), cellI, edgeI, face0, face1);
 192
 193         // Walk circumference of hex, cutting edges only
 194         loop.setSize(4);
 195         loopWeights.setSize(4);
 196
 197         success = walkHex(cellI, face0, edgeI, loop, loopWeights);
 198     }
 199     else
 200     {
 201         success = geomCellLooper::cut
 202         (
 203             refDir,
 204             cellI,
 205             vertIsCut,
 206             edgeIsCut,
 207             edgeWeight,
 208
 209             loop,
 210             loopWeights
 211         );
 212     }
 213
 214     if (debug)
 215     {
 216         if (loop.empty())
 217         {
 218             WarningIn("hexCellLooper")
 219                 << "could not cut cell " << cellI << endl;
 220
 221             fileName cutsFile("hexCellLooper_" + name(cellI) + ".obj");
 222
 223             Pout<< "hexCellLooper : writing cell to " << cutsFile << endl;
 224
 225             OFstream cutsStream(cutsFile);
 226
 227             meshTools::writeOBJ
 228             (
 229                 cutsStream,
 230                 mesh().cells(),
 231                 mesh().faces(),
 232                 mesh().points(),
 233                 labelList(1, cellI)
 234             );
 235
 236             return false;
 237         }
 238
 239         // Check for duplicate cuts.
 240         labelHashSet loopSet(loop.size());
 241
 242         forAll(loop, elemI)
 243         {
 244             label elem = loop[elemI];
 245
 246             if (loopSet.found(elem))
 247             {
 248                 FatalErrorIn("hexCellLooper::walkHex") << " duplicate cut"
 249                     << abort(FatalError);
 250             }
 251             loopSet.insert(elem);
 252         }
 253
 254
 255         face faceVerts(loop.size());
 256         pointField facePoints(loop.size());
 257
 258         makeFace(loop, loopWeights, faceVerts, facePoints);
 259
 260         if ((faceVerts.mag(facePoints) < SMALL) || (loop.size() < 3))
 261         {
 262             FatalErrorIn("hexCellLooper::walkHex") << "Face:" << faceVerts
 263                 << " on points:" << facePoints << endl
 264                 << UIndirectList<point>(facePoints, faceVerts)()
 265                 << abort(FatalError);
 266         }
 267     }
 268     return success;
 269 }
 270
 271
 272 // No shortcuts for cutting with arbitrary plane.
 273 bool Foam::hexCellLooper::cut
 274 (
 275     const plane& cutPlane,
 276     const label cellI,
 277     const boolList& vertIsCut,
 278     const boolList& edgeIsCut,
 279     const scalarField& edgeWeight,
 280
 281     labelList& loop,
 282     scalarField& loopWeights
 283 ) const
 284 {
 285     return
 286         geomCellLooper::cut
 287         (
 288             cutPlane,
 289             cellI,
 290             vertIsCut,
 291             edgeIsCut,
 292             edgeWeight,
 293
 294             loop,
 295             loopWeights
 296         );
 297 }
 298
 299
 300 // ************************************************************************* //