src/dynamicMesh/meshCut/cellLooper/geomCellLooper.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 "geomCellLooper.H"
  28 #include "polyMesh.H"
  29 #include "DynamicList.H"
  30 #include "plane.H"
  31 #include "meshTools.H"
  32 #include "SortableList.H"
  33 #include "triSurfaceTools.H"
  34 #include "HashSet.H"
  35 #include "ListOps.H"
  36 #include "transform.H"
  37
  38 #include "addToRunTimeSelectionTable.H"
  39
  40
  41 // Extension factor of edges to make sure we catch intersections through
  42 // edge endpoints
  43 const Foam::scalar Foam::geomCellLooper::pointEqualTol_ = 1E-3;
  44
  45
  46 // Snap cuts through edges onto edge endpoints. Fraction of edge length.
  47 Foam::scalar Foam::geomCellLooper::snapTol_ = 0.1;
  48
  49
  50 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  51
  52 namespace Foam
  53 {
  54
  55 defineTypeNameAndDebug(geomCellLooper, 0);
  56 addToRunTimeSelectionTable(cellLooper, geomCellLooper, word);
  57
  58 }
  59
  60
  61 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  62
  63 Foam::scalar Foam::geomCellLooper::minEdgeLen(const label vertI) const
  64 {
  65     scalar minLen = GREAT;
  66
  67     const labelList& pEdges = mesh().pointEdges()[vertI];
  68
  69     forAll(pEdges, pEdgeI)
  70     {
  71         const edge& e = mesh().edges()[pEdges[pEdgeI]];
  72
  73         minLen = min(minLen, e.mag(mesh().points()));
  74     }
  75     return minLen;
  76 }
  77
  78
  79 // Cut edge with plane. Return true and set weight to fraction between
  80 // edge-start and edge-end
  81 bool Foam::geomCellLooper::cutEdge
  82 (
  83     const plane& cutPlane,
  84     const label edgeI,
  85     scalar& weight
  86 ) const
  87 {
  88     const pointField& pts = mesh().points();
  89
  90     const edge& e = mesh().edges()[edgeI];
  91
  92     scalar s = cutPlane.normalIntersect(pts[e.start()], e.vec(pts));
  93
  94     if ((s > -pointEqualTol_) && (s < 1 + pointEqualTol_))
  95     {
  96         weight = s;
  97
  98         return true;
  99     }
 100     else
 101     {
 102         // Make sure we don't use this value
 103         weight = -GREAT;
 104
 105         return false;
 106     }
 107 }
 108
 109
 110 // If edge close enough to endpoint snap to endpoint.
 111 Foam::label Foam::geomCellLooper::snapToVert
 112 (
 113     const scalar tol,
 114     const label edgeI,
 115     const scalar weight
 116 ) const
 117 {
 118     const edge& e = mesh().edges()[edgeI];
 119
 120     if (weight < tol)
 121     {
 122         return e.start();
 123     }
 124     else if (weight > (1-tol))
 125     {
 126         return e.end();
 127     }
 128     else
 129     {
 130         return -1;
 131     }
 132 }
 133
 134
 135 void Foam::geomCellLooper::getBase(const vector& n, vector& e0, vector& e1)
 136  const
 137 {
 138     // Guess for vector normal to n.
 139     vector base(1,0,0);
 140
 141     scalar nComp = n & base;
 142
 143     if (mag(nComp) > 0.8)
 144     {
 145         // Was bad guess. Try with different vector.
 146
 147         base.x() = 0;
 148         base.y() = 1;
 149
 150         nComp = n & base;
 151
 152         if (mag(nComp) > 0.8)
 153         {
 154             base.y() = 0;
 155             base.z() = 1;
 156
 157             nComp = n & base;
 158         }
 159     }
 160
 161
 162     // Use component normal to n as base vector.
 163     e0 = base - nComp*n;
 164
 165     e0 /= mag(e0) + VSMALL;
 166
 167     e1 = n ^ e0;
 168
 169     //Pout<< "Coord system:" << endl
 170     //    << "    n  : " << n << ' ' << mag(n) << endl
 171     //    << "    e0 : " << e0 << ' ' << mag(e0) << endl
 172     //    << "    e1 : " << e1 << ' ' << mag(e1) << endl
 173     //    << endl;
 174 }
 175
 176
 177 // Return true if the cut edge at loop[index] is preceded by cuts through
 178 // the edge end points.
 179 bool Foam::geomCellLooper::edgeEndsCut
 180 (
 181     const labelList& loop,
 182     const label index
 183 ) const
 184 {
 185     label edgeI = getEdge(loop[index]);
 186
 187     const edge& e = mesh().edges()[edgeI];
 188
 189     const label prevCut = loop[loop.rcIndex(index)];
 190     const label nextCut = loop[loop.fcIndex(index)];
 191
 192     if (!isEdge(prevCut) && !isEdge(nextCut))
 193     {
 194         // Cuts before and after are both vertices. Check if both
 195         // the edge endpoints
 196         label v0 = getVertex(prevCut);
 197         label v1 = getVertex(nextCut);
 198
 199         if
 200         (
 201             (v0 == e[0] && v1 == e[1])
 202          || (v0 == e[1] && v1 == e[0])
 203         )
 204         {
 205             return true;
 206         }
 207     }
 208     return false;
 209 }
 210
 211
 212 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 213
 214 // Construct from components
 215 Foam::geomCellLooper::geomCellLooper(const polyMesh& mesh)
 216 :
 217     cellLooper(mesh)
 218 {}
 219
 220
 221 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 222
 223 Foam::geomCellLooper::~geomCellLooper()
 224 {}
 225
 226
 227 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 228
 229 bool Foam::geomCellLooper::cut
 230 (
 231     const vector& refDir,
 232     const label cellI,
 233     const boolList& vertIsCut,
 234     const boolList& edgeIsCut,
 235     const scalarField& edgeWeight,
 236
 237     labelList& loop,
 238     scalarField& loopWeights
 239 ) const
 240 {
 241     // Cut through cell centre normal to refDir.
 242     return cut
 243     (
 244         plane(mesh().cellCentres()[cellI], refDir),
 245         cellI,
 246         vertIsCut,
 247         edgeIsCut,
 248         edgeWeight,
 249         loop,
 250         loopWeights
 251     );
 252 }
 253
 254
 255 bool Foam::geomCellLooper::cut
 256 (
 257     const plane& cutPlane,
 258     const label cellI,
 259     const boolList&,
 260     const boolList&,
 261     const scalarField&,
 262
 263     labelList& loop,
 264     scalarField& loopWeights
 265 ) const
 266 {
 267     const pointField& points = mesh().points();
 268     const edgeList& edges = mesh().edges();
 269
 270     // Find all cuts through edges.
 271     // Special cases:
 272     // - edge cut close to endpoint. Snap to endpoint.
 273     // - edge endpoint close to plane (but edge not cut). Snap to endpoint.
 274     // - both endpoints close to plane. Edge parallel to plane. No need to
 275     //   cut to edge.
 276     // Note: any snap-to-point check must be based on all edges using a point
 277     // since otherwise cut through close to point but on neighbouring edge
 278     // might not be snapped.
 279
 280     // Size overly big.
 281     label nEstCuts = 2*mesh().cells()[cellI].size();
 282
 283     DynamicList<label> localLoop(nEstCuts);
 284     DynamicList<scalar> localLoopWeights(nEstCuts);
 285
 286     // Points checked. Used to make sure we don't cut edge and edge endpoints
 287     // at the same time.
 288     labelHashSet checkedPoints(nEstCuts);
 289
 290     const labelList& cellEdges = mesh().cellEdges()[cellI];
 291
 292     forAll(cellEdges, i)
 293     {
 294         label edgeI = cellEdges[i];
 295
 296         const edge& e = edges[edgeI];
 297
 298         bool useStart = false;
 299
 300         bool useEnd = false;
 301
 302         //
 303         // Check distance of endpoints to cutPlane
 304         //
 305
 306         if (!checkedPoints.found(e.start()))
 307         {
 308             checkedPoints.insert(e.start());
 309
 310             scalar typStartLen = pointEqualTol_ * minEdgeLen(e.start());
 311
 312             // Check distance of startPt to plane.
 313             if (cutPlane.distance(points[e.start()]) < typStartLen)
 314             {
 315                 // Use point.
 316                 localLoop.append(vertToEVert(e.start()));
 317                 localLoopWeights.append(-GREAT);
 318
 319                 useStart = true;
 320             }
 321         }
 322         if (!checkedPoints.found(e.end()))
 323         {
 324             checkedPoints.insert(e.end());
 325
 326             scalar typEndLen = pointEqualTol_ * minEdgeLen(e.end());
 327
 328             // Check distance of endPt to plane.
 329             if (cutPlane.distance(points[e.end()]) < typEndLen)
 330             {
 331                 // Use point.
 332                 localLoop.append(vertToEVert(e.end()));
 333                 localLoopWeights.append(-GREAT);
 334
 335                 useEnd = true;
 336             }
 337         }
 338
 339         //
 340         // Check cut of edge
 341         //
 342
 343         if (!useEnd && !useStart)
 344         {
 345             // Edge end points not close to plane so edge not close to
 346             // plane. Cut edge.
 347             scalar cutWeight;
 348
 349             if (cutEdge(cutPlane, edgeI, cutWeight))
 350             {
 351                 // Snap edge cut to vertex.
 352                 label cutVertI = snapToVert(snapTol_, edgeI, cutWeight);
 353
 354                 if (cutVertI == -1)
 355                 {
 356                     // Proper cut of edge
 357                     localLoop.append(edgeToEVert(edgeI));
 358                     localLoopWeights.append(cutWeight);
 359                 }
 360                 else
 361                 {
 362                     // Cut through edge end point. Might be duplicate
 363                     // since connected edge might also be snapped to same
 364                     // endpoint so only insert if unique.
 365                     label cut = vertToEVert(cutVertI);
 366
 367                     if (findIndex(localLoop, cut) == -1)
 368                     {
 369                         localLoop.append(vertToEVert(cutVertI));
 370                         localLoopWeights.append(-GREAT);
 371                     }
 372                 }
 373             }
 374         }
 375     }
 376
 377     if (localLoop.size() <= 2)
 378     {
 379         return false;
 380     }
 381
 382     localLoop.shrink();
 383     localLoopWeights.shrink();
 384
 385
 386     // Get points on loop and centre of loop
 387     pointField loopPoints(localLoop.size());
 388     point ctr(vector::zero);
 389
 390     forAll(localLoop, i)
 391     {
 392         loopPoints[i] = coord(localLoop[i], localLoopWeights[i]);
 393         ctr += loopPoints[i];
 394     }
 395     ctr /= localLoop.size();
 396
 397
 398     // Get base vectors of coordinate system normal to refDir
 399     vector e0, e1;
 400     getBase(cutPlane.normal(), e0, e1);
 401
 402
 403     // Get sorted angles from point on loop to centre of loop.
 404     SortableList<scalar> sortedAngles(localLoop.size());
 405
 406     forAll(sortedAngles, i)
 407     {
 408         vector toCtr(loopPoints[i] - ctr);
 409         toCtr /= mag(toCtr);
 410
 411         sortedAngles[i] = pseudoAngle(e0, e1, toCtr);
 412     }
 413     sortedAngles.sort();
 414
 415     loop.setSize(localLoop.size());
 416     loopWeights.setSize(loop.size());
 417
 418
 419     // Fill loop and loopweights according to sorted angles
 420
 421     const labelList& indices = sortedAngles.indices();
 422
 423     forAll(indices, i)
 424     {
 425         loop[i] = localLoop[indices[i]];
 426         loopWeights[i] = localLoopWeights[indices[i]];
 427     }
 428
 429
 430     // Check for cut edges along already cut vertices.
 431     bool filterLoop = false;
 432
 433     forAll(loop, i)
 434     {
 435         label cut = loop[i];
 436
 437         if (isEdge(cut) && edgeEndsCut(loop, i))
 438         {
 439             filterLoop = true;
 440         }
 441     }
 442
 443     if (filterLoop)
 444     {
 445         // Found edges in loop where both end points are also in loop. This
 446         // is superfluous information so we can remove it.
 447
 448         labelList filteredLoop(loop.size());
 449         scalarField filteredLoopWeights(loopWeights.size());
 450         label filterI = 0;
 451
 452         forAll(loop, i)
 453         {
 454             label cut = loop[i];
 455
 456             if (isEdge(cut) && edgeEndsCut(loop, i))
 457             {
 458                 // Superfluous cut. Edge end points cut and edge itself as well.
 459             }
 460             else
 461             {
 462                 filteredLoop[filterI] = loop[i];
 463                 filteredLoopWeights[filterI] = loopWeights[i];
 464                 filterI++;
 465             }
 466         }
 467         filteredLoop.setSize(filterI);
 468         filteredLoopWeights.setSize(filterI);
 469
 470         loop.transfer(filteredLoop);
 471         loopWeights.transfer(filteredLoopWeights);
 472     }
 473
 474     if (debug&2)
 475     {
 476         Pout<< "cell:" << cellI << endl;
 477         forAll(loop, i)
 478         {
 479             Pout<< "At angle:" << sortedAngles[i] << endl
 480                 << "    cut:";
 481
 482             writeCut(Pout, loop[i], loopWeights[i]);
 483
 484             Pout<< "  coord:" << coord(loop[i], loopWeights[i]);
 485
 486             Pout<< endl;
 487         }
 488     }
 489
 490     return true;
 491 }
 492
 493
 494 // ************************************************************************* //