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