src/mesh/cfMesh/utilities/meshes/partTetMesh/partTetMeshAddressing.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 "demandDrivenData.H"
  29 #include "polyMeshGenModifier.H"
  30 #include "partTetMesh.H"
  31 #include "tetMatcher.H"
  32 #include "polyMeshGenAddressing.H"
  33 #include "helperFunctionsPar.H"
  34
  35 //#define DEBUGSmooth
  36
  37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  38
  39 namespace Foam
  40 {
  41
  42 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  43
  44 void partTetMesh::createPointsAndTets
  45 (
  46     const List<direction>& useCell,
  47     const boolList& lockedPoints
  48 )
  49 {
  50     const pointFieldPMG& points = origMesh_.points();
  51     const faceListPMG& faces = origMesh_.faces();
  52     const cellListPMG& cells = origMesh_.cells();
  53     const labelList& owner = origMesh_.owner();
  54     const labelList& neighbour = origMesh_.neighbour();
  55     const PtrList<boundaryPatch>& boundaries = origMesh_.boundaries();
  56     const PtrList<processorBoundaryPatch>& procBoundaries =
  57         origMesh_.procBoundaries();
  58     const label nInternalFaces = origMesh_.nInternalFaces();
  59
  60     //- check how many neighbours of a face are marked for smoothing
  61     labelList usedFace(faces.size(), 0);
  62
  63     //- mark faces
  64     forAll(faces, faceI)
  65     {
  66         if( useCell[owner[faceI]] )
  67             ++usedFace[faceI];
  68
  69         if( neighbour[faceI] < 0 )
  70             continue;
  71
  72         if( useCell[neighbour[faceI]] )
  73             ++usedFace[faceI];
  74     }
  75
  76     //- send data at processor boundaries
  77     forAll(procBoundaries, patchI)
  78     {
  79         const label start = procBoundaries[patchI].patchStart();
  80         const label size = procBoundaries[patchI].patchSize();
  81
  82         labelLongList dataToSend;
  83         for(label faceI=0;faceI<size;++faceI)
  84         {
  85             if( usedFace[start+faceI] )
  86                 dataToSend.append(faceI);
  87         }
  88
  89         OPstream toOtherProc
  90         (
  91             Pstream::blocking,
  92             procBoundaries[patchI].neiProcNo(),
  93             dataToSend.byteSize()
  94         );
  95
  96         toOtherProc << dataToSend;
  97     }
  98
  99     //- receive data at proc boundaries
 100     forAll(procBoundaries, patchI)
 101     {
 102         labelLongList receivedData;
 103
 104         IPstream fromOtherProc
 105         (
 106             Pstream::blocking,
 107             procBoundaries[patchI].neiProcNo()
 108         );
 109
 110         fromOtherProc >> receivedData;
 111
 112         const label start = procBoundaries[patchI].patchStart();
 113         forAll(receivedData, faceI)
 114             ++usedFace[start+receivedData[faceI]];
 115     }
 116
 117     const vectorField& faceCentres = origMesh_.addressingData().faceCentres();
 118     const vectorField& cellCentres = origMesh_.addressingData().cellCentres();
 119
 120     labelLongList nodeLabelForPoint(points.size(), -1);
 121     labelLongList nodeLabelForFace(faces.size(), -1);
 122     labelLongList nodeLabelForCell(cells.size(), -1);
 123
 124     points_.clear();
 125     smoothVertex_.clear();
 126
 127     //- create BOUNDARY points
 128     forAll(boundaries, patchI)
 129     {
 130         const boundaryPatch& patch = boundaries[patchI];
 131         const label start = patch.patchStart();
 132         const label end = start + patch.patchSize();
 133
 134         for(label faceI=start;faceI<end;++faceI)
 135         {
 136             if( !usedFace[faceI] )
 137                 continue;
 138
 139             const face& f = faces[faceI];
 140
 141             if( f.size() > 3 )
 142             {
 143                 //- create face centre
 144                 nodeLabelForFace[faceI] = points_.size();
 145                 points_.append(faceCentres[faceI]);
 146                 smoothVertex_.append(FACECENTRE);
 147             }
 148
 149             //- add face corners
 150             forAll(f, pI)
 151             {
 152                 const label pointI = f[pI];
 153                 if( nodeLabelForPoint[pointI] == -1 )
 154                 {
 155                     nodeLabelForPoint[pointI] = points_.size();
 156                     points_.append(points[pointI]);
 157
 158                     smoothVertex_.append(BOUNDARY);
 159                 }
 160             }
 161         }
 162     }
 163
 164     //- create points at processor boundaries
 165     forAll(procBoundaries, patchI)
 166     {
 167         const processorBoundaryPatch& patch = procBoundaries[patchI];
 168         const label start = patch.patchStart();
 169         const label end = start + patch.patchSize();
 170
 171         for(label faceI=start;faceI<end;++faceI)
 172         {
 173             if( !usedFace[faceI] )
 174                 continue;
 175
 176             const face& f = faces[faceI];
 177
 178             if( f.size() > 3 )
 179             {
 180                 //- create face centre
 181                 nodeLabelForFace[faceI] = points_.size();
 182                 points_.append(faceCentres[faceI]);
 183                 smoothVertex_.append(FACECENTRE);
 184             }
 185
 186             //- add face corners
 187             const direction vType = usedFace[faceI]==2?SMOOTH:NONE;
 188             forAll(f, pI)
 189             {
 190                 const label pointI = f[pI];
 191                 if( nodeLabelForPoint[pointI] == -1 )
 192                 {
 193                     nodeLabelForPoint[pointI] = points_.size();
 194                     points_.append(points[pointI]);
 195
 196                     smoothVertex_.append(vType);
 197                 }
 198                 else if( vType == NONE )
 199                 {
 200                      smoothVertex_[nodeLabelForPoint[pointI]] = NONE;
 201                 }
 202             }
 203         }
 204     }
 205
 206     //- create points for internal faces
 207     for(label faceI=0;faceI<nInternalFaces;++faceI)
 208     {
 209         if( usedFace[faceI] )
 210         {
 211             const face& f = faces[faceI];
 212
 213             if( f.size() > 3 )
 214             {
 215                 //- create face centre
 216                 nodeLabelForFace[faceI] = points_.size();
 217                 points_.append(faceCentres[faceI]);
 218                 smoothVertex_.append(FACECENTRE);
 219             }
 220
 221             //- add face corners
 222             const direction vType = usedFace[faceI]==2?SMOOTH:NONE;
 223             forAll(f, pI)
 224             {
 225                 const label pointI = f[pI];
 226                 if( nodeLabelForPoint[pointI] == -1 )
 227                 {
 228                     nodeLabelForPoint[pointI] = points_.size();
 229                     points_.append(points[pointI]);
 230
 231                     smoothVertex_.append(vType);
 232                 }
 233                 else if( vType == NONE )
 234                 {
 235                      smoothVertex_[nodeLabelForPoint[pointI]] = NONE;
 236                 }
 237             }
 238         }
 239     }
 240
 241     //- create tets
 242     tetMatcher tet;
 243     forAll(useCell, cI)
 244         if( useCell[cI] )
 245         {
 246             const cell& c = cells[cI];
 247
 248             if( tet.matchShape(false, faces, owner, cI, cells[cI]) )
 249             {
 250                 const labelList& tVrt = tet.vertLabels();
 251
 252                 //- add tet
 253                 tets_.append
 254                 (
 255                     partTet
 256                     (
 257                         nodeLabelForPoint[tVrt[0]],
 258                         nodeLabelForPoint[tVrt[1]],
 259                         nodeLabelForPoint[tVrt[2]],
 260                         nodeLabelForPoint[tVrt[3]]
 261                     )
 262                 );
 263
 264                 continue;
 265             }
 266
 267             nodeLabelForCell[cI] = points_.size();
 268             const label centreLabel = points_.size();
 269             points_.append(cellCentres[cI]);
 270             smoothVertex_.append(CELLCENTRE);
 271
 272             forAll(c, fI)
 273             {
 274                 const face& f = faces[c[fI]];
 275
 276                 if( owner[c[fI]] == cI )
 277                 {
 278                     if( f.size() == 3 )
 279                     {
 280                         partTet tet
 281                         (
 282                             nodeLabelForPoint[f[0]],
 283                             nodeLabelForPoint[f[2]],
 284                             nodeLabelForPoint[f[1]],
 285                             centreLabel
 286                         );
 287
 288                         # ifdef DEBUGSmooth
 289                         Info << "1.1 Tet " << tets_.size() << " is "
 290                             << tet << endl;
 291                         # endif
 292
 293                         tets_.append(tet);
 294                     }
 295                     else
 296                     {
 297                         forAll(f, pI)
 298                         {
 299                             partTet tet
 300                             (
 301                                 nodeLabelForPoint[f[pI]],
 302                                 nodeLabelForPoint[f.prevLabel(pI)],
 303                                 nodeLabelForFace[c[fI]],
 304                                 centreLabel
 305                             );
 306
 307                             # ifdef DEBUGSmooth
 308                             Info << "1.2 Tet " << tets_.size() << " is "
 309                                 << tet << endl;
 310                             # endif
 311
 312                             tets_.append(tet);
 313                         }
 314                     }
 315                 }
 316                 else
 317                 {
 318                     if( f.size() == 3 )
 319                     {
 320                         partTet tet
 321                         (
 322                             nodeLabelForPoint[f[0]],
 323                             nodeLabelForPoint[f[1]],
 324                             nodeLabelForPoint[f[2]],
 325                             centreLabel
 326                         );
 327
 328                         # ifdef DEBUGSmooth
 329                         Info << "2.1 Tet " << tets_.size() << " is "
 330                             << tet << endl;
 331                         # endif
 332
 333                         tets_.append(tet);
 334                     }
 335                     else
 336                     {
 337                         forAll(f, pI)
 338                         {
 339                             partTet tet
 340                             (
 341                                 nodeLabelForPoint[f[pI]],
 342                                 nodeLabelForPoint[f.nextLabel(pI)],
 343                                 nodeLabelForFace[c[fI]],
 344                                 centreLabel
 345                             );
 346
 347                             # ifdef DEBUGSmooth
 348                             Info << "2.2 Tet " << tets_.size() << " is "
 349                                 << tet << endl;
 350                             # endif
 351
 352                             tets_.append(tet);
 353                         }
 354                     }
 355                 }
 356             }
 357         }
 358
 359     //- create node labels in origMesh_
 360     nodeLabelInOrigMesh_.setSize(points_.size());
 361     nodeLabelInOrigMesh_ = -1;
 362     forAll(nodeLabelForPoint, pI)
 363         if( nodeLabelForPoint[pI] != -1 )
 364         {
 365             //- lock mesh vertices
 366             if( lockedPoints[pI] )
 367                 smoothVertex_[nodeLabelForPoint[pI]] |= LOCKED;
 368
 369             nodeLabelInOrigMesh_[nodeLabelForPoint[pI]] = pI;
 370         }
 371
 372
 373     //- create pointTets_
 374     pointTets_.reverseAddressing(points_.size(), tets_);
 375
 376     //- create addressing for parallel runs
 377     if( Pstream::parRun() )
 378     {
 379         createParallelAddressing
 380         (
 381             nodeLabelForPoint,
 382             nodeLabelForFace,
 383             nodeLabelForCell
 384         );
 385
 386         createBufferLayers();
 387     }
 388
 389     # ifdef DEBUGSmooth
 390     forAll(nodeLabelInOrigMesh_, pI)
 391         if(
 392             (nodeLabelInOrigMesh_[pI] != -1) &&
 393             (mag(points_[pI] - points[nodeLabelInOrigMesh_[pI]]) > SMALL)
 394         )
 395             FatalErrorIn
 396             (
 397                 "void partTetMesh::createPointsAndTets"
 398                 "(const boolList& useCell)"
 399             ) << "Node " << pI << " is dislocated" << abort(FatalError);
 400     # endif
 401 }
 402
 403 void partTetMesh::createSMOOTHPointsOrdering() const
 404 {
 405     internalPointsOrderPtr_ = new VRWGraph();
 406     VRWGraph& internalPointsOrder = *internalPointsOrderPtr_;
 407
 408     internalPointsOrder.setSize(0);
 409     labelLongList order(points_.size(), -1);
 410     boolList helper(points_.size());
 411
 412     bool found;
 413     do
 414     {
 415         found = false;
 416         helper = false;
 417         labelLongList selectedPoints;
 418
 419         forAll(points_, nodeI)
 420         {
 421             if( smoothVertex_[nodeI] & SMOOTH )
 422             {
 423                 if( helper[nodeI] )
 424                     continue;
 425                 if( order[nodeI] != -1 )
 426                     continue;
 427
 428                 //- find neighbouring FACECENTRE and CELLCENTRE points
 429                 DynList<label, 64> neiCentrePoints, neiSmoothPoints;
 430                 forAllRow(pointTets_, nodeI, ptI)
 431                 {
 432                     const partTet& tet = tets_[pointTets_(nodeI, ptI)];
 433
 434                     for(label i=0;i<4;++i)
 435                         if( smoothVertex_[tet[i]] & (FACECENTRE+CELLCENTRE) )
 436                         {
 437                             neiCentrePoints.appendIfNotIn(tet[i]);
 438                         }
 439                         else if( smoothVertex_[tet[i]] & SMOOTH )
 440                         {
 441                             neiSmoothPoints.appendIfNotIn(tet[i]);
 442                         }
 443                 }
 444
 445                 //- find neighbouring SMOOTH points
 446                 forAll(neiCentrePoints, ncI)
 447                 {
 448                     const label centreI = neiCentrePoints[ncI];
 449
 450                     forAllRow(pointTets_, centreI, ptI)
 451                     {
 452                         const partTet& tet = tets_[pointTets_(centreI, ptI)];
 453
 454                         for(label i=0;i<4;++i)
 455                             if( smoothVertex_[tet[i]] & SMOOTH )
 456                                 neiSmoothPoints.appendIfNotIn(tet[i]);
 457                     }
 458                 }
 459
 460                 //- select the point and mark neighbouring SMOOTH points
 461                 selectedPoints.append(nodeI);
 462                 order[nodeI] = internalPointsOrder.size();
 463
 464                 forAll(neiSmoothPoints, i)
 465                     helper[neiSmoothPoints[i]] = true;
 466             }
 467         }
 468
 469         if( selectedPoints.size() != 0 )
 470         {
 471             internalPointsOrder.appendList(selectedPoints);
 472             found = true;
 473         }
 474
 475     } while( found );
 476 }
 477
 478 void partTetMesh::createBOUNDARYPointsOrdering() const
 479 {
 480     boundaryPointsOrderPtr_ = new VRWGraph();
 481     VRWGraph& boundaryPointsOrder = *boundaryPointsOrderPtr_;
 482
 483     boundaryPointsOrder.setSize(0);
 484     labelLongList order(points_.size(), -1);
 485     boolList helper(points_.size());
 486
 487     bool found;
 488     do
 489     {
 490         found = false;
 491         helper = false;
 492
 493         labelLongList selectedPoints;
 494         forAll(points_, nodeI)
 495         {
 496             if( smoothVertex_[nodeI] & BOUNDARY )
 497             {
 498                 if( helper[nodeI] )
 499                     continue;
 500                 if( order[nodeI] != -1 )
 501                     continue;
 502
 503                 //- find neighbouring FACECENTRE and CELLCENTRE points
 504                 DynList<label, 64> neiCentrePoints, neiSmoothPoints;
 505                 forAllRow(pointTets_, nodeI, ptI)
 506                 {
 507                     const partTet& tet = tets_[pointTets_(nodeI, ptI)];
 508
 509                     for(label i=0;i<4;++i)
 510                         if( smoothVertex_[tet[i]] & (FACECENTRE+CELLCENTRE) )
 511                         {
 512                             neiCentrePoints.appendIfNotIn(tet[i]);
 513                         }
 514                         else if( smoothVertex_[tet[i]] & BOUNDARY )
 515                         {
 516                             neiSmoothPoints.appendIfNotIn(tet[i]);
 517                         }
 518                 }
 519
 520                 //- find neighbouring BOUNDARY points
 521                 forAll(neiCentrePoints, ncI)
 522                 {
 523                     const label centreI = neiCentrePoints[ncI];
 524
 525                     forAllRow(pointTets_, centreI, ptI)
 526                     {
 527                         const partTet& tet = tets_[pointTets_(centreI, ptI)];
 528
 529                         for(label i=0;i<4;++i)
 530                             if( smoothVertex_[tet[i]] & BOUNDARY )
 531                                 neiSmoothPoints.appendIfNotIn(tet[i]);
 532                     }
 533                 }
 534
 535                 //- select the point and mark neighbouring  BOUNDARY points
 536                 selectedPoints.append(nodeI);
 537                 order[nodeI] = boundaryPointsOrder.size();
 538
 539                 forAll(neiSmoothPoints, i)
 540                     helper[neiSmoothPoints[i]] = true;
 541             }
 542         }
 543
 544         if( selectedPoints.size() != 0 )
 545         {
 546             boundaryPointsOrder.appendList(selectedPoints);
 547             found = true;
 548         }
 549
 550     } while( found );
 551 }
 552
 553 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 554
 555 } // End namespace Foam
 556
 557 // ************************************************************************* //