src/mesh/cfMesh/utilities/tetrahedra/tetCreatorOctree/tetCreatorOctreeTetsAroundSplitEdges.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 "tetCreatorOctree.H"
  29 #include "demandDrivenData.H"
  30 #include "meshOctree.H"
  31
  32 //#define DEBUGTets
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 void tetCreatorOctree::createTetsAroundSplitEdges()
  42 {
  43     Info << "Creating tets around split edges " << endl;
  44
  45     const labelList& cubeLabel = *cubeLabelPtr_;
  46     const meshOctree& octree = octreeCheck_.octree();
  47     const VRWGraph& nodeLabels = octreeCheck_.nodeLabels();
  48     const FRWGraph<label, 8>& pLeaves = octreeCheck_.nodeLeaves();
  49     const VRWGraph& subNodeLabels = *subNodeLabelsPtr_;
  50     const VRWGraph& faceCentreLabel = *faceCentreLabelPtr_;
  51
  52     # ifdef DEBUGTets
  53     Info << "Number of octree nodes " << octreeCheck_.numberOfNodes() << endl;
  54     # endif
  55
  56     //- find maximum refinement level of octree leaves attached to each vertex
  57     List<direction> nodeLevel(octreeCheck_.numberOfNodes());
  58
  59     forAll(pLeaves, nodeI)
  60     {
  61         direction level(0);
  62
  63         for(label plI=0;plI<8;++plI)
  64         {
  65             const label leafI = pLeaves(nodeI, plI);
  66
  67             if( leafI < 0 )
  68                 continue;
  69
  70             level = Foam::max(level, octree.returnLeaf(leafI).level());
  71         }
  72
  73         nodeLevel[nodeI] = level;
  74     }
  75
  76     //- start creating tets around split edges
  77     label helpNodes[2][8];
  78     label faceCentres[4];
  79
  80     forAllReverse(sortedLeaves_, levelI)
  81     {
  82         const labelLongList& curLevelLeaves = sortedLeaves_[levelI];
  83
  84         const direction level = direction(levelI);
  85
  86         forAll(curLevelLeaves, leafI)
  87         {
  88             const label curLabel = curLevelLeaves[leafI];
  89
  90             if( cubeLabel[curLabel] == -1 )
  91                 continue;
  92
  93             //- start checking edges
  94             for(label eI=0;eI<12;++eI)
  95             {
  96                 const label start =
  97                     nodeLabels
  98                     (
  99                         curLabel,
 100                         meshOctreeCubeCoordinates::edgeNodes_[eI][0]
 101                     );
 102                 const label end =
 103                     nodeLabels
 104                     (
 105                         curLabel,
 106                         meshOctreeCubeCoordinates::edgeNodes_[eI][1]
 107                     );
 108
 109                 if( (nodeLevel[start] == level) && (nodeLevel[end] == level) )
 110                     continue;
 111
 112                 //- the edge has at least one vertex at different ref level
 113                 bool create(true);
 114
 115                 FixedList<label, 4> edgeCubes;
 116                 const label fI = 2*(eI/4)+1;
 117
 118                 const label* fNodes =
 119                     meshOctreeCubeCoordinates::faceNodes_[fI];
 120
 121                 //- store octree leaves at this edge
 122                 //- they are all adjacent to the start point
 123                 for(label i=0;i<4;++i)
 124                     edgeCubes[i] = pLeaves(start, fNodes[i]);
 125
 126                 # ifdef DEBUGTets
 127                 Info << "Cube " << curLabel << " has nodes ";
 128                 forAllRow(nodeLabels, curLabel, i)
 129                     Info << nodeLabels(curLabel, i) << " ";
 130                 Info << endl;
 131                 Info << "Creating tets around edge " << eI << endl;
 132                 Info << "Edge nodes are " << start << " and " << end << endl;
 133                 Info << "Edge cubes " << edgeCubes << endl;
 134                 # endif
 135
 136                 forAll(edgeCubes, i)
 137                 {
 138                     const label cLabel = edgeCubes[i];
 139
 140                     if(
 141                         (cLabel == -1) ||
 142                         (cLabel < curLabel) ||
 143                         (octree.returnLeaf(cLabel).level() != level)
 144                     )
 145                     {
 146                         create = false;
 147                         break;
 148                     }
 149
 150                     # ifdef DEBUGTets
 151                     Info << "Edge cube " << i << " is " << cLabel << endl;
 152                     # endif
 153
 154                     for(label j=0;j<8;++j)
 155                     {
 156                         if( nodeLabels(cLabel,j) == start )
 157                         {
 158                             if( subNodeLabels.sizeOfRow(cLabel) != 0 )
 159                             {
 160                                 helpNodes[0][i] = subNodeLabels(cLabel,j);
 161                             }
 162                             else
 163                             {
 164                                 helpNodes[0][i] = -1;
 165                             }
 166
 167                             helpNodes[0][i+4] = cubeLabel[cLabel];
 168                         }
 169
 170                         if( nodeLabels(cLabel, j) == end )
 171                         {
 172                             if( subNodeLabels.sizeOfRow(cLabel) != 0 )
 173                             {
 174                                 helpNodes[1][i+4] = subNodeLabels(cLabel,j);
 175                             }
 176                             else
 177                             {
 178                                 helpNodes[1][i+4] = -1;
 179                             }
 180
 181                             helpNodes[1][i] = cubeLabel[cLabel];
 182                         }
 183                     }
 184
 185                     if( faceCentreLabel.sizeOfRow(cLabel) != 0 )
 186                     {
 187                         const label helpFace = eI/4;
 188
 189                         faceCentres[i] =
 190                             faceCentreLabel
 191                             (
 192                                 cLabel,
 193                                 faceCentreHelper_[helpFace][i]
 194                             );
 195                     }
 196                     else
 197                     {
 198                         faceCentres[i] = -1;
 199                     }
 200                 }
 201
 202                 if( !create )
 203                     continue;
 204
 205                 # ifdef DEBUGTets
 206                 for(label n=0;n<4;++n)
 207                 {
 208                     Info << "Face centre " << n << "  "
 209                         << faceCentres[n] << endl;
 210                     Info << "Hex 0 " << helpNodes[0][n] << " and "
 211                         << helpNodes[0][n+4] << endl;
 212                     Info << "Hex 1 " << helpNodes[1][n] << " and "
 213                         << helpNodes[1][n+4] << endl;
 214                 }
 215                 # endif
 216
 217                 if( nodeLevel[start] > level )
 218                 {
 219                     for(label k=0;k<4;++k)
 220                     {
 221                         //- add 4 tets
 222                         checkAndAppendTet
 223                         (
 224                             partTet
 225                             (
 226                                 start,
 227                                 helpNodes[0][k],
 228                                 helpNodes[0][(k+1)%4],
 229                                 tetPoints_.size()
 230                             )
 231                         );
 232
 233                         //- 2. tet
 234                         checkAndAppendTet
 235                         (
 236                             partTet
 237                             (
 238                                 tetPoints_.size(),
 239                                 helpNodes[0][k],
 240                                 helpNodes[0][(k+1)%4],
 241                                 faceCentres[k]
 242                             )
 243                         );
 244
 245                         //- 3. tet
 246                         checkAndAppendTet
 247                         (
 248                             partTet
 249                             (
 250                                 faceCentres[k],
 251                                 helpNodes[0][k],
 252                                 helpNodes[0][k+4],
 253                                 tetPoints_.size()
 254                             )
 255                         );
 256
 257                         //- 4. tet
 258                         checkAndAppendTet
 259                         (
 260                             partTet
 261                             (
 262                                 faceCentres[(k+3)%4],
 263                                 helpNodes[0][k+4],
 264                                 helpNodes[0][k],
 265                                 tetPoints_.size()
 266                             )
 267                         );
 268                     }
 269                 }
 270                 else
 271                 {
 272                     for(label k=0;k<4;++k)
 273                     {
 274                         checkAndAppendTet
 275                         (
 276                             partTet
 277                             (
 278                                 faceCentres[(k+3)%4],
 279                                 helpNodes[0][k+4],
 280                                 start,
 281                                 tetPoints_.size()
 282                             )
 283                         );
 284
 285                         checkAndAppendTet
 286                         (
 287                             partTet
 288                             (
 289                                 faceCentres[k],
 290                                 start,
 291                                 helpNodes[0][k+4],
 292                                 tetPoints_.size()
 293                             )
 294                         );
 295                     }
 296                 }
 297
 298                 if( nodeLevel[end] > level )
 299                 {
 300                     for(label k=0;k<4;++k)
 301                     {
 302                         //- add 4 tets
 303                         checkAndAppendTet
 304                         (
 305                             partTet
 306                             (
 307                                 tetPoints_.size(),
 308                                 helpNodes[1][k+4],
 309                                 helpNodes[1][((k+1)%4)+4],
 310                                 end
 311                             )
 312                         );
 313
 314                         // 2. tet
 315                         checkAndAppendTet
 316                         (
 317                             partTet
 318                             (
 319                                 faceCentres[k],
 320                                 helpNodes[1][k+4],
 321                                 helpNodes[1][((k+1)%4)+4],
 322                                 tetPoints_.size()
 323                             )
 324                         );
 325
 326                         //- 3. tet
 327                         checkAndAppendTet
 328                         (
 329                             partTet
 330                             (
 331                                 tetPoints_.size(),
 332                                 helpNodes[1][k],
 333                                 faceCentres[k],
 334                                 helpNodes[1][k+4]
 335                             )
 336                         );
 337
 338                         //- 4. tet
 339                         checkAndAppendTet
 340                         (
 341                             partTet
 342                             (
 343                                 faceCentres[(k+3)%4],
 344                                 helpNodes[1][k],
 345                                 tetPoints_.size(),
 346                                 helpNodes[1][k+4]
 347                             )
 348                         );
 349                     }
 350                 }
 351                 else
 352                 {
 353                     for(label k=0;k<4;++k)
 354                     {
 355                         checkAndAppendTet
 356                         (
 357                             partTet
 358                             (
 359                                 faceCentres[(k+3)%4],
 360                                 helpNodes[1][k],
 361                                 tetPoints_.size(),
 362                                 end
 363                             )
 364                         );
 365
 366                         checkAndAppendTet
 367                         (
 368                             partTet
 369                             (
 370                                 helpNodes[1][k],
 371                                 faceCentres[k],
 372                                 tetPoints_.size(),
 373                                 end
 374                             )
 375                         );
 376                     }
 377                 }
 378
 379                 //- add the edge centre
 380                 tetPoints_.append(0.5 * (tetPoints_[start] + tetPoints_[end]));
 381             }
 382         }
 383     }
 384
 385     # ifdef DEBUGTets
 386     Info << "Created tets " << tets_ << endl;
 387     # endif
 388 }
 389
 390 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 391
 392 } // End namespace Foam
 393
 394 // ************************************************************************* //