src/mesh/cfMesh/meshLibrary/utilities/tetrahedra/tetCreatorOctree/tetCreatorOctreePointsAndAddressing.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 "meshOctree.H"
  30 #include "triSurface.H"
  31
  32 //#define DEBUGTets
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  40
  41 void tetCreatorOctree::selectElements()
  42 {
  43     const List<direction>& boxType = octreeCheck_.boxType();
  44     const meshOctree& octree = octreeCheck_.octree();
  45     const boundBox& rootBox = octree.rootBox();
  46
  47     //- store nodeLabels first
  48     const VRWGraph& nodeLabels = octreeCheck_.nodeLabels();
  49     tetPoints_.setSize(octreeCheck_.numberOfNodes());
  50
  51     forAll(nodeLabels, leafI)
  52         if( nodeLabels.sizeOfRow(leafI) != 0 )
  53         {
  54             const meshOctreeCubeBasic& oc = octree.returnLeaf(leafI);
  55             FixedList<point, 8> lv;
  56             oc.vertices(rootBox, lv);
  57
  58             forAll(lv, vI)
  59                 tetPoints_[nodeLabels(leafI, vI)] = lv[vI];
  60         }
  61
  62     //- create cubeLabel list
  63     if( !cubeLabelPtr_ )
  64         cubeLabelPtr_ = new labelList();
  65     labelList& cubeLabel = *cubeLabelPtr_;
  66     cubeLabel.setSize(octree.numberOfLeaves());
  67     cubeLabel = -1;
  68
  69     forAll(boxType, leafI)
  70         if( boxType[leafI] & meshOctreeAddressing::MESHCELL )
  71         {
  72             const meshOctreeCubeBasic& oc = octree.returnLeaf(leafI);
  73             cubeLabel[leafI] = tetPoints_.size();
  74             tetPoints_.append(oc.centre(rootBox));
  75         }
  76 }
  77
  78 void tetCreatorOctree::createPointsAndAddressing()
  79 {
  80     selectElements();
  81
  82     const meshOctree& octree = octreeCheck_.octree();
  83     const boundBox& rootBox = octree.rootBox();
  84     const VRWGraph& nodeLabels = octreeCheck_.nodeLabels();
  85
  86     if( !subNodeLabelsPtr_ )
  87         subNodeLabelsPtr_ = new VRWGraph(octree.numberOfLeaves());
  88     VRWGraph& subNodeLabels = *subNodeLabelsPtr_;
  89
  90     //- store nodeLabels
  91     direction maxLevel(0);
  92     forAll(nodeLabels, leafI)
  93         if( octree.returnLeaf(leafI).level() > maxLevel )
  94             maxLevel = octree.returnLeaf(leafI).level();
  95
  96     sortedLeaves_.setSize(maxLevel+1);
  97     forAll(sortedLeaves_, levelI)
  98         sortedLeaves_[levelI].clear();
  99
 100     forAll(nodeLabels, leafI)
 101     {
 102         const meshOctreeCubeBasic& oc = octree.returnLeaf(leafI);
 103         sortedLeaves_[oc.level()].append(leafI);
 104     }
 105
 106     //- create subNodeLabels
 107     const FRWGraph<label, 8>& pointLeaves = octreeCheck_.nodeLeaves();
 108
 109     forAll(pointLeaves, pointI)
 110     {
 111         bool validLeaf[8];
 112         direction levelI(0);
 113         for(label i=0;i<8;++i)
 114         {
 115             const label pointLeafI = pointLeaves(pointI, i);
 116
 117             if( pointLeafI == -1 )
 118             {
 119                 validLeaf[i] = false;
 120             }
 121             else
 122             {
 123                 validLeaf[i] = true;
 124                 for(label j=i+1;j<8;++j)
 125                     if( pointLeafI == pointLeaves(pointI, j) )
 126                     {
 127                         validLeaf[i] = false;
 128                         validLeaf[j] = false;
 129                     }
 130
 131                 const direction level =
 132                     octree.returnLeaf(pointLeafI).level();
 133
 134                 if( level > levelI )
 135                     levelI = level;
 136             }
 137         }
 138
 139         for(label plI=0;plI<8;++plI)
 140             if( validLeaf[plI] )
 141             {
 142                 const label pointLeafI = pointLeaves(pointI, plI);
 143
 144                 const meshOctreeCubeBasic& lc =
 145                     octree.returnLeaf(pointLeafI);
 146
 147                 if( lc.level() < levelI )
 148                 {
 149                     if( subNodeLabels.sizeOfRow(pointLeafI) != 8 )
 150                     {
 151                         subNodeLabels.setRowSize(pointLeafI, 8);
 152                         forAllRow(subNodeLabels, pointLeafI, k)
 153                             subNodeLabels(pointLeafI, k) = -1;
 154                     }
 155
 156                     subNodeLabels(pointLeafI, (7-plI)) = tetPoints_.size();
 157                     FixedList<point, 8> lv;
 158                     lc.vertices(rootBox, lv);
 159
 160                     tetPoints_.append
 161                     (
 162                         0.5 * (lv[7-plI] + lc.centre(rootBox))
 163                     );
 164                 }
 165             }
 166     }
 167
 168     createFaceCentreLabels();
 169
 170     # ifdef DEBUGTets
 171     forAll(nodeLabels, leafI)
 172     {
 173         forAllRow(nodeLabels, leafI, nlI)
 174             if( leafI != pointLeaves(nodeLabels(leafI, nlI), (7-nlI)) )
 175                 FatalError << "Shit" << abort(FatalError);
 176     }
 177     # endif
 178 }
 179
 180 void tetCreatorOctree::createFaceCentreLabels()
 181 {
 182     Info << "Creating face centre labels " << endl;
 183     const labelList& cubeLabel = *cubeLabelPtr_;
 184     const VRWGraph& nodeLabels = octreeCheck_.nodeLabels();
 185     const FRWGraph<label, 8>& pointLeaves = octreeCheck_.nodeLeaves();
 186     const meshOctree& octree = octreeCheck_.octree();
 187
 188
 189     # ifdef DEBUGTets
 190     Info << "Node labels " << nodeLabels << endl;
 191     # endif
 192
 193     List<direction> nodeLevel(pointLeaves.size(), direction(0));
 194     forAll(nodeLabels, leafI)
 195     {
 196         const direction level = octree.returnLeaf(leafI).level();
 197
 198         forAllRow(nodeLabels, leafI, nlI)
 199         {
 200             const label nLabel = nodeLabels(leafI, nlI);
 201             # ifdef DEBUGTets
 202             Info << "Node label[" << leafI << "][" << nlI << "] "
 203                 << nLabel << endl;
 204             # endif
 205
 206             if( nodeLevel[nLabel] < level )
 207                 nodeLevel[nLabel] = level;
 208         }
 209     }
 210
 211     if( !faceCentreLabelPtr_ )
 212         faceCentreLabelPtr_ = new VRWGraph(cubeLabel.size());
 213     VRWGraph& faceCentreLabel = *faceCentreLabelPtr_;
 214
 215     forAll(cubeLabel, cubeI)
 216         if( cubeLabel[cubeI] != -1 )
 217         {
 218             const direction level = octree.returnLeaf(cubeI).level();
 219
 220             for(label i=0;i<6;++i)
 221             {
 222                 if(
 223                     (faceCentreLabel.sizeOfRow(cubeI) != 0) &&
 224                     (faceCentreLabel(cubeI, i) != -1)
 225                 )
 226                     continue;
 227
 228                 FixedList<label, 4> faceNodes;
 229                 forAll(faceNodes, fnI)
 230                     faceNodes[fnI] =
 231                         meshOctreeCubeCoordinates::faceNodes_[i][fnI];
 232
 233                 label highLevelNode(-1);
 234                 for(label j=0;j<4;++j)
 235                     if( nodeLevel[nodeLabels(cubeI, faceNodes[j])] > level )
 236                     {
 237                         highLevelNode = j;
 238                         break;
 239                     }
 240
 241                 if( highLevelNode == -1 )
 242                     continue;
 243
 244                 DynList<label> neighbours;
 245                 octree.findNeighboursInDirection(cubeI, i, neighbours);
 246
 247                 if( (neighbours.size() != 1) || (neighbours[0] == -1) )
 248                     continue;
 249
 250                 if( faceCentreLabel.sizeOfRow(cubeI) == 0 )
 251                 {
 252                     faceCentreLabel.setRowSize(cubeI, 6);
 253                     forAllRow(faceCentreLabel, cubeI, colI)
 254                         faceCentreLabel(cubeI, colI) = -1;
 255                 }
 256                 const label cNei = neighbours[0];
 257                 if( faceCentreLabel.sizeOfRow(cNei) == 0 )
 258                 {
 259                     faceCentreLabel.setRowSize(cNei, 6);
 260                     forAllRow(faceCentreLabel, cNei, colI)
 261                         faceCentreLabel(cNei, colI) = -1;
 262                 }
 263
 264                 faceCentreLabel(cubeI, i) = tetPoints_.size();
 265                 switch( i )
 266                 {
 267                     case 0:
 268                     {
 269                         faceCentreLabel(cNei, 1) = tetPoints_.size();
 270                     } break;
 271                     case 1:
 272                     {
 273                         faceCentreLabel(cNei, 0) = tetPoints_.size();
 274                     } break;
 275                     case 2:
 276                     {
 277                         faceCentreLabel(cNei, 3) = tetPoints_.size();
 278                     } break;
 279                     case 3:
 280                     {
 281                         faceCentreLabel(cNei, 2) = tetPoints_.size();
 282                     } break;
 283                     case 4:
 284                     {
 285                         faceCentreLabel(cNei, 5) = tetPoints_.size();
 286                     } break;
 287                     case 5:
 288                     {
 289                         faceCentreLabel(cNei, 4) = tetPoints_.size();
 290                     } break;
 291                 };
 292
 293                 point p(vector::zero);
 294                 for(label j=0;j<4;++j)
 295                     p += tetPoints_[nodeLabels(cubeI, faceNodes[j])];
 296                 p /= 4;
 297                 tetPoints_.append(p);
 298             }
 299         }
 300 }
 301
 302 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 303
 304 } // End namespace Foam
 305
 306 // ************************************************************************* //