src/sampling/cuttingPlane/cuttingPlane.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 "cuttingPlane.H"
  28 #include "primitiveMesh.H"
  29 #include "linePointRef.H"
  30 #include "meshTools.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 // set values for what is close to zero and what is considered to
  35 // be positive (and not just rounding noise)
  36 //! @cond localScope
  37 const Foam::scalar zeroish  = Foam::SMALL;
  38 const Foam::scalar positive = Foam::SMALL * 1E3;
  39 //! @endcond localScope
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 // Find cut cells
  44 void Foam::cuttingPlane::calcCutCells
  45 (
  46     const primitiveMesh& mesh,
  47     const scalarField& dotProducts,
  48     const UList<label>& cellIdLabels
  49 )
  50 {
  51     const labelListList& cellEdges = mesh.cellEdges();
  52     const edgeList& edges = mesh.edges();
  53
  54     label listSize = cellEdges.size();
  55     if (&cellIdLabels)
  56     {
  57         listSize = cellIdLabels.size();
  58     }
  59
  60     cutCells_.setSize(listSize);
  61     label cutcellI(0);
  62
  63     // Find the cut cells by detecting any cell that uses points with
  64     // opposing dotProducts.
  65     for (label listI = 0; listI < listSize; ++listI)
  66     {
  67         label cellI = listI;
  68         if (&cellIdLabels)
  69         {
  70             cellI = cellIdLabels[listI];
  71         }
  72
  73         const labelList& cEdges = cellEdges[cellI];
  74
  75         label nCutEdges = 0;
  76
  77         forAll(cEdges, i)
  78         {
  79             const edge& e = edges[cEdges[i]];
  80
  81             if
  82             (
  83                 (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
  84              || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
  85             )
  86             {
  87                 nCutEdges++;
  88
  89                 if (nCutEdges > 2)
  90                 {
  91                     cutCells_[cutcellI++] = cellI;
  92
  93                     break;
  94                 }
  95             }
  96         }
  97     }
  98
  99     // Set correct list size
 100     cutCells_.setSize(cutcellI);
 101 }
 102
 103
 104 // Determine for each edge the intersection point. Calculates
 105 // - cutPoints_ : coordinates of all intersection points
 106 // - edgePoint  : per edge -1 or the index into cutPoints
 107 void Foam::cuttingPlane::intersectEdges
 108 (
 109     const primitiveMesh& mesh,
 110     const scalarField& dotProducts,
 111     List<label>& edgePoint
 112 )
 113 {
 114     // Use the dotProducts to find out the cut edges.
 115     const edgeList& edges = mesh.edges();
 116     const pointField& points = mesh.points();
 117
 118     // Per edge -1 or the label of the intersection point
 119     edgePoint.setSize(edges.size());
 120
 121     DynamicList<point> dynCuttingPoints(4*cutCells_.size());
 122
 123     forAll(edges, edgeI)
 124     {
 125         const edge& e = edges[edgeI];
 126
 127         if
 128         (
 129             (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
 130          || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
 131         )
 132         {
 133             // Edge is cut
 134             edgePoint[edgeI] = dynCuttingPoints.size();
 135
 136             const point& p0 = points[e[0]];
 137             const point& p1 = points[e[1]];
 138
 139             scalar alpha = lineIntersect(linePointRef(p0, p1));
 140
 141             if (alpha < zeroish)
 142             {
 143                 dynCuttingPoints.append(p0);
 144             }
 145             else if (alpha >= 1.0)
 146             {
 147                 dynCuttingPoints.append(p1);
 148             }
 149             else
 150             {
 151                 dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
 152             }
 153         }
 154         else
 155         {
 156             edgePoint[edgeI] = -1;
 157         }
 158     }
 159
 160     this->storedPoints().transfer(dynCuttingPoints);
 161 }
 162
 163
 164 // Coming from startEdgeI cross the edge to the other face
 165 // across to the next cut edge.
 166 bool Foam::cuttingPlane::walkCell
 167 (
 168     const primitiveMesh& mesh,
 169     const UList<label>& edgePoint,
 170     const label cellI,
 171     const label startEdgeI,
 172     DynamicList<label>& faceVerts
 173 )
 174 {
 175     label faceI = -1;
 176     label edgeI = startEdgeI;
 177
 178     label nIter = 0;
 179
 180     faceVerts.clear();
 181     do
 182     {
 183         faceVerts.append(edgePoint[edgeI]);
 184
 185         // Cross edge to other face
 186         faceI = meshTools::otherFace(mesh, cellI, faceI, edgeI);
 187
 188         // Find next cut edge on face.
 189         const labelList& fEdges = mesh.faceEdges()[faceI];
 190
 191         label nextEdgeI = -1;
 192
 193         //Note: here is where we should check for whether there are more
 194         // than 2 intersections with the face (warped/non-convex face).
 195         // If so should e.g. decompose the cells on both faces and redo
 196         // the calculation.
 197
 198         forAll(fEdges, i)
 199         {
 200             label edge2I = fEdges[i];
 201
 202             if (edge2I != edgeI && edgePoint[edge2I] != -1)
 203             {
 204                 nextEdgeI = edge2I;
 205                 break;
 206             }
 207         }
 208
 209         if (nextEdgeI == -1)
 210         {
 211             // Did not find another cut edge on faceI. Do what?
 212             WarningIn("Foam::cuttingPlane::walkCell")
 213                 << "Did not find closed walk along surface of cell " << cellI
 214                 << " starting from edge " << startEdgeI
 215                 << " in " << nIter << " iterations." << nl
 216                 << "Collected cutPoints so far:" << faceVerts
 217                 << endl;
 218
 219             return false;
 220         }
 221
 222         edgeI = nextEdgeI;
 223
 224         nIter++;
 225
 226         if (nIter > 1000)
 227         {
 228             WarningIn("Foam::cuttingPlane::walkCell")
 229                 << "Did not find closed walk along surface of cell " << cellI
 230                 << " starting from edge " << startEdgeI
 231                 << " in " << nIter << " iterations." << nl
 232                 << "Collected cutPoints so far:" << faceVerts
 233                 << endl;
 234             return false;
 235         }
 236
 237     } while (edgeI != startEdgeI);
 238
 239
 240     if (faceVerts.size() >= 3)
 241     {
 242         return true;
 243     }
 244     else
 245     {
 246         WarningIn("Foam::cuttingPlane::walkCell")
 247             << "Did not find closed walk along surface of cell " << cellI
 248             << " starting from edge " << startEdgeI << nl
 249             << "Collected cutPoints so far:" << faceVerts
 250             << endl;
 251
 252         return false;
 253     }
 254 }
 255
 256
 257 // For every cut cell determine a walk through all? its cuts.
 258 void Foam::cuttingPlane::walkCellCuts
 259 (
 260     const primitiveMesh& mesh,
 261     const UList<label>& edgePoint
 262 )
 263 {
 264     const pointField& cutPoints = this->points();
 265
 266     // use dynamic lists to handle triangulation and/or missed cuts
 267     DynamicList<face>  dynCutFaces(cutCells_.size());
 268     DynamicList<label> dynCutCells(cutCells_.size());
 269
 270     // scratch space for calculating the face vertices
 271     DynamicList<label> faceVerts(10);
 272
 273     forAll(cutCells_, i)
 274     {
 275         label cellI = cutCells_[i];
 276
 277         // Find the starting edge to walk from.
 278         const labelList& cEdges = mesh.cellEdges()[cellI];
 279
 280         label startEdgeI = -1;
 281
 282         forAll(cEdges, cEdgeI)
 283         {
 284             label edgeI = cEdges[cEdgeI];
 285
 286             if (edgePoint[edgeI] != -1)
 287             {
 288                 startEdgeI = edgeI;
 289                 break;
 290             }
 291         }
 292
 293         // Check for the unexpected ...
 294         if (startEdgeI == -1)
 295         {
 296             FatalErrorIn("Foam::cuttingPlane::walkCellCuts")
 297                 << "Cannot find cut edge for cut cell " << cellI
 298                 << abort(FatalError);
 299         }
 300
 301         // Walk from starting edge around the circumference of the cell.
 302         bool okCut = walkCell
 303         (
 304             mesh,
 305             edgePoint,
 306             cellI,
 307             startEdgeI,
 308             faceVerts
 309         );
 310
 311         if (okCut)
 312         {
 313             face f(faceVerts);
 314
 315             // Orient face to point in the same direction as the plane normal
 316             if ((f.normal(cutPoints) && normal()) < 0)
 317             {
 318                 f = f.reverseFace();
 319             }
 320
 321             // the cut faces are usually quite ugly, so always triangulate
 322             label nTri = f.triangles(cutPoints, dynCutFaces);
 323             while (nTri--)
 324             {
 325                 dynCutCells.append(cellI);
 326             }
 327         }
 328     }
 329
 330     this->storedFaces().transfer(dynCutFaces);
 331     cutCells_.transfer(dynCutCells);
 332 }
 333
 334
 335 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 336
 337 // Construct without cutting
 338 Foam::cuttingPlane::cuttingPlane(const plane& pln)
 339 :
 340     plane(pln)
 341 {}
 342
 343
 344 // Construct from plane and mesh reference, restricted to a list of cells
 345 Foam::cuttingPlane::cuttingPlane
 346 (
 347     const plane& pln,
 348     const primitiveMesh& mesh,
 349     const UList<label>& cellIdLabels
 350 )
 351 :
 352     plane(pln)
 353 {
 354     reCut(mesh, cellIdLabels);
 355 }
 356
 357
 358
 359 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 360
 361 // recut mesh with existing planeDesc
 362 void Foam::cuttingPlane::reCut
 363 (
 364     const primitiveMesh& mesh,
 365     const UList<label>& cellIdLabels
 366 )
 367 {
 368     MeshStorage::clear();
 369     cutCells_.clear();
 370
 371     scalarField dotProducts = (mesh.points() - refPoint()) & normal();
 372
 373     // Determine cells that are (probably) cut.
 374     calcCutCells(mesh, dotProducts, cellIdLabels);
 375
 376     // Determine cutPoints and return list of edge cuts.
 377     // per edge -1 or the label of the intersection point
 378     labelList edgePoint;
 379     intersectEdges(mesh, dotProducts, edgePoint);
 380
 381     // Do topological walk around cell to find closed loop.
 382     walkCellCuts(mesh, edgePoint);
 383 }
 384
 385
 386 // remap action on triangulation
 387 void Foam::cuttingPlane::remapFaces
 388 (
 389     const UList<label>& faceMap
 390 )
 391 {
 392     // recalculate the cells cut
 393     if (&faceMap && faceMap.size())
 394     {
 395         MeshStorage::remapFaces(faceMap);
 396
 397         List<label> newCutCells(faceMap.size());
 398         forAll(faceMap, faceI)
 399         {
 400             newCutCells[faceI] = cutCells_[faceMap[faceI]];
 401         }
 402         cutCells_.transfer(newCutCells);
 403     }
 404 }
 405
 406 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 407
 408 void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
 409 {
 410     // Check for assignment to self
 411     if (this == &rhs)
 412     {
 413         FatalErrorIn ("Foam::cuttingPlane::operator=(const cuttingPlane&)")
 414             << "Attempted assignment to self"
 415             << abort(FatalError);
 416     }
 417
 418     static_cast<MeshStorage&>(*this) = rhs;
 419     static_cast<plane&>(*this) = rhs;
 420     cutCells_ = rhs.cutCells();
 421 }
 422
 423
 424 // ************************************************************************* //