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-2010 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
  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 localScope
  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 UList<label>& 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 UList<label>& 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 UList<label>& edgePoint
 261 )
 262 {
 263     const pointField& cutPoints = this->points();
 264
 265     // use dynamic lists to handle triangulation and/or missed cuts
 266     DynamicList<face>  dynCutFaces(cutCells_.size());
 267     DynamicList<label> dynCutCells(cutCells_.size());
 268
 269     // scratch space for calculating the face vertices
 270     DynamicList<label> faceVerts(10);
 271
 272     forAll(cutCells_, i)
 273     {
 274         label cellI = cutCells_[i];
 275
 276         // Find the starting edge to walk from.
 277         const labelList& cEdges = mesh.cellEdges()[cellI];
 278
 279         label startEdgeI = -1;
 280
 281         forAll(cEdges, cEdgeI)
 282         {
 283             label edgeI = cEdges[cEdgeI];
 284
 285             if (edgePoint[edgeI] != -1)
 286             {
 287                 startEdgeI = edgeI;
 288                 break;
 289             }
 290         }
 291
 292         // Check for the unexpected ...
 293         if (startEdgeI == -1)
 294         {
 295             FatalErrorIn("Foam::cuttingPlane::walkCellCuts")
 296                 << "Cannot find cut edge for cut cell " << cellI
 297                 << abort(FatalError);
 298         }
 299
 300         // Walk from starting edge around the circumference of the cell.
 301         bool okCut = walkCell
 302         (
 303             mesh,
 304             edgePoint,
 305             cellI,
 306             startEdgeI,
 307             faceVerts
 308         );
 309
 310         if (okCut)
 311         {
 312             face f(faceVerts);
 313
 314             // Orient face to point in the same direction as the plane normal
 315             if ((f.normal(cutPoints) & normal()) < 0)
 316             {
 317                 f = f.reverseFace();
 318             }
 319
 320             // the cut faces are usually quite ugly, so always triangulate
 321             label nTri = f.triangles(cutPoints, dynCutFaces);
 322             while (nTri--)
 323             {
 324                 dynCutCells.append(cellI);
 325             }
 326         }
 327     }
 328
 329     this->storedFaces().transfer(dynCutFaces);
 330     cutCells_.transfer(dynCutCells);
 331 }
 332
 333
 334 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 335
 336 // Construct without cutting
 337 Foam::cuttingPlane::cuttingPlane(const plane& pln)
 338 :
 339     plane(pln)
 340 {}
 341
 342
 343 // Construct from plane and mesh reference, restricted to a list of cells
 344 Foam::cuttingPlane::cuttingPlane
 345 (
 346     const plane& pln,
 347     const primitiveMesh& mesh,
 348     const UList<label>& cellIdLabels
 349 )
 350 :
 351     plane(pln)
 352 {
 353     reCut(mesh, cellIdLabels);
 354 }
 355
 356
 357
 358 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 359
 360 // recut mesh with existing planeDesc
 361 void Foam::cuttingPlane::reCut
 362 (
 363     const primitiveMesh& mesh,
 364     const UList<label>& cellIdLabels
 365 )
 366 {
 367     MeshStorage::clear();
 368     cutCells_.clear();
 369
 370     scalarField dotProducts = (mesh.points() - refPoint()) & normal();
 371
 372     // Determine cells that are (probably) cut.
 373     calcCutCells(mesh, dotProducts, cellIdLabels);
 374
 375     // Determine cutPoints and return list of edge cuts.
 376     // per edge -1 or the label of the intersection point
 377     labelList edgePoint;
 378     intersectEdges(mesh, dotProducts, edgePoint);
 379
 380     // Do topological walk around cell to find closed loop.
 381     walkCellCuts(mesh, edgePoint);
 382 }
 383
 384
 385 // remap action on triangulation
 386 void Foam::cuttingPlane::remapFaces
 387 (
 388     const UList<label>& faceMap
 389 )
 390 {
 391     // recalculate the cells cut
 392     if (&faceMap && faceMap.size())
 393     {
 394         MeshStorage::remapFaces(faceMap);
 395
 396         List<label> newCutCells(faceMap.size());
 397         forAll(faceMap, faceI)
 398         {
 399             newCutCells[faceI] = cutCells_[faceMap[faceI]];
 400         }
 401         cutCells_.transfer(newCutCells);
 402     }
 403 }
 404
 405 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 406
 407 void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
 408 {
 409     // Check for assignment to self
 410     if (this == &rhs)
 411     {
 412         FatalErrorIn ("Foam::cuttingPlane::operator=(const cuttingPlane&)")
 413             << "Attempted assignment to self"
 414             << abort(FatalError);
 415     }
 416
 417     static_cast<MeshStorage&>(*this) = rhs;
 418     static_cast<plane&>(*this) = rhs;
 419     cutCells_ = rhs.cutCells();
 420 }
 421
 422
 423 // ************************************************************************* //