applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Manipulates mesh elements.
  26
  27     Actions are:
  28         (boundary)points:
  29             - move
  30
  31         (boundary)edges:
  32             - split and move introduced point
  33
  34         (boundary)faces:
  35             - split(triangulate) and move introduced point
  36
  37         edges:
  38             - collapse
  39
  40         cells:
  41             - split into polygonal base pyramids around newly introduced mid
  42               point
  43
  44     Is a bit of a loose collection of mesh change drivers.
  45
  46 \*---------------------------------------------------------------------------*/
  47
  48 #include "argList.H"
  49 #include "objectRegistry.H"
  50 #include "foamTime.H"
  51 #include "polyMesh.H"
  52 #include "directTopoChange.H"
  53 #include "mapPolyMesh.H"
  54 #include "boundaryCutter.H"
  55 #include "cellSplitter.H"
  56 #include "edgeCollapser.H"
  57 #include "meshTools.H"
  58 #include "Pair.H"
  59
  60 using namespace Foam;
  61
  62 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  63
  64 // Locate point on patch. Returns (mesh) point label.
  65 label findPoint(const primitivePatch& pp, const point& nearPoint)
  66 {
  67     const pointField& points = pp.points();
  68     const labelList& meshPoints = pp.meshPoints();
  69
  70     // Find nearest and next nearest
  71     scalar minDistSqr = GREAT;
  72     label minI = -1;
  73
  74     scalar almostMinDistSqr = GREAT;
  75     label almostMinI = -1;
  76
  77     forAll(meshPoints, i)
  78     {
  79         label pointI = meshPoints[i];
  80
  81         scalar distSqr = magSqr(nearPoint - points[pointI]);
  82
  83         if (distSqr < minDistSqr)
  84         {
  85             almostMinDistSqr = minDistSqr;
  86             almostMinI = minI;
  87
  88             minDistSqr = distSqr;
  89             minI = pointI;
  90         }
  91         else if (distSqr < almostMinDistSqr)
  92         {
  93             almostMinDistSqr = distSqr;
  94             almostMinI = pointI;
  95         }
  96     }
  97
  98
  99     // Decide if nearPoint unique enough.
 100     Info<< "Found to point " << nearPoint << nl
 101         << "    nearest point      : " << minI
 102         << " distance " <<  Foam::sqrt(minDistSqr)
 103         << " at " << points[minI] << nl
 104         << "    next nearest point : " << almostMinI
 105         << " distance " <<  Foam::sqrt(almostMinDistSqr)
 106         << " at " << points[almostMinI] << endl;
 107
 108     if (almostMinDistSqr < 4*minDistSqr)
 109     {
 110         Info<< "Next nearest too close to nearest. Aborting" << endl;
 111
 112         return -1;
 113     }
 114     else
 115     {
 116         return minI;
 117     }
 118 }
 119
 120
 121 // Locate edge on patch. Return mesh edge label.
 122 label findEdge
 123 (
 124     const primitiveMesh& mesh,
 125     const primitivePatch& pp,
 126     const point& nearPoint
 127 )
 128 {
 129     const pointField& localPoints = pp.localPoints();
 130     const pointField& points = pp.points();
 131     const labelList& meshPoints = pp.meshPoints();
 132     const edgeList& edges = pp.edges();
 133
 134     // Find nearest and next nearest
 135     scalar minDist = GREAT;
 136     label minI = -1;
 137
 138     scalar almostMinDist = GREAT;
 139     label almostMinI = -1;
 140
 141     forAll(edges, edgeI)
 142     {
 143         const edge& e = edges[edgeI];
 144
 145         pointHit pHit(e.line(localPoints).nearestDist(nearPoint));
 146
 147         if (pHit.hit())
 148         {
 149             if (pHit.distance() < minDist)
 150             {
 151                 almostMinDist = minDist;
 152                 almostMinI = minI;
 153
 154                 minDist = pHit.distance();
 155                 minI = meshTools::findEdge
 156                 (
 157                     mesh,
 158                     meshPoints[e[0]],
 159                     meshPoints[e[1]]
 160                 );
 161             }
 162             else if (pHit.distance() < almostMinDist)
 163             {
 164                 almostMinDist = pHit.distance();
 165                 almostMinI = meshTools::findEdge
 166                 (
 167                     mesh,
 168                     meshPoints[e[0]],
 169                     meshPoints[e[1]]
 170                 );
 171             }
 172         }
 173     }
 174
 175     if (minI == -1)
 176     {
 177         Info<< "Did not find edge close to point " << nearPoint << endl;
 178
 179         return -1;
 180     }
 181
 182
 183     // Decide if nearPoint unique enough.
 184     Info<< "Found to point " << nearPoint << nl
 185         << "    nearest edge      : " << minI
 186         << " distance " << minDist << " endpoints "
 187         << mesh.edges()[minI].line(points) << nl
 188         << "    next nearest edge : " << almostMinI
 189         << " distance " << almostMinDist << " endpoints "
 190         << mesh.edges()[almostMinI].line(points) << nl
 191         << endl;
 192
 193     if (almostMinDist < 2*minDist)
 194     {
 195         Info<< "Next nearest too close to nearest. Aborting" << endl;
 196
 197         return -1;
 198     }
 199     else
 200     {
 201         return minI;
 202     }
 203 }
 204
 205
 206 // Find face on patch. Return mesh face label.
 207 label findFace
 208 (
 209     const primitiveMesh& mesh,
 210     const primitivePatch& pp,
 211     const point& nearPoint
 212 )
 213 {
 214     const pointField& points = pp.points();
 215
 216     // Find nearest and next nearest
 217     scalar minDist = GREAT;
 218     label minI = -1;
 219
 220     scalar almostMinDist = GREAT;
 221     label almostMinI = -1;
 222
 223     forAll(pp, patchFaceI)
 224     {
 225         pointHit pHit(pp[patchFaceI].nearestPoint(nearPoint, points));
 226
 227         if (pHit.hit())
 228         {
 229             if (pHit.distance() < minDist)
 230             {
 231                 almostMinDist = minDist;
 232                 almostMinI = minI;
 233
 234                 minDist = pHit.distance();
 235                 minI = patchFaceI + mesh.nInternalFaces();
 236             }
 237             else if (pHit.distance() < almostMinDist)
 238             {
 239                 almostMinDist = pHit.distance();
 240                 almostMinI = patchFaceI + mesh.nInternalFaces();
 241             }
 242         }
 243     }
 244
 245     if (minI == -1)
 246     {
 247         Info<< "Did not find face close to point " << nearPoint << endl;
 248
 249         return -1;
 250     }
 251
 252
 253     // Decide if nearPoint unique enough.
 254     Info<< "Found to point " << nearPoint << nl
 255         << "    nearest face      : " << minI
 256         << " distance " << minDist
 257         << " to face centre " << mesh.faceCentres()[minI] << nl
 258         << "    next nearest face : " << almostMinI
 259         << " distance " << almostMinDist
 260         << " to face centre " << mesh.faceCentres()[almostMinI] << nl
 261         << endl;
 262
 263     if (almostMinDist < 2*minDist)
 264     {
 265         Info<< "Next nearest too close to nearest. Aborting" << endl;
 266
 267         return -1;
 268     }
 269     else
 270     {
 271         return minI;
 272     }
 273 }
 274
 275
 276 // Find cell with cell centre close to given point.
 277 label findCell(const primitiveMesh& mesh, const point& nearPoint)
 278 {
 279     label cellI = mesh.findCell(nearPoint);
 280
 281     if (cellI != -1)
 282     {
 283         scalar distToCcSqr = magSqr(nearPoint - mesh.cellCentres()[cellI]);
 284
 285         const labelList& cPoints = mesh.cellPoints()[cellI];
 286
 287         label minI = -1;
 288         scalar minDistSqr = GREAT;
 289
 290         forAll(cPoints, i)
 291         {
 292             label pointI = cPoints[i];
 293
 294             scalar distSqr = magSqr(nearPoint - mesh.points()[pointI]);
 295
 296             if (distSqr < minDistSqr)
 297             {
 298                 minDistSqr = distSqr;
 299                 minI = pointI;
 300             }
 301         }
 302
 303         // Decide if nearPoint unique enough.
 304         Info<< "Found to point " << nearPoint << nl
 305             << "    nearest cell       : " << cellI
 306             << " distance " << Foam::sqrt(distToCcSqr)
 307             << " to cell centre " << mesh.cellCentres()[cellI] << nl
 308             << "    nearest mesh point : " << minI
 309             << " distance " << Foam::sqrt(minDistSqr)
 310             << " to " << mesh.points()[minI] << nl
 311             << endl;
 312
 313         if (minDistSqr < 4*distToCcSqr)
 314         {
 315             Info<< "Mesh point too close to nearest cell centre. Aborting"
 316                 << endl;
 317
 318             cellI = -1;
 319         }
 320     }
 321
 322     return cellI;
 323 }
 324
 325
 326
 327 // Main program:
 328
 329 int main(int argc, char *argv[])
 330 {
 331     argList::validOptions.insert("overwrite", "");
 332
 333 #   include "setRootCase.H"
 334 #   include "createTime.H"
 335     runTime.functionObjects().off();
 336 #   include "createPolyMesh.H"
 337     const word oldInstance = mesh.pointsInstance();
 338
 339     bool overwrite = args.optionFound("overwrite");
 340
 341     Info<< "Reading modifyMeshDict\n" << endl;
 342
 343     IOdictionary dict
 344     (
 345         IOobject
 346         (
 347             "modifyMeshDict",
 348             runTime.system(),
 349             mesh,
 350             IOobject::MUST_READ,
 351             IOobject::NO_WRITE
 352         )
 353     );
 354
 355     // Read all from the dictionary.
 356     List<Pair<point> > pointsToMove(dict.lookup("pointsToMove"));
 357     List<Pair<point> > edgesToSplit(dict.lookup("edgesToSplit"));
 358     List<Pair<point> > facesToTriangulate
 359     (
 360         dict.lookup("facesToTriangulate")
 361     );
 362
 363     bool cutBoundary =
 364     (
 365         pointsToMove.size()
 366      || edgesToSplit.size()
 367      || facesToTriangulate.size()
 368     );
 369
 370     List<Pair<point> > edgesToCollapse(dict.lookup("edgesToCollapse"));
 371
 372     bool collapseEdge = edgesToCollapse.size();
 373
 374     List<Pair<point> > cellsToPyramidise(dict.lookup("cellsToSplit"));
 375
 376     bool cellsToSplit = cellsToPyramidise.size();
 377
 378     //List<Tuple<pointField,point> >
 379     //  cellsToCreate(dict.lookup("cellsToCreate"));
 380
 381     Info<< "Read from " << dict.name() << nl
 382         << "  Boundary cutting module:" << nl
 383         << "    points to move      :" << pointsToMove.size() << nl
 384         << "    edges to split      :" << edgesToSplit.size() << nl
 385         << "    faces to triangulate:" << facesToTriangulate.size() << nl
 386         << "  Cell splitting module:" << nl
 387         << "    cells to split      :" << cellsToPyramidise.size() << nl
 388         << "  Edge collapsing module:" << nl
 389         << "    edges to collapse   :" << edgesToCollapse.size() << nl
 390         //<< "    cells to create     :" << cellsToCreate.size() << nl
 391         << endl;
 392
 393     if
 394     (
 395         (cutBoundary && collapseEdge)
 396      || (cutBoundary && cellsToSplit)
 397      || (collapseEdge && cellsToSplit)
 398     )
 399     {
 400         FatalErrorIn(args.executable())
 401             << "Used more than one mesh modifying module "
 402             << "(boundary cutting, cell splitting, edge collapsing)" << nl
 403             << "Please do them in separate passes." << exit(FatalError);
 404     }
 405
 406
 407
 408     // Get calculating engine for all of outside
 409     const SubList<face> outsideFaces
 410     (
 411         mesh.faces(),
 412         mesh.nFaces() - mesh.nInternalFaces(),
 413         mesh.nInternalFaces()
 414     );
 415
 416     primitivePatch allBoundary(outsideFaces, mesh.points());
 417
 418
 419     // Look up mesh labels and convert to input for boundaryCutter.
 420
 421     bool validInputs = true;
 422
 423
 424     Info<< nl << "Looking up points to move ..." << nl << endl;
 425     Map<point> pointToPos(pointsToMove.size());
 426     forAll(pointsToMove, i)
 427     {
 428         const Pair<point>& pts = pointsToMove[i];
 429
 430         label pointI = findPoint(allBoundary, pts.first());
 431
 432         if (pointI == -1 || !pointToPos.insert(pointI, pts.second()))
 433         {
 434             Info<< "Could not insert mesh point " << pointI
 435                 << " for input point " << pts.first() << nl
 436                 << "Perhaps the point is already marked for moving?" << endl;
 437             validInputs = false;
 438         }
 439     }
 440
 441
 442     Info<< nl << "Looking up edges to split ..." << nl << endl;
 443     Map<List<point> > edgeToCuts(edgesToSplit.size());
 444     forAll(edgesToSplit, i)
 445     {
 446         const Pair<point>& pts = edgesToSplit[i];
 447
 448         label edgeI = findEdge(mesh, allBoundary, pts.first());
 449
 450         if
 451         (
 452             edgeI == -1
 453         || !edgeToCuts.insert(edgeI, List<point>(1, pts.second()))
 454         )
 455         {
 456             Info<< "Could not insert mesh edge " << edgeI
 457                 << " for input point " << pts.first() << nl
 458                 << "Perhaps the edge is already marked for cutting?" << endl;
 459
 460             validInputs = false;
 461         }
 462     }
 463
 464
 465     Info<< nl << "Looking up faces to triangulate ..." << nl << endl;
 466     Map<point> faceToDecompose(facesToTriangulate.size());
 467     forAll(facesToTriangulate, i)
 468     {
 469         const Pair<point>& pts = facesToTriangulate[i];
 470
 471         label faceI = findFace(mesh, allBoundary, pts.first());
 472
 473         if (faceI == -1 || !faceToDecompose.insert(faceI, pts.second()))
 474         {
 475             Info<< "Could not insert mesh face " << faceI
 476                 << " for input point " << pts.first() << nl
 477                 << "Perhaps the face is already marked for splitting?" << endl;
 478
 479             validInputs = false;
 480         }
 481     }
 482
 483
 484
 485     Info<< nl << "Looking up cells to convert to pyramids around"
 486         << " cell centre ..." << nl << endl;
 487     Map<point> cellToPyrCentre(cellsToPyramidise.size());
 488     forAll(cellsToPyramidise, i)
 489     {
 490         const Pair<point>& pts = cellsToPyramidise[i];
 491
 492         label cellI = findCell(mesh, pts.first());
 493
 494         if (cellI == -1 || !cellToPyrCentre.insert(cellI, pts.second()))
 495         {
 496             Info<< "Could not insert mesh cell " << cellI
 497                 << " for input point " << pts.first() << nl
 498                 << "Perhaps the cell is already marked for splitting?" << endl;
 499
 500             validInputs = false;
 501         }
 502     }
 503
 504
 505     Info<< nl << "Looking up edges to collapse ..." << nl << endl;
 506     Map<point> edgeToPos(edgesToCollapse.size());
 507     forAll(edgesToCollapse, i)
 508     {
 509         const Pair<point>& pts = edgesToCollapse[i];
 510
 511         label edgeI = findEdge(mesh, allBoundary, pts.first());
 512
 513         if (edgeI == -1 || !edgeToPos.insert(edgeI, pts.second()))
 514         {
 515             Info<< "Could not insert mesh edge " << edgeI
 516                 << " for input point " << pts.first() << nl
 517                 << "Perhaps the edge is already marked for collaping?" << endl;
 518
 519             validInputs = false;
 520         }
 521     }
 522
 523
 524
 525     if (!validInputs)
 526     {
 527         Info<< nl << "There was a problem in one of the inputs in the"
 528             << " dictionary. Not modifying mesh." << endl;
 529     }
 530     else if (cellToPyrCentre.size())
 531     {
 532         Info<< nl << "All input cells located. Modifying mesh." << endl;
 533
 534         // Mesh change engine
 535         cellSplitter cutter(mesh);
 536
 537         // Topo change container
 538         directTopoChange meshMod(mesh);
 539
 540         // Insert commands into meshMod
 541         cutter.setRefinement(cellToPyrCentre, meshMod);
 542
 543         // Do changes
 544         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 545
 546         if (morphMap().hasMotionPoints())
 547         {
 548             mesh.movePoints(morphMap().preMotionPoints());
 549         }
 550
 551         cutter.updateMesh(morphMap());
 552
 553         if (!overwrite)
 554         {
 555             runTime++;
 556         }
 557         else
 558         {
 559             mesh.setInstance(oldInstance);
 560         }
 561
 562         // Write resulting mesh
 563         Info << "Writing modified mesh to time " << runTime.timeName() << endl;
 564         mesh.write();
 565     }
 566     else if (edgeToPos.size())
 567     {
 568         Info<< nl << "All input edges located. Modifying mesh." << endl;
 569
 570         // Mesh change engine
 571         edgeCollapser cutter(mesh);
 572
 573         pointField newPoints(mesh.points());
 574
 575         // Get new positions and construct collapse network
 576         forAllConstIter(Map<point>, edgeToPos, iter)
 577         {
 578             label edgeI = iter.key();
 579             const edge& e = mesh.edges()[edgeI];
 580
 581             cutter.collapseEdge(edgeI, e[0]);
 582             newPoints[e[0]] = iter();
 583         }
 584
 585         // Move master point to destination.
 586         mesh.movePoints(newPoints);
 587
 588         // Topo change container
 589         directTopoChange meshMod(mesh);
 590
 591         // Insert
 592         cutter.setRefinement(meshMod);
 593
 594         // Do changes
 595         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 596
 597         if (morphMap().hasMotionPoints())
 598         {
 599             mesh.movePoints(morphMap().preMotionPoints());
 600         }
 601
 602         // Not implemented yet:
 603         //cutter.updateMesh(morphMap());
 604
 605
 606         if (!overwrite)
 607         {
 608             runTime++;
 609         }
 610         else
 611         {
 612             mesh.setInstance(oldInstance);
 613         }
 614
 615         // Write resulting mesh
 616         Info << "Writing modified mesh to time " << runTime.timeName() << endl;
 617         mesh.write();
 618     }
 619     else
 620     {
 621         Info<< nl << "All input points located. Modifying mesh." << endl;
 622
 623         // Mesh change engine
 624         boundaryCutter cutter(mesh);
 625
 626         // Topo change container
 627         directTopoChange meshMod(mesh);
 628
 629         // Insert commands into meshMod
 630         cutter.setRefinement
 631         (
 632             pointToPos,
 633             edgeToCuts,
 634             Map<labelPair>(0),  // Faces to split diagonally
 635             faceToDecompose,    // Faces to triangulate
 636             meshMod
 637         );
 638
 639         // Do changes
 640         autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
 641
 642         if (morphMap().hasMotionPoints())
 643         {
 644             mesh.movePoints(morphMap().preMotionPoints());
 645         }
 646
 647         cutter.updateMesh(morphMap());
 648
 649         if (!overwrite)
 650         {
 651             runTime++;
 652         }
 653         else
 654         {
 655             mesh.setInstance(oldInstance);
 656         }
 657
 658         // Write resulting mesh
 659         Info << "Writing modified mesh to time " << runTime.timeName() << endl;
 660         mesh.write();
 661     }
 662
 663
 664     Info << nl << "End" << endl;
 665
 666     return 0;
 667 }
 668
 669
 670 // ************************************************************************* //