src/mesh/cfMesh/utilities/surfaceTools/surfaceMorpherCells/surfaceMorpherCellsMorphInternalFaces.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 #include "surfaceMorpherCells.H"
  29 #include "demandDrivenData.H"
  30 #include "helperFunctions.H"
  31 #include "HashSet.H"
  32
  33 //#define DEBUGMorph
  34
  35 #ifdef DEBUGMorph
  36 #include "polyMeshGenAddressing.H"
  37 #endif
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 namespace Foam
  42 {
  43
  44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  45
  46 void surfaceMorpherCells::findBoundaryVertices()
  47 {
  48     const faceListPMG& faces = mesh_.faces();
  49
  50     boundaryVertex_.setSize(mesh_.points().size());
  51     boundaryVertex_ = false;
  52
  53     const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
  54     forAll(boundaries, patchI)
  55     {
  56         const label start = boundaries[patchI].patchStart();
  57         const label end = start + boundaries[patchI].patchSize();
  58
  59         for(label faceI=start;faceI<end;++faceI)
  60         {
  61             const face& f = faces[faceI];
  62
  63             forAll(f, pI)
  64             {
  65                 # ifdef DEBUGMorph
  66                 if( (f[pI] >= boundaryVertex_.size()) || (f[pI] < 0) )
  67                 {
  68                     Info << f << endl;
  69                     FatalError << "Wrong label " << f[pI] << " in face "
  70                         << faceI << abort(FatalError);
  71                 }
  72                 # endif
  73
  74                 boundaryVertex_[f[pI]] = true;
  75             }
  76         }
  77     }
  78
  79     if( Pstream::parRun() )
  80     {
  81         const PtrList<processorBoundaryPatch>& procBoundaries =
  82             mesh_.procBoundaries();
  83
  84         bool changed;
  85         do
  86         {
  87             changed = false;
  88
  89             //- send data about boundary vertices to other processors
  90             forAll(procBoundaries, patchI)
  91             {
  92                 const label start = procBoundaries[patchI].patchStart();
  93                 const label end = start + procBoundaries[patchI].patchSize();
  94
  95                 //- create information about bnd nodes which must be exchanged
  96                 //- with other processors
  97                 labelHashSet addToSend;
  98                 labelLongList dts;
  99                 for(label faceI=start;faceI<end;++faceI)
 100                 {
 101                     const face& f = faces[faceI];
 102
 103                     forAll(f, pI)
 104                         if( boundaryVertex_[f[pI]] && !addToSend.found(f[pI]) )
 105                         {
 106                             addToSend.insert(f[pI]);
 107                             dts.append(faceI-start);
 108                             dts.append((f.size()-pI)%f.size());
 109                         }
 110                 }
 111
 112                 labelList bndVertsToSend(dts.size());
 113                 forAll(dts, i)
 114                     bndVertsToSend[i] = dts[i];
 115
 116                 //- send the list of other processor
 117                 OPstream toOtherProc
 118                 (
 119                     Pstream::blocking,
 120                     procBoundaries[patchI].neiProcNo(),
 121                     bndVertsToSend.byteSize()
 122                 );
 123                 toOtherProc << bndVertsToSend;
 124             }
 125
 126             //- receive boundary vertices from other processor and
 127             //- continue sending and receiving as long as it makes some change
 128             forAll(procBoundaries, patchI)
 129             {
 130                 labelList receivedBndNodes;
 131                 IPstream fromOtherProc
 132                 (
 133                     Pstream::blocking,
 134                     procBoundaries[patchI].neiProcNo()
 135                 );
 136                 fromOtherProc >> receivedBndNodes;
 137                 const label start = procBoundaries[patchI].patchStart();
 138
 139                 label entryI(0);
 140                 while( entryI < receivedBndNodes.size() )
 141                 {
 142                     const label fI = receivedBndNodes[entryI++];
 143                     const label pI = receivedBndNodes[entryI++];
 144
 145                     const face& f = faces[start+fI];
 146
 147                     if( !boundaryVertex_[f[pI]] )
 148                     {
 149                         boundaryVertex_[f[pI]] = true;
 150                         changed = true;
 151                     }
 152                 }
 153             }
 154
 155             reduce(changed, maxOp<bool>());
 156         } while( changed );
 157     }
 158 }
 159
 160 void surfaceMorpherCells::findBoundaryCells()
 161 {
 162     const labelList& owner = mesh_.owner();
 163
 164     cellFlags_.setSize(mesh_.cells().size());
 165     cellFlags_ = NONE;
 166
 167     const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
 168
 169     forAll(boundaries, patchI)
 170     {
 171         const label start = boundaries[patchI].patchStart();
 172         const label end = start + boundaries[patchI].patchSize();
 173
 174         for(label faceI=start;faceI<end;++faceI)
 175             cellFlags_[owner[faceI]] = BOUNDARY;
 176     }
 177 }
 178
 179 bool surfaceMorpherCells::morphInternalFaces()
 180 {
 181     Info << "Morphing internal faces" << endl;
 182
 183     # ifdef DEBUGMorph
 184     Serr << "Zip check 1" << endl;
 185     labelHashSet zipCellsBefore;
 186     mesh_.addressingData().checkCellsZipUp(true, &zipCellsBefore);
 187     if( zipCellsBefore.size() )
 188     {
 189         Serr << Pstream::myProcNo() << "Open cells found!" << endl;
 190         Serr << "Cells " << zipCellsBefore << " are not zipped!!" << endl;
 191     }
 192     mesh_.clearAddressingData();
 193     # endif
 194
 195     bool changed(false);
 196
 197     newBoundaryFaces_.setSize(0);
 198     newBoundaryOwners_.setSize(0);
 199     newBoundaryPatches_.setSize(0);
 200
 201     const label nIntFaces = mesh_.nInternalFaces();
 202     const faceListPMG& faces = mesh_.faces();
 203     const labelList& owner = mesh_.owner();
 204     const labelList& neighbour = mesh_.neighbour();
 205
 206     //- copy boundary faces
 207     const PtrList<boundaryPatch>& boundaries = mesh_.boundaries();
 208
 209     forAll(boundaries, patchI)
 210     {
 211         const label start = boundaries[patchI].patchStart();
 212         const label end = start + boundaries[patchI].patchSize();
 213
 214         for(label faceI=start;faceI<end;++faceI)
 215         {
 216             newBoundaryFaces_.appendList(faces[faceI]);
 217             newBoundaryOwners_.append(owner[faceI]);
 218             newBoundaryPatches_.append(0);
 219         }
 220     }
 221
 222     //- start morphing internal faces
 223     for(label faceI=0;faceI<nIntFaces;++faceI)
 224     {
 225         # ifdef DEBUGMorph
 226         Info << "Morphing internal face " << faceI << endl;
 227         # endif
 228
 229         const face& f = faces[faceI];
 230
 231         DynList<bool> removeFaceVertex(f.size(), false);
 232
 233         face newF(f.size());
 234         label i(0);
 235
 236         forAll(f, pI)
 237             if(
 238                 boundaryVertex_[f.prevLabel(pI)] &&
 239                 boundaryVertex_[f[pI]] &&
 240                 boundaryVertex_[f.nextLabel(pI)]
 241             )
 242             {
 243                 removeFaceVertex[pI] = true;
 244
 245                 # ifdef DEBUGMorph
 246                 Info << "Removing vertex " << f[pI] << " from face "
 247                     << f << endl;
 248                 # endif
 249             }
 250             else
 251             {
 252                 newF[i++] = f[pI];
 253             }
 254
 255         if( i < f.size() )
 256         {
 257             changed = true;
 258
 259             //- store shrinked face
 260             newF.setSize(i);
 261
 262             # ifdef DEBUGMorph
 263             Info << "Removing edges " << removeEdge << endl;
 264             Info << "Face label " << faceI << " owner " << owner[faceI]
 265                 << " neighbour " << neighbour[faceI] << endl;
 266             Info << "Internal face " << f << " has been shrinked into "
 267                 << newF << endl;
 268             # endif
 269
 270             //- create new boundary faces from the removed part
 271             label mat(1);
 272             DynList<direction> nodeMaterial(f.size(), direction(0));
 273             DynList<DynList<edge>, 2> edgeMats;
 274             forAll(nodeMaterial, nI)
 275                 if( !nodeMaterial[nI] && removeFaceVertex[nI] )
 276                 {
 277                     edgeMats.append(DynList<edge>());
 278                     DynList<label> front;
 279                     front.append(nI);
 280
 281                     do
 282                     {
 283                         DynList<label> newFront;
 284                         forAll(front, pI)
 285                         {
 286                             const label fLabel = front[pI];
 287                             if( nodeMaterial[fLabel] )
 288                                 continue;
 289                             nodeMaterial[fLabel] = mat;
 290                             edgeMats[mat-1].appendIfNotIn(f.faceEdge(fLabel));
 291                             edgeMats[mat-1].appendIfNotIn
 292                             (
 293                                 f.faceEdge(f.rcIndex(fLabel))
 294                             );
 295
 296                             if( removeFaceVertex[f.rcIndex(fLabel)] )
 297                                 newFront.append(f.rcIndex(fLabel));
 298                             if( removeFaceVertex[f.fcIndex(fLabel)] )
 299                                 newFront.append(f.fcIndex(fLabel));
 300                         }
 301                         front = newFront;
 302                     }
 303                     while( front.size() != 0 );
 304
 305                     ++mat;
 306                 }
 307
 308             forAll(edgeMats, matI)
 309             {
 310                 edgeMats[matI].shrink();
 311                 const face rf = help::removeEdgesFromFace(f, edgeMats[matI]);
 312                 const face newBf = help::createFaceFromRemovedPart(f, rf);
 313
 314                 # ifdef DEBUGMorph
 315                 Info << "Adding face " << newBf << " as boundary faces" << endl;
 316                 Info << "Owner and neighbour are " << owner[faceI] << " and "
 317                     << neighbour[faceI] << endl;
 318                 # endif
 319
 320                 newBoundaryFaces_.appendList(newBf);
 321                 newBoundaryOwners_.append(owner[faceI]);
 322                 newBoundaryPatches_.append(0);
 323
 324                 newBoundaryFaces_.appendList(newBf.reverseFace());
 325                 newBoundaryOwners_.append(neighbour[faceI]);
 326                 newBoundaryPatches_.append(0);
 327             }
 328
 329             const_cast<face&>(faces[faceI]) = newF;
 330         }
 331     }
 332
 333     //- treat processor boundaries
 334     if( Pstream::parRun() )
 335     {
 336         const PtrList<processorBoundaryPatch>& procBoundaries =
 337             mesh_.procBoundaries();
 338
 339         forAll(procBoundaries, patchI)
 340         {
 341             const label start = procBoundaries[patchI].patchStart();
 342             const label end = start + procBoundaries[patchI].patchSize();
 343             const bool isOwner = procBoundaries[patchI].owner();
 344
 345             for(label faceI=start;faceI<end;++faceI)
 346             {
 347                 face copy;
 348                 if( isOwner )
 349                 {
 350                     copy = faces[faceI];
 351                 }
 352                 else
 353                 {
 354                     copy = faces[faceI].reverseFace();
 355                 }
 356
 357                 const face& f = copy;
 358
 359                 DynList<bool> removeFaceVertex(f.size(), false);
 360
 361                 face newF(f.size());
 362                 label i(0);
 363
 364                 forAll(f, pI)
 365                     if(
 366                         boundaryVertex_[f.prevLabel(pI)] &&
 367                         boundaryVertex_[f[pI]] &&
 368                         boundaryVertex_[f.nextLabel(pI)]
 369                     )
 370                     {
 371                         removeFaceVertex[pI] = true;
 372
 373                         # ifdef DEBUGMorph
 374                         Info << "Removing vertex " << f[pI] << " from face "
 375                             << f << endl;
 376                         # endif
 377                     }
 378                     else
 379                     {
 380                         newF[i++] = f[pI];
 381                     }
 382
 383                 if( i < f.size() )
 384                 {
 385                     changed = true;
 386
 387                     //- store shrinked face
 388                     newF.setSize(i);
 389
 390                     # ifdef DEBUGMorph
 391                     Info << "Removing edges " << removeEdge << endl;
 392                     Info << "Face label " << faceI << " owner " << owner[faceI]
 393                         << " neighbour " << neighbour[faceI] << endl;
 394                     Info << "Internal face " << f << " has been shrinked into "
 395                         << newF << endl;
 396                     # endif
 397
 398                     //- create new boundary faces from the removed part
 399                     label mat(1);
 400                     DynList<direction> nodeMaterial(f.size(), direction(0));
 401                     DynList< DynList<edge> > edgeMats;
 402                     forAll(nodeMaterial, nI)
 403                         if( !nodeMaterial[nI] && removeFaceVertex[nI] )
 404                         {
 405                             edgeMats.append(DynList<edge>());
 406                             DynList<label> front;
 407                             front.append(nI);
 408
 409                             do
 410                             {
 411                                 DynList<label> newFront;
 412                                 forAll(front, pI)
 413                                 {
 414                                     const label fLabel = front[pI];
 415                                     if( nodeMaterial[fLabel] ) continue;
 416                                     nodeMaterial[fLabel] = mat;
 417                                     edgeMats[mat-1].appendIfNotIn
 418                                     (
 419                                         f.faceEdge(fLabel)
 420                                     );
 421                                     edgeMats[mat-1].appendIfNotIn
 422                                     (
 423                                         f.faceEdge(f.rcIndex(fLabel))
 424                                     );
 425
 426                                     if( removeFaceVertex[f.rcIndex(fLabel)] )
 427                                         newFront.append(f.rcIndex(fLabel));
 428                                     if( removeFaceVertex[f.fcIndex(fLabel)] )
 429                                         newFront.append(f.fcIndex(fLabel));
 430                                 }
 431                                 front = newFront;
 432                             }
 433                             while( front.size() != 0 );
 434
 435                             ++mat;
 436                         }
 437
 438                     forAll(edgeMats, matI)
 439                     {
 440                         edgeMats[matI].shrink();
 441                         const face rf =
 442                             help::removeEdgesFromFace(f, edgeMats[matI]);
 443                         const face newBf =
 444                             help::createFaceFromRemovedPart(f, rf);
 445
 446                         # ifdef DEBUGMorph
 447                         Info << "Adding face " << newBf
 448                         << " as boundary faces" << endl;
 449                         Info << "Owner and neighbour are "
 450                             << owner[faceI] << endl;
 451                         # endif
 452
 453                         if( isOwner )
 454                         {
 455                             newBoundaryFaces_.appendList(newBf);
 456                         }
 457                         else
 458                         {
 459                             newBoundaryFaces_.appendList(newBf.reverseFace());
 460                         }
 461                         newBoundaryOwners_.append(owner[faceI]);
 462                         newBoundaryPatches_.append(0);
 463                     }
 464
 465                     if( isOwner || (newF.size() < 3) )
 466                     {
 467                         const_cast<face&>(faces[faceI]) = newF;
 468                     }
 469                     else
 470                     {
 471                         const_cast<face&>(faces[faceI]) = newF.reverseFace();
 472                     }
 473                 }
 474             }
 475         }
 476     }
 477
 478     polyMeshGenModifier meshModifier(mesh_);
 479
 480     if( Pstream::parRun() )
 481     {
 482         reduce(changed, maxOp<bool>());
 483     }
 484
 485     if( changed )
 486     {
 487         //- replace boundary of the mesh
 488         replaceMeshBoundary();
 489
 490         # ifdef DEBUGMorph
 491         Serr << "Zip check 2" << endl;
 492         const cellListPMG& cells = mesh_.cells();
 493         forAll(cells, cI)
 494         {
 495             const cell& c = cells[cI];
 496
 497             DynList<edge> edges;
 498             DynList<direction> nAppearances;
 499
 500             forAll(c, fI)
 501             {
 502                 const face& f = faces[c[fI]];
 503                 forAll(f, eI)
 504                 {
 505                     const label pos = edges.containsAtPosition(f.faceEdge(eI));
 506
 507                     if( pos == -1 )
 508                     {
 509                         edges.append(f.faceEdge(eI));
 510                         nAppearances.append(1);
 511                     }
 512                     else
 513                     {
 514                         ++nAppearances[pos];
 515                     }
 516                 }
 517             }
 518
 519             forAll(nAppearances, eI)
 520                 if( nAppearances[eI] != 2 )
 521                 {
 522                     Warning << "Cell " << cI << " is not closed" << endl;
 523                     Serr << "Cell faces " << c << endl;
 524                     forAll(c, fI)
 525                         Serr << "Face " << c[fI] << " is "
 526                             << faces[c[fI]] << endl;
 527                 }
 528         }
 529         # endif
 530
 531         //- remove faces which do not exist any more
 532         boolList removeFace(faces.size(), false);
 533         bool removeFaces(false);
 534         for(register label faceI=0;faceI<nIntFaces;++faceI)
 535             if( faces[faceI].size() < 3 )
 536             {
 537                 removeFace[faceI] = true;
 538                 removeFaces = true;
 539             }
 540
 541         if( Pstream::parRun() )
 542         {
 543             const PtrList<processorBoundaryPatch>& procBoundaries =
 544                 mesh_.procBoundaries();
 545
 546             forAll(procBoundaries, patchI)
 547             {
 548                 const label start = procBoundaries[patchI].patchStart();
 549                 const label end = start + procBoundaries[patchI].patchSize();
 550                 for(label faceI=start;faceI<end;++faceI)
 551                     if( faces[faceI].size() < 3 )
 552                     {
 553                         removeFace[faceI] = true;
 554                         removeFaces = true;
 555                     }
 556             }
 557
 558             reduce(removeFaces, maxOp<bool>());
 559         }
 560
 561         if( removeFaces )
 562             meshModifier.removeFaces(removeFace);
 563     }
 564
 565     # ifdef DEBUGMorph
 566     Serr << "Zip check 3" << endl;
 567     labelHashSet zipCells;
 568     mesh_.addressingData().checkCellsZipUp(true, &zipCells);
 569     if( zipCells.size() )
 570     {
 571         Serr << Pstream::myProcNo() << "Open cells appeared!" << endl;
 572         //mesh_.write();
 573         Serr << "Cells " << zipCells << " are not zipped!!" << endl;
 574     }
 575     mesh_.clearAddressingData();
 576     # endif
 577
 578     Info << "Finished morphing internal faces" << endl;
 579
 580     return changed;
 581 }
 582
 583 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 584
 585 } // End namespace Foam
 586
 587 // ************************************************************************* //