src/dynamicMesh/polyTopoChange/polyTopoChange/faceCollapser.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 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "faceCollapser.H"
  30 #include "polyMesh.H"
  31 #include "polyTopoChange.H"
  32 #include "polyModifyPoint.H"
  33 #include "polyModifyFace.H"
  34 #include "polyRemoveFace.H"
  35 #include "SortableList.H"
  36 #include "meshTools.H"
  37 #include "OFstream.H"
  38
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 // Insert labelList into labelHashSet. Optional excluded element.
  43 void Foam::faceCollapser::insert
  44 (
  45     const labelList& elems,
  46     const label excludeElem,
  47     labelHashSet& set
  48 )
  49 {
  50     forAll(elems, i)
  51     {
  52         if (elems[i] != excludeElem)
  53         {
  54             set.insert(elems[i]);
  55         }
  56     }
  57 }
  58
  59
  60 // Find edge amongst candidate edges. FatalError if none.
  61 Foam::label Foam::faceCollapser::findEdge
  62 (
  63     const edgeList& edges,
  64     const labelList& edgeLabels,
  65     const label v0,
  66     const label v1
  67 )
  68 {
  69     forAll(edgeLabels, i)
  70     {
  71         label edgeI = edgeLabels[i];
  72
  73         const edge& e = edges[edgeI];
  74
  75         if
  76         (
  77             (e[0] == v0 && e[1] == v1)
  78          || (e[0] == v1 && e[1] == v0)
  79         )
  80         {
  81             return edgeI;
  82         }
  83     }
  84
  85     FatalErrorIn("findEdge") << "Cannot find edge between vertices " << v0
  86         << " and " << v1 << " in edge labels " << edgeLabels
  87         << abort(FatalError);
  88
  89     return -1;
  90 }
  91
  92
  93 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  94
  95 // Replace vertices in face
  96 void Foam::faceCollapser::filterFace
  97 (
  98     const Map<labelList>& splitEdges,
  99     const label faceI,
 100     polyTopoChange& meshMod
 101 ) const
 102 {
 103     const face& f = mesh_.faces()[faceI];
 104     const labelList& fEdges = mesh_.faceEdges()[faceI];
 105
 106     // Space for replaced vertices and split edges.
 107     DynamicList<label> newFace(10 * f.size());
 108
 109     forAll(f, fp)
 110     {
 111         label v0 = f[fp];
 112
 113         newFace.append(v0);
 114
 115         // Look ahead one to get edge.
 116         label fp1 = f.fcIndex(fp);
 117
 118         label v1 = f[fp1];
 119
 120         // Get split on edge if any.
 121         label edgeI = findEdge(mesh_.edges(), fEdges, v0, v1);
 122
 123         Map<labelList>::const_iterator edgeFnd =
 124             splitEdges.find(edgeI);
 125
 126         if (edgeFnd != splitEdges.end())
 127         {
 128             // edgeI has been split (by introducing new vertices).
 129             // Insert new vertices in face in correct order
 130             // (splitEdges was constructed to be from edge.start() to end())
 131
 132             const labelList& extraVerts = edgeFnd();
 133
 134             if (v0 == mesh_.edges()[edgeI].start())
 135             {
 136                 forAll(extraVerts, i)
 137                 {
 138                     newFace.append(extraVerts[i]);
 139                 }
 140             }
 141             else
 142             {
 143                 forAllReverse(extraVerts, i)
 144                 {
 145                     newFace.append(extraVerts[i]);
 146                 }
 147             }
 148         }
 149     }
 150     face newF(newFace.shrink());
 151
 152     //Pout<< "Modifying face:" << faceI << " from " << f << " to " << newFace
 153     //    << endl;
 154
 155     if (newF != f)
 156     {
 157         label nei = -1;
 158
 159         label patchI = -1;
 160
 161         if (mesh_.isInternalFace(faceI))
 162         {
 163             nei = mesh_.faceNeighbour()[faceI];
 164         }
 165         else
 166         {
 167             patchI = mesh_.boundaryMesh().whichPatch(faceI);
 168         }
 169
 170         // Get current zone info
 171         label zoneID = mesh_.faceZones().whichZone(faceI);
 172
 173         bool zoneFlip = false;
 174
 175         if (zoneID >= 0)
 176         {
 177             const faceZone& fZone = mesh_.faceZones()[zoneID];
 178
 179             zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
 180         }
 181
 182         meshMod.setAction
 183         (
 184             polyModifyFace
 185             (
 186                 newF,                       // modified face
 187                 faceI,                      // label of face being modified
 188                 mesh_.faceOwner()[faceI],   // owner
 189                 nei,                        // neighbour
 190                 false,                      // face flip
 191                 patchI,                     // patch for face
 192                 false,                      // remove from zone
 193                 zoneID,                     // zone for face
 194                 zoneFlip                    // face flip in zone
 195             )
 196         );
 197     }
 198 }
 199
 200
 201 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 202
 203 // Construct from mesh
 204 Foam::faceCollapser::faceCollapser(const polyMesh& mesh)
 205 :
 206     mesh_(mesh)
 207 {}
 208
 209
 210 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 211
 212 void Foam::faceCollapser::setRefinement
 213 (
 214     const labelList& faceLabels,
 215     const labelList& fpStart,
 216     const labelList& fpEnd,
 217     polyTopoChange& meshMod
 218 ) const
 219 {
 220     const pointField& points = mesh_.points();
 221     const edgeList& edges = mesh_.edges();
 222     const faceList& faces = mesh_.faces();
 223     const labelListList& edgeFaces = mesh_.edgeFaces();
 224
 225
 226     // From split edge to newly introduced point(s). Can be more than one per
 227     // edge!
 228     Map<labelList> splitEdges(faceLabels.size());
 229
 230     // Mark faces in any way affect by modifying points/edges. Used later on
 231     // to prevent having to redo all faces.
 232     labelHashSet affectedFaces(4*faceLabels.size());
 233
 234
 235     //
 236     // Add/remove vertices and construct mapping
 237     //
 238
 239     forAll(faceLabels, i)
 240     {
 241         const label faceI = faceLabels[i];
 242
 243         const face& f = faces[faceI];
 244
 245         const label fpA = fpStart[i];
 246         const label fpB = fpEnd[i];
 247
 248         const point& pA = points[f[fpA]];
 249         const point& pB = points[f[fpB]];
 250
 251         Pout<< "Face:" << f << " collapsed to fp:" << fpA << ' '  << fpB
 252             << " with points:" << pA << ' ' << pB
 253             << endl;
 254
 255         // Create line from fpA to fpB
 256         linePointRef lineAB(pA, pB);
 257
 258         // Get projections of all vertices onto line.
 259
 260         // Distance(squared) to pA for every point on face.
 261         SortableList<scalar> dist(f.size());
 262
 263         dist[fpA] = 0;
 264         dist[fpB] = magSqr(pB - pA);
 265
 266         // Step from fpA to fpB
 267         // ~~~~~~~~~~~~~~~~~~~~
 268         // (by incrementing)
 269
 270         label fpMin1 = fpA;
 271         label fp = f.fcIndex(fpMin1);
 272
 273         while (fp != fpB)
 274         {
 275             // See where fp sorts. Make sure it is above fpMin1!
 276             pointHit near = lineAB.nearestDist(points[f[fp]]);
 277
 278             scalar w = magSqr(near.rawPoint() - pA);
 279
 280             if (w <= dist[fpMin1])
 281             {
 282                 // Offset.
 283                 w = dist[fpMin1] + 1E-6*(dist[fpB] - dist[fpA]);
 284
 285                 point newPoint
 286                 (
 287                     pA + Foam::sqrt(w / (dist[fpB] - dist[fpA]))*(pB - pA)
 288                 );
 289
 290                 Pout<< "Adapting position of vertex " << f[fp] << " on face "
 291                     << f << " from " << near.rawPoint() << " to " << newPoint
 292                     << endl;
 293
 294                 near.setPoint(newPoint);
 295             }
 296
 297             // Responsability of caller to make sure polyModifyPoint is only
 298             // called once per point. (so max only one collapse face per
 299             // edge)
 300             meshMod.setAction
 301             (
 302                 polyModifyPoint
 303                 (
 304                     f[fp],
 305                     near.rawPoint(),
 306                     false,
 307                     -1,
 308                     true
 309                 )
 310             );
 311
 312             dist[fp] = w;
 313
 314             // Step to next vertex.
 315             fpMin1 = fp;
 316             fp = f.fcIndex(fpMin1);
 317         }
 318
 319         // Step from fpA to fpB
 320         // ~~~~~~~~~~~~~~~~~~~~
 321         // (by decrementing)
 322
 323         fpMin1 = fpA;
 324         fp = f.rcIndex(fpMin1);
 325
 326         while (fp != fpB)
 327         {
 328             // See where fp sorts. Make sure it is below fpMin1!
 329             pointHit near = lineAB.nearestDist(points[f[fp]]);
 330
 331             scalar w = magSqr(near.rawPoint() - pA);
 332
 333             if (w <= dist[fpMin1])
 334             {
 335                 // Offset.
 336                 w = dist[fpMin1] + 1E-6*(dist[fpB] - dist[fpA]);
 337
 338                 point newPoint
 339                 (
 340                     pA + Foam::sqrt(w / (dist[fpB] - dist[fpA]))*(pB - pA)
 341                 );
 342
 343                 Pout<< "Adapting position of vertex " << f[fp] << " on face "
 344                     << f << " from " << near.rawPoint() << " to " << newPoint
 345                     << endl;
 346
 347                 near.setPoint(newPoint);
 348             }
 349
 350             // Responsability of caller to make sure polyModifyPoint is only
 351             // called once per point. (so max only one collapse face per
 352             // edge)
 353             meshMod.setAction
 354             (
 355                 polyModifyPoint
 356                 (
 357                     f[fp],
 358                     near.rawPoint(),
 359                     false,
 360                     -1,
 361                     true
 362                 )
 363             );
 364
 365             dist[fp] = w;
 366
 367             // Step to previous vertex.
 368             fpMin1 = fp;
 369             fp = f.rcIndex(fpMin1);
 370         }
 371
 372         dist.sort();
 373
 374         // Check that fpB sorts latest.
 375         if (dist.indices()[dist.size()-1] != fpB)
 376         {
 377             OFstream str("conflictingFace.obj");
 378             meshTools::writeOBJ(str, faceList(1, f), points);
 379
 380             FatalErrorIn("faceCollapser::setRefinement")
 381                 << "Trying to collapse face:" << faceI << " vertices:" << f
 382                 << " to edges between vertices " << f[fpA] << " and "
 383                 << f[fpB] << " but " << f[fpB] << " does not seem to be the"
 384                 << " vertex furthest away from " << f[fpA] << endl
 385                 << "Dumped conflicting face to obj file conflictingFace.obj"
 386                 << abort(FatalError);
 387         }
 388
 389
 390         // From fp to index in sort:
 391         Pout<< "Face:" << f << " fpA:" << fpA << " fpB:" << fpB << nl;
 392
 393         labelList sortedFp(f.size());
 394         forAll(dist.indices(), i)
 395         {
 396             label fp = dist.indices()[i];
 397
 398             Pout<< "   fp:" << fp << " distance:" << dist[i] << nl;
 399
 400             sortedFp[fp] = i;
 401         }
 402
 403         const labelList& fEdges = mesh_.faceEdges()[faceI];
 404
 405         // Now look up all edges in the face and see if they get extra
 406         // vertices inserted and build an edge-to-intersected-points table.
 407
 408         // Order of inserted points is in edge order (from e.start to
 409         // e.end)
 410
 411         forAll(f, fp)
 412         {
 413             label fp1 = f.fcIndex(fp);
 414
 415             // Get index in sorted list
 416             label sorted0 = sortedFp[fp];
 417             label sorted1 = sortedFp[fp1];
 418
 419             // Get indices in increasing order
 420             DynamicList<label> edgePoints(f.size());
 421
 422             // Whether edgePoints are from fp to fp1
 423             bool fpToFp1;
 424
 425             if (sorted0 < sorted1)
 426             {
 427                 fpToFp1 = true;
 428
 429                 for (label j = sorted0+1; j < sorted1; j++)
 430                 {
 431                     edgePoints.append(f[dist.indices()[j]]);
 432                 }
 433             }
 434             else
 435             {
 436                 fpToFp1 = false;
 437
 438                 for (label j = sorted1+1; j < sorted0; j++)
 439                 {
 440                     edgePoints.append(f[dist.indices()[j]]);
 441                 }
 442             }
 443
 444             if (edgePoints.size())
 445             {
 446                 edgePoints.shrink();
 447
 448                 label edgeI = findEdge(edges, fEdges, f[fp], f[fp1]);
 449
 450                 const edge& e = edges[edgeI];
 451
 452                 if (fpToFp1 == (f[fp] == e.start()))
 453                 {
 454                     splitEdges.insert(edgeI, edgePoints);
 455                 }
 456                 else
 457                 {
 458                     reverse(edgePoints);
 459                     splitEdges.insert(edgeI, edgePoints);
 460                 }
 461
 462                 // Mark all faces affected
 463                 insert(edgeFaces[edgeI], faceI, affectedFaces);
 464             }
 465         }
 466     }
 467
 468     for
 469     (
 470         Map<labelList>::const_iterator iter = splitEdges.begin();
 471         iter != splitEdges.end();
 472         ++iter
 473     )
 474     {
 475         Pout<< "Split edge:" << iter.key()
 476             << " verts:" << mesh_.edges()[iter.key()]
 477             << " in:" << nl;
 478         const labelList& edgePoints = iter();
 479
 480         forAll(edgePoints, i)
 481         {
 482             Pout<< "    " << edgePoints[i] << nl;
 483         }
 484     }
 485
 486
 487     //
 488     // Remove faces.
 489     //
 490
 491     forAll(faceLabels, i)
 492     {
 493         const label faceI = faceLabels[i];
 494
 495         meshMod.setAction(polyRemoveFace(faceI));
 496
 497         // Update list of faces we still have to modify
 498         affectedFaces.erase(faceI);
 499     }
 500
 501
 502     //
 503     // Modify faces affected (but not removed)
 504     //
 505
 506     for
 507     (
 508         labelHashSet::const_iterator iter = affectedFaces.begin();
 509         iter != affectedFaces.end();
 510         ++iter
 511     )
 512     {
 513         filterFace(splitEdges, iter.key(), meshMod);
 514     }
 515 }
 516
 517
 518 // ************************************************************************* //