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