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-2008 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 labelList& 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 Foam::labelList Foam::cuttingPlane::intersectEdges
 108 (
 109     const primitiveMesh& mesh,
 110     const scalarField& dotProducts
 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     labelList edgePoint(edges.size(), -1);
 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             const point& p0 = points[e[0]];
 134             const point& p1 = points[e[1]];
 135
 136             scalar alpha = lineIntersect(linePointRef(p0, p1));
 137
 138             if (alpha < zeroish)
 139             {
 140                 dynCuttingPoints.append(p0);
 141             }
 142             else if (alpha > 1.0)
 143             {
 144                 dynCuttingPoints.append(p1);
 145             }
 146             else
 147             {
 148                 dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
 149             }
 150
 151             edgePoint[edgeI] = dynCuttingPoints.size() - 1;
 152         }
 153     }
 154
 155     dynCuttingPoints.shrink();
 156     cutPoints_.transfer(dynCuttingPoints);
 157     dynCuttingPoints.clear();
 158
 159     return edgePoint;
 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 labelList& 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     do
 180     {
 181         faceVerts.append(edgePoint[edgeI]);
 182
 183         // Cross edge to other face
 184         faceI = meshTools::otherFace(mesh, cellI, faceI, edgeI);
 185
 186         // Find next cut edge on face.
 187         const labelList& fEdges = mesh.faceEdges()[faceI];
 188
 189         label nextEdgeI = -1;
 190
 191         //Note: here is where we should check for whether there are more
 192         // than 2 intersections with the face (warped/non-convex face).
 193         // If so should e.g. decompose the cells on both faces and redo
 194         // the calculation.
 195
 196         forAll(fEdges, i)
 197         {
 198             label edge2I = fEdges[i];
 199
 200             if (edge2I != edgeI && edgePoint[edge2I] != -1)
 201             {
 202                 nextEdgeI = edge2I;
 203                 break;
 204             }
 205         }
 206
 207         if (nextEdgeI == -1)
 208         {
 209             // Did not find another cut edge on faceI. Do what?
 210             WarningIn("Foam::cuttingPlane::walkCell")
 211                 << "Did not find closed walk along surface of cell " << cellI
 212                 << " starting from edge " << startEdgeI
 213                 << " in " << nIter << " iterations." << nl
 214                 << "Collected cutPoints so far:" << faceVerts
 215                 << endl;
 216
 217             return false;
 218         }
 219
 220         edgeI = nextEdgeI;
 221
 222         nIter++;
 223
 224         if (nIter > 1000)
 225         {
 226             WarningIn("Foam::cuttingPlane::walkCell")
 227                 << "Did not find closed walk along surface of cell " << cellI
 228                 << " starting from edge " << startEdgeI
 229                 << " in " << nIter << " iterations." << nl
 230                 << "Collected cutPoints so far:" << faceVerts
 231                 << endl;
 232             return false;
 233         }
 234
 235     } while (edgeI != startEdgeI);
 236
 237
 238     if (faceVerts.size() >= 3)
 239     {
 240         return true;
 241     }
 242     else
 243     {
 244         WarningIn("Foam::cuttingPlane::walkCell")
 245             << "Did not find closed walk along surface of cell " << cellI
 246             << " starting from edge " << startEdgeI << nl
 247             << "Collected cutPoints so far:" << faceVerts
 248             << endl;
 249
 250         return false;
 251     }
 252 }
 253
 254
 255 // For every cut cell determine a walk through all? its cuts.
 256 void Foam::cuttingPlane::walkCellCuts
 257 (
 258     const primitiveMesh& mesh,
 259     const labelList& edgePoint
 260 )
 261 {
 262     cutFaces_.setSize(cutCells_.size());
 263     label cutFaceI = 0;
 264
 265     forAll(cutCells_, i)
 266     {
 267         label cellI = cutCells_[i];
 268
 269         // Find the starting edge to walk from.
 270         const labelList& cEdges = mesh.cellEdges()[cellI];
 271
 272         label startEdgeI = -1;
 273
 274         forAll(cEdges, cEdgeI)
 275         {
 276             label edgeI = cEdges[cEdgeI];
 277
 278             if (edgePoint[edgeI] != -1)
 279             {
 280                 startEdgeI = edgeI;
 281                 break;
 282             }
 283         }
 284
 285         // Check for the unexpected ...
 286         if (startEdgeI == -1)
 287         {
 288             FatalErrorIn("Foam::cuttingPlane::walkCellCuts")
 289                 << "Cannot find cut edge for cut cell " << cellI
 290                 << abort(FatalError);
 291         }
 292
 293         // Walk from starting edge around the circumference of the cell.
 294         DynamicList<label> faceVerts(2*mesh.faces()[cellI].size());
 295         bool okCut = walkCell
 296         (
 297             mesh,
 298             edgePoint,
 299             cellI,
 300             startEdgeI,
 301             faceVerts
 302         );
 303
 304         if (okCut)
 305         {
 306             faceVerts.shrink();
 307
 308             face cutFace(faceVerts);
 309
 310             // Orient face.
 311             if ((cutFace.normal(cutPoints_) && normal()) < 0)
 312             {
 313                 cutFace = cutFace.reverseFace();
 314             }
 315
 316             cutFaces_[cutFaceI++] = cutFace;
 317         }
 318     }
 319
 320     cutFaces_.setSize(cutFaceI);
 321 }
 322
 323
 324 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 325
 326 // Construct without cutting
 327 Foam::cuttingPlane::cuttingPlane(const plane& newPlane)
 328 :
 329     plane(newPlane)
 330 {}
 331
 332
 333 // Construct from components
 334 Foam::cuttingPlane::cuttingPlane
 335 (
 336     const primitiveMesh& mesh,
 337     const plane& newPlane
 338 )
 339 :
 340     plane(newPlane)
 341 {
 342     reCut(mesh);
 343 }
 344
 345 // Construct from mesh reference and plane, restricted to a list of cells
 346 Foam::cuttingPlane::cuttingPlane
 347 (
 348     const primitiveMesh& mesh,
 349     const plane& newPlane,
 350     const labelList& cellIdLabels
 351 )
 352 :
 353     plane(newPlane)
 354 {
 355     reCut(mesh, cellIdLabels);
 356 }
 357
 358
 359
 360 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 361
 362 // reCut mesh with existing planeDesc
 363 void Foam::cuttingPlane::reCut
 364 (
 365     const primitiveMesh& mesh
 366 )
 367 {
 368     cutCells_.clear();
 369     cutPoints_.clear();
 370     cutFaces_.clear();
 371
 372     scalarField dotProducts = (mesh.points() - refPoint()) & normal();
 373
 374     //// Perturb points cuts so edges are cut properly.
 375     //const tmp<scalarField> tdotProducts = stabilise(rawDotProducts, SMALL);
 376     //const scalarField& dotProducts = tdotProducts();
 377
 378     // Determine cells that are (probably) cut.
 379     calcCutCells(mesh, dotProducts);
 380
 381     // Determine cutPoints and return list of edge cuts. (per edge -1 or
 382     // the label of the intersection point (in cutPoints_)
 383     labelList edgePoint(intersectEdges(mesh, dotProducts));
 384
 385     // Do topological walk around cell to find closed loop.
 386     walkCellCuts(mesh, edgePoint);
 387 }
 388
 389
 390 // recut mesh with existing planeDesc
 391 void Foam::cuttingPlane::reCut
 392 (
 393     const primitiveMesh& mesh,
 394     const labelList& cellIdLabels
 395 )
 396 {
 397     cutCells_.clear();
 398     cutPoints_.clear();
 399     cutFaces_.clear();
 400
 401     scalarField dotProducts = (mesh.points() - refPoint()) & normal();
 402
 403     // Determine cells that are (probably) cut.
 404     calcCutCells(mesh, dotProducts, cellIdLabels);
 405
 406     // Determine cutPoints and return list of edge cuts. (per edge -1 or
 407     // the label of the intersection point (in cutPoints_)
 408     labelList edgePoint(intersectEdges(mesh, dotProducts));
 409
 410     // Do topological walk around cell to find closed loop.
 411     walkCellCuts(mesh, edgePoint);
 412 }
 413
 414
 415
 416 // Return plane used
 417 const Foam::plane& Foam::cuttingPlane::planeDesc() const
 418 {
 419     return static_cast<const plane&>(*this);
 420 }
 421
 422
 423 // Return vectorField of cutting points
 424 const Foam::pointField& Foam::cuttingPlane::points() const
 425 {
 426     return cutPoints_;
 427 }
 428
 429
 430 // Return unallocFaceList of points in cells
 431 const Foam::faceList& Foam::cuttingPlane::faces() const
 432 {
 433     return cutFaces_;
 434 }
 435
 436
 437 // Return labelList of cut cells
 438 const Foam::labelList& Foam::cuttingPlane::cells() const
 439 {
 440     return cutCells_;
 441 }
 442
 443
 444 bool Foam::cuttingPlane::cut()
 445 {
 446     if (cutCells_.size() > 0)
 447     {
 448         return true;
 449     }
 450     else
 451     {
 452         return false;
 453     }
 454 }
 455
 456
 457 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 458
 459 void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
 460 {
 461     // Check for assignment to self
 462     if (this == &rhs)
 463     {
 464         FatalErrorIn ("Foam::cuttingPlane::operator=(const cuttingPlane&)")
 465             << "Attempted assignment to self"
 466             << abort(FatalError);
 467     }
 468
 469     static_cast<plane&>(*this) = rhs;
 470
 471     cutCells_ = rhs.cells();
 472     cutPoints_ = rhs.points();
 473     cutFaces_ = rhs.faces();
 474 }
 475
 476
 477 // ************************************************************************* //