src/mesh/cfMesh/utilities/faceDecomposition/decomposeFaces.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | cfMesh: A library for mesh generation
   4    \\    /   O peration     |
   5     \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
   6      \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of cfMesh.
  10
  11     cfMesh 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     cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28
  29 #include "decomposeFaces.H"
  30 #include "boolList.H"
  31
  32 # ifdef USE_OMP
  33 #include <omp.h>
  34 # endif
  35
  36 //#define DEBUGDec
  37
  38 # ifdef DEBUGDec
  39 #include "polyMeshGenAddressing.H"
  40 # endif
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46
  47 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  48 //- Constructor
  49 decomposeFaces::decomposeFaces(polyMeshGen& mesh)
  50 :
  51     mesh_(mesh),
  52     newFacesForFace_(),
  53     done_(false)
  54 {}
  55
  56 //- Destructor
  57 decomposeFaces::~decomposeFaces()
  58 {}
  59
  60 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  61
  62 void decomposeFaces::decomposeMeshFaces(const boolList& decomposeFace)
  63 {
  64     done_ = false;
  65
  66     newFacesForFace_.setSize(mesh_.faces().size());
  67     forAll(newFacesForFace_, fI)
  68         newFacesForFace_.setRowSize(fI, 0);
  69
  70     const label nIntFaces = mesh_.nInternalFaces();
  71     label nFaces(0), nPoints = mesh_.points().size();
  72     point p;
  73
  74     pointFieldPMG& points = mesh_.points();
  75     forAll(decomposeFace, fI)
  76         if( decomposeFace[fI] )
  77             ++nFaces;
  78
  79     points.setSize(nPoints + nFaces);
  80
  81     polyMeshGenModifier meshModifier(mesh_);
  82     faceListPMG& faces = meshModifier.facesAccess();
  83
  84     if( decomposeFace.size() != faces.size() )
  85         FatalErrorIn
  86         (
  87             "void decomposeFaces::decomposeMeshFaces(const boolList&)"
  88         ) << "Incorrect size of the decomposeFace list!" << abort(FatalError);
  89
  90     nFaces = 0;
  91     VRWGraph newFaces;
  92
  93     //- decompose internal faces
  94     for(label faceI=0;faceI<nIntFaces;++faceI)
  95     {
  96         const face& f = faces[faceI];
  97
  98         if( decomposeFace[faceI] )
  99         {
 100             # ifdef DEBUGDec
 101             Info << "Decomposing internal face " << faceI << " with nodes "
 102                 << f << endl;
 103             # endif
 104
 105             FixedList<label, 3> newF;
 106
 107             forAll(f, pI)
 108             {
 109                 newF[0] = f[pI];
 110                 newF[1] = f.nextLabel(pI);
 111                 newF[2] = nPoints;
 112
 113                 # ifdef DEBUGDec
 114                 Info << "Storing face " << newF << " with label "
 115                     << nFaces << endl;
 116                 # endif
 117
 118                 newFaces.appendList(newF);
 119                 newFacesForFace_.append(faceI, nFaces++);
 120             }
 121
 122             p = f.centre(points);
 123             points[nPoints] = p;
 124             ++nPoints;
 125         }
 126         else
 127         {
 128             # ifdef DEBUGDec
 129             Info << "Storing internal face " << faceI << " with nodes "
 130                 << f << " as new face " << faceI << endl;
 131             # endif
 132
 133             newFaces.appendList(f);
 134             newFacesForFace_.append(faceI, nFaces++);
 135         }
 136     }
 137
 138     //- decompose boundary faces
 139     PtrList<boundaryPatch>& boundaries = meshModifier.boundariesAccess();
 140     forAll(boundaries, patchI)
 141     {
 142         const label start = boundaries[patchI].patchStart();
 143         const label end = start + boundaries[patchI].patchSize();
 144
 145         boundaries[patchI].patchStart() = nFaces;
 146
 147         for(label bfI=start;bfI<end;++bfI)
 148         {
 149             const face& f = faces[bfI];
 150             if( decomposeFace[bfI] )
 151             {
 152                 # ifdef DEBUGDec
 153                 Info << "Decomposing boundary face " << bfI
 154                     << " with nodes " << f << endl;
 155                 # endif
 156
 157                 FixedList<label, 3> newF;
 158
 159                 forAll(f, pI)
 160                 {
 161                     newF[0] = f[pI];
 162                     newF[1] = f.nextLabel(pI);
 163                     newF[2] = nPoints;
 164
 165                     # ifdef DEBUGDec
 166                     Info << "Storing face " << newF << " with label "
 167                         << nFaces << endl;
 168                     # endif
 169
 170                     newFaces.appendList(newF);
 171                     newFacesForFace_.append(bfI, nFaces++);
 172                 }
 173
 174                 p = f.centre(points);
 175                 points[nPoints++] = p;
 176             }
 177             else
 178             {
 179                 # ifdef DEBUGDec
 180                 Info << "Storing boundary face " << bfI << " in patch"
 181                     << patchI << " as new face " << bfI << endl;
 182                 # endif
 183
 184                 newFaces.appendList(f);
 185                 newFacesForFace_.append(bfI, nFaces++);
 186             }
 187         }
 188
 189         boundaries[patchI].patchSize() =
 190             nFaces - boundaries[patchI].patchStart();
 191     }
 192
 193     //- decompose processor faces
 194     if( Pstream::parRun() )
 195     {
 196         PtrList<processorBoundaryPatch>& procBoundaries =
 197             meshModifier.procBoundariesAccess();
 198
 199         forAll(procBoundaries, patchI)
 200         {
 201             const label start = procBoundaries[patchI].patchStart();
 202             const label end = start + procBoundaries[patchI].patchSize();
 203
 204             const bool own = procBoundaries[patchI].owner();
 205
 206             procBoundaries[patchI].patchStart() = nFaces;
 207
 208             for(label bfI=start;bfI<end;++bfI)
 209             {
 210                 face f;
 211                 if( own )
 212                 {
 213                     f = faces[bfI];
 214                 }
 215                 else
 216                 {
 217                     f = faces[bfI].reverseFace();
 218                 }
 219
 220                 if( decomposeFace[bfI] )
 221                 {
 222                     # ifdef DEBUGDec
 223                     Info << "Decomposing processor boundary face " << bfI
 224                         << " with nodes " << f << endl;
 225                     # endif
 226
 227                     face newF(3);
 228
 229                     forAll(f, pI)
 230                     {
 231                         newF[0] = f[pI];
 232                         newF[1] = f.nextLabel(pI);
 233                         newF[2] = nPoints;
 234
 235                         # ifdef DEBUGDec
 236                         Info << "Storing face " << newF << " with label "
 237                             << nFaces << endl;
 238                         # endif
 239
 240                         if( own )
 241                         {
 242                             newFaces.appendList(newF);
 243                         }
 244                         else
 245                         {
 246                             newFaces.appendList(newF.reverseFace());
 247                         }
 248                         newFacesForFace_.append(bfI, nFaces++);
 249                     }
 250
 251                     p = f.centre(points);
 252                     points[nPoints++] = p;
 253                 }
 254                 else
 255                 {
 256                     # ifdef DEBUGDec
 257                     Info << "Storing boundary face " << bfI << " in patch"
 258                         << patchI << " as new face " << bfI << endl;
 259                     # endif
 260
 261                     if( own )
 262                     {
 263                         newFaces.appendList(f);
 264                     }
 265                     else
 266                     {
 267                         newFaces.appendList(f.reverseFace());
 268                     }
 269                     newFacesForFace_.append(bfI, nFaces++);
 270                 }
 271             }
 272
 273             procBoundaries[patchI].patchSize() =
 274                 nFaces - procBoundaries[patchI].patchStart();
 275         }
 276     }
 277
 278     //- store the faces back into their list
 279     faces.setSize(nFaces);
 280     forAll(faces, faceI)
 281     {
 282         face& f = faces[faceI];
 283         f.setSize(newFaces.sizeOfRow(faceI));
 284
 285         forAll(f, pI)
 286             f[pI] = newFaces(faceI, pI);
 287     }
 288     newFaces.setSize(0);
 289
 290     //- update subsets
 291     mesh_.updateFaceSubsets(newFacesForFace_);
 292
 293     //- change the mesh
 294     cellListPMG& cells = meshModifier.cellsAccess();
 295
 296     # ifdef USE_OMP
 297     # pragma omp parallel for schedule(dynamic, 40)
 298     # endif
 299     forAll(cells, cellI)
 300     {
 301         cell& c = cells[cellI];
 302
 303         DynList<label> newC;
 304
 305         forAll(c, fJ)
 306         {
 307             const label faceI = c[fJ];
 308             forAllRow(newFacesForFace_, faceI, nfI)
 309             newC.append(newFacesForFace_(faceI, nfI));
 310         }
 311
 312         # ifdef DEBUGDec
 313         Info << "Cell " << cellI << " with faces " << c
 314             << " is changed into " << newC << endl;
 315         # endif
 316
 317         c.setSize(newC.size());
 318         forAll(newC, fJ)
 319         c[fJ] = newC[fJ];
 320     }
 321
 322     meshModifier.clearAll();
 323
 324     done_ = true;
 325
 326     # ifdef DEBUGDec
 327     Info << "New cells are " << cells << endl;
 328     mesh_.write();
 329     # endif
 330 }
 331
 332 void decomposeFaces::decomposeConcaveInternalFaces
 333 (
 334     const boolList& concaveVertex
 335 )
 336 {
 337     if( Pstream::parRun() )
 338     {
 339         FatalErrorIn
 340         (
 341             "void decomposeFaces::decomposeConcaveInternalFaces"
 342             "(const boolList& concaveVertex)"
 343         ) << "This procedure is not parallelised!" << exit(FatalError);
 344     }
 345
 346     done_ = false;
 347
 348     newFacesForFace_.setSize(mesh_.faces().size());
 349     forAll(newFacesForFace_, fI)
 350         newFacesForFace_.setRowSize(fI, 0);
 351
 352     const label nIntFaces = mesh_.nInternalFaces();
 353
 354     polyMeshGenModifier meshModifier(mesh_);
 355     pointFieldPMG& points = meshModifier.pointsAccess();
 356     faceListPMG& faces = meshModifier.facesAccess();
 357
 358     if( concaveVertex.size() != mesh_.points().size() )
 359         FatalErrorIn
 360         (
 361             "void decomposeFaces::decomposeMeshFaces(const boolList&)"
 362         ) << "Incorrect size of the concaveVertex list!" << abort(FatalError);
 363
 364     VRWGraph newFaces;
 365     DynList<label> newF;
 366     newF.setSize(3);
 367
 368     # ifdef DEBUGDec
 369     const label id = mesh_.addFaceSubset("decomposedFaces");
 370     # endif
 371
 372     //- decompose internal faces
 373     for(register label faceI=0;faceI<nIntFaces;++faceI)
 374     {
 375         const face& f = faces[faceI];
 376
 377         DynList<label> concavePos;
 378         forAll(f, pI)
 379             if( concaveVertex[f[pI]] )
 380             {
 381                 concavePos.append(pI);
 382             }
 383
 384         if( concavePos.size() == 1 )
 385         {
 386             # ifdef DEBUGDec
 387             Info << "1. Decomposing internal face " << faceI << " with nodes "
 388                 << f << endl;
 389             mesh_.addFaceToSubset(id, faceI);
 390             # endif
 391
 392             newF[0] = f[concavePos[0]];
 393
 394             for(label pI=1;pI<(f.size()-1);++pI)
 395             {
 396                 const label pJ = (concavePos[0] + pI) % f.size();
 397                 newF[1] = f[pJ];
 398                 newF[2] = f.nextLabel(pJ);
 399
 400                 # ifdef DEBUGDec
 401                 Info << "Storing face " << newF << " with label "
 402                     << newFaces.size() << endl;
 403                 # endif
 404
 405                 newFacesForFace_.append(faceI, newFaces.size());
 406                 newFaces.appendList(newF);
 407             }
 408         }
 409         else if( concavePos.size() > 1 )
 410         {
 411             # ifdef DEBUGDec
 412             Info << "2. Decomposing internal face " << faceI << " with nodes "
 413                 << f << endl;
 414             mesh_.addFaceToSubset(id, faceI);
 415             # endif
 416
 417             newF[0] = points.size();
 418             forAll(f, pI)
 419             {
 420                 newF[1] = f[pI];
 421                 newF[2] = f.nextLabel(pI);
 422
 423                 # ifdef DEBUGDec
 424                 Info << "2. Storing face " << newF << " with label "
 425                     << newFaces.size() << endl;
 426                 # endif
 427
 428                 newFacesForFace_.append(faceI, newFaces.size());
 429                 newFaces.appendList(newF);
 430             }
 431
 432             const point fCent = f.centre(points);
 433             points.append(fCent);
 434         }
 435         else
 436         {
 437             # ifdef DEBUGDec
 438             Info << "Storing internal face " << faceI << " with nodes "
 439                 << f << " as new face " << newFaces.size() << endl;
 440             # endif
 441
 442             newFacesForFace_.append(faceI, newFaces.size());
 443             newFaces.appendList(f);
 444         }
 445     }
 446
 447     PtrList<boundaryPatch>& boundaries = meshModifier.boundariesAccess();
 448     forAll(boundaries, patchI)
 449     {
 450         const label start = boundaries[patchI].patchStart();
 451         const label end = start + boundaries[patchI].patchSize();
 452
 453         //- set new patch start
 454         boundaries[patchI].patchStart() = newFaces.size();
 455
 456         //- store faces into newFaces
 457         for(label bfI=start;bfI<end;++bfI)
 458         {
 459             newFacesForFace_.append(bfI, newFaces.size());
 460             newFaces.appendList(faces[bfI]);
 461         }
 462     }
 463
 464     //- copy new faces into the faceListPMG
 465     faces.setSize(newFaces.size());
 466     forAll(newFaces, faceI)
 467     {
 468         faces[faceI].setSize(newFaces.sizeOfRow(faceI));
 469
 470         forAllRow(newFaces, faceI, pI)
 471             faces[faceI][pI] = newFaces(faceI, pI);
 472     }
 473
 474     newFaces.setSize(0);
 475
 476     //- update cells
 477     cellListPMG& cells = meshModifier.cellsAccess();
 478
 479     # ifdef USE_OMP
 480     # pragma omp parallel for schedule(dynamic, 40)
 481     # endif
 482     forAll(cells, cellI)
 483     {
 484         cell& c = cells[cellI];
 485
 486         DynList<label, 24> newC;
 487
 488         forAll(c, fJ)
 489         {
 490             const label faceI = c[fJ];
 491             forAllRow(newFacesForFace_, faceI, nfI)
 492             newC.append(newFacesForFace_(faceI, nfI));
 493         }
 494
 495         # ifdef DEBUGDec
 496         Info << "Cell " << cellI << " with faces " << c
 497             << " is changed into " << newC << endl;
 498         # endif
 499
 500         c.setSize(newC.size());
 501         forAll(newC, fJ)
 502             c[fJ] = newC[fJ];
 503     }
 504
 505     meshModifier.clearAll();
 506
 507     //- update subsets
 508     mesh_.updateFaceSubsets(newFacesForFace_);
 509
 510     # ifdef DEBUGDec
 511     Info << "New cells are " << cells << endl;
 512     mesh_.write();
 513     ::exit(1);
 514     # endif
 515
 516     meshModifier.removeUnusedVertices();
 517
 518     done_ = true;
 519 }
 520
 521 const VRWGraph& decomposeFaces::newFacesForFace() const
 522 {
 523     if( !done_ )
 524         WarningIn
 525         (
 526             "const VRWGraph& decomposeFaces::newFacesForFace() const"
 527         ) << "Decomposition is not yet performed!" << endl;
 528
 529     return newFacesForFace_;
 530 }
 531
 532 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
 533
 534 } // End namespace Foam
 535
 536 // ************************************************************************* //