src/sampling/cuttingPlane/cuttingPlane.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "objectRegistry.H"
  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 localSmall = Foam::SMALL;
  38 const Foam::scalar localSmallPositive = 1e-3*Foam::SMALL;
  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.empty())
  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.empty())
  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                 (
  84                     dotProducts[e[0]] < localSmall
  85                  && dotProducts[e[1]] > localSmallPositive
  86                 )
  87              || (
  88                     dotProducts[e[1]] < localSmall
  89                  && dotProducts[e[0]] > localSmallPositive
  90                 )
  91             )
  92             {
  93                 nCutEdges++;
  94
  95                 if (nCutEdges > 2)
  96                 {
  97                     cutCells_[cutcellI++] = cellI;
  98
  99                     break;
 100                 }
 101             }
 102         }
 103     }
 104
 105     // Set correct list size
 106     cutCells_.setSize(cutcellI);
 107 }
 108
 109
 110 // Determine for each edge the intersection point. Calculates
 111 // - cutPoints_ : coordinates of all intersection points
 112 // - edgePoint  : per edge -1 or the index into cutPoints
 113 void Foam::cuttingPlane::intersectEdges
 114 (
 115     const primitiveMesh& mesh,
 116     const scalarField& dotProducts,
 117     List<label>& edgePoint
 118 )
 119 {
 120     // Use the dotProducts to find out the cut edges.
 121     const edgeList& edges = mesh.edges();
 122     const pointField& points = mesh.points();
 123
 124     // Per edge -1 or the label of the intersection point
 125     edgePoint.setSize(edges.size());
 126
 127     DynamicList<point> dynCuttingPoints(4*cutCells_.size());
 128
 129     forAll(edges, edgeI)
 130     {
 131         const edge& e = edges[edgeI];
 132
 133         if
 134         (
 135             (
 136                 dotProducts[e[0]] < localSmall
 137              && dotProducts[e[1]] > localSmallPositive
 138             )
 139          || (
 140                 dotProducts[e[1]] < localSmall
 141              && dotProducts[e[0]] > localSmallPositive
 142             )
 143         )
 144         {
 145             // Edge is cut
 146             edgePoint[edgeI] = dynCuttingPoints.size();
 147
 148             const point& p0 = points[e[0]];
 149             const point& p1 = points[e[1]];
 150
 151             scalar alpha = lineIntersect(linePointRef(p0, p1));
 152
 153             if (alpha < localSmall)
 154             {
 155                 dynCuttingPoints.append(p0);
 156             }
 157             else if (alpha >= 1.0)
 158             {
 159                 dynCuttingPoints.append(p1);
 160             }
 161             else
 162             {
 163                 dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
 164             }
 165         }
 166         else
 167         {
 168             edgePoint[edgeI] = -1;
 169         }
 170     }
 171
 172     this->storedPoints().transfer(dynCuttingPoints);
 173 }
 174
 175
 176 // Coming from startEdgeI cross the edge to the other face
 177 // across to the next cut edge.
 178 bool Foam::cuttingPlane::walkCell
 179 (
 180     const primitiveMesh& mesh,
 181     const UList<label>& edgePoint,
 182     const label cellI,
 183     const label startEdgeI,
 184     DynamicList<label>& faceVerts
 185 )
 186 {
 187     label faceI = -1;
 188     label edgeI = startEdgeI;
 189
 190     label nIter = 0;
 191
 192     faceVerts.clear();
 193     do
 194     {
 195         faceVerts.append(edgePoint[edgeI]);
 196
 197         // Cross edge to other face
 198         faceI = meshTools::otherFace(mesh, cellI, faceI, edgeI);
 199
 200         // Find next cut edge on face.
 201         const labelList& fEdges = mesh.faceEdges()[faceI];
 202
 203         label nextEdgeI = -1;
 204
 205         //Note: here is where we should check for whether there are more
 206         // than 2 intersections with the face (warped/non-convex face).
 207         // If so should e.g. decompose the cells on both faces and redo
 208         // the calculation.
 209
 210         forAll(fEdges, i)
 211         {
 212             label edge2I = fEdges[i];
 213
 214             if (edge2I != edgeI && edgePoint[edge2I] != -1)
 215             {
 216                 nextEdgeI = edge2I;
 217                 break;
 218             }
 219         }
 220
 221         if (nextEdgeI == -1)
 222         {
 223             // Did not find another cut edge on faceI. Do what?
 224             WarningIn("Foam::cuttingPlane::walkCell")
 225                 << "Did not find closed walk along surface of cell " << cellI
 226                 << " starting from edge " << startEdgeI
 227                 << " in " << nIter << " iterations." << nl
 228                 << "Collected cutPoints so far:" << faceVerts
 229                 << endl;
 230
 231             return false;
 232         }
 233
 234         edgeI = nextEdgeI;
 235
 236         nIter++;
 237
 238         if (nIter > 1000)
 239         {
 240             WarningIn("Foam::cuttingPlane::walkCell")
 241                 << "Did not find closed walk along surface of cell " << cellI
 242                 << " starting from edge " << startEdgeI
 243                 << " in " << nIter << " iterations." << nl
 244                 << "Collected cutPoints so far:" << faceVerts
 245                 << endl;
 246             return false;
 247         }
 248
 249     } while (edgeI != startEdgeI);
 250
 251
 252     if (faceVerts.size() >= 3)
 253     {
 254         return true;
 255     }
 256     else
 257     {
 258         WarningIn("Foam::cuttingPlane::walkCell")
 259             << "Did not find closed walk along surface of cell " << cellI
 260             << " starting from edge " << startEdgeI << nl
 261             << "Collected cutPoints so far:" << faceVerts
 262             << endl;
 263
 264         return false;
 265     }
 266 }
 267
 268
 269 // For every cut cell determine a walk through all? its cuts.
 270 void Foam::cuttingPlane::walkCellCuts
 271 (
 272     const primitiveMesh& mesh,
 273     const UList<label>& edgePoint
 274 )
 275 {
 276     const pointField& cutPoints = this->points();
 277
 278     // use dynamic lists to handle triangulation and/or missed cuts
 279     DynamicList<face>  dynCutFaces(cutCells_.size());
 280     DynamicList<label> dynCutCells(cutCells_.size());
 281
 282     // scratch space for calculating the face vertices
 283     DynamicList<label> faceVerts(10);
 284
 285     forAll(cutCells_, i)
 286     {
 287         label cellI = cutCells_[i];
 288
 289         // Find the starting edge to walk from.
 290         const labelList& cEdges = mesh.cellEdges()[cellI];
 291
 292         label startEdgeI = -1;
 293
 294         forAll(cEdges, cEdgeI)
 295         {
 296             label edgeI = cEdges[cEdgeI];
 297
 298             if (edgePoint[edgeI] != -1)
 299             {
 300                 startEdgeI = edgeI;
 301                 break;
 302             }
 303         }
 304
 305         // Check for the unexpected ...
 306         if (startEdgeI == -1)
 307         {
 308             FatalErrorIn("Foam::cuttingPlane::walkCellCuts")
 309                 << "Cannot find cut edge for cut cell " << cellI
 310                 << abort(FatalError);
 311         }
 312
 313         // Walk from starting edge around the circumference of the cell.
 314         bool okCut = walkCell
 315         (
 316             mesh,
 317             edgePoint,
 318             cellI,
 319             startEdgeI,
 320             faceVerts
 321         );
 322
 323         if (okCut)
 324         {
 325             face f(faceVerts);
 326
 327             // Orient face to point in the same direction as the plane normal
 328             if ((f.normal(cutPoints) && normal()) < 0)
 329             {
 330                 f = f.reverseFace();
 331             }
 332
 333             // the cut faces are usually quite ugly, so always triangulate
 334             label nTri = f.triangles(cutPoints, dynCutFaces);
 335             while (nTri--)
 336             {
 337                 dynCutCells.append(cellI);
 338             }
 339         }
 340     }
 341
 342     this->storedFaces().transfer(dynCutFaces);
 343     cutCells_.transfer(dynCutCells);
 344 }
 345
 346
 347 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 348
 349 // Construct without cutting
 350 Foam::cuttingPlane::cuttingPlane(const plane& pln)
 351 :
 352     plane(pln)
 353 {}
 354
 355
 356 // Construct from plane and mesh reference, restricted to a list of cells
 357 Foam::cuttingPlane::cuttingPlane
 358 (
 359     const plane& pln,
 360     const primitiveMesh& mesh,
 361     const UList<label>& cellIdLabels
 362 )
 363 :
 364     plane(pln)
 365 {
 366     reCut(mesh, cellIdLabels);
 367 }
 368
 369
 370
 371 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 372
 373 // recut mesh with existing planeDesc
 374 void Foam::cuttingPlane::reCut
 375 (
 376     const primitiveMesh& mesh,
 377     const UList<label>& cellIdLabels
 378 )
 379 {
 380     MeshStorage::clear();
 381     cutCells_.clear();
 382
 383     scalarField dotProducts = (mesh.points() - refPoint()) & normal();
 384
 385     // Determine cells that are (probably) cut.
 386     calcCutCells(mesh, dotProducts, cellIdLabels);
 387
 388     // Determine cutPoints and return list of edge cuts.
 389     // per edge -1 or the label of the intersection point
 390     labelList edgePoint;
 391     intersectEdges(mesh, dotProducts, edgePoint);
 392
 393     // Do topological walk around cell to find closed loop.
 394     walkCellCuts(mesh, edgePoint);
 395 }
 396
 397
 398 // remap action on triangulation
 399 void Foam::cuttingPlane::remapFaces
 400 (
 401     const UList<label>& faceMap
 402 )
 403 {
 404     // recalculate the cells cut
 405     if (!faceMap.empty())
 406     {
 407         MeshStorage::remapFaces(faceMap);
 408
 409         List<label> newCutCells(faceMap.size());
 410         forAll(faceMap, faceI)
 411         {
 412             newCutCells[faceI] = cutCells_[faceMap[faceI]];
 413         }
 414         cutCells_.transfer(newCutCells);
 415     }
 416 }
 417
 418 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 419
 420 void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
 421 {
 422     // Check for assignment to self
 423     if (this == &rhs)
 424     {
 425         FatalErrorIn ("Foam::cuttingPlane::operator=(const cuttingPlane&)")
 426             << "Attempted assignment to self"
 427             << abort(FatalError);
 428     }
 429
 430     static_cast<MeshStorage&>(*this) = rhs;
 431     static_cast<plane&>(*this) = rhs;
 432     cutCells_ = rhs.cutCells();
 433 }
 434
 435
 436 // ************************************************************************* //