src/meshLibrary/utilities/tetrahedra/tetCreatorOctree/tetCreatorOctreeTetsAroundEdges.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::createTetsAroundEdges()
  42 {
  43     Info << "Creating tets around edges" << endl;
  44
  45     const labelList& cubeLabel = *cubeLabelPtr_;
  46     const VRWGraph& nodeLabels = octreeCheck_.nodeLabels();
  47     const FRWGraph<label, 8>& pLeaves = octreeCheck_.nodeLeaves();
  48     const VRWGraph& subNodeLabels = *subNodeLabelsPtr_;
  49     const VRWGraph& faceCentreLabel = *faceCentreLabelPtr_;
  50     const meshOctree& octree = octreeCheck_.octree();
  51     const FixedList<FixedList<meshOctreeCubeCoordinates, 8>, 8>& vlPos =
  52         octree.positionsOfLeavesAtNodes();
  53
  54     //- find maximum refinement level of octree leaves attached to each vertex
  55     List<direction> nodeLevel(octreeCheck_.numberOfNodes());
  56
  57     forAll(pLeaves, nodeI)
  58     {
  59         direction level(0);
  60
  61         for(label plI=0;plI<8;++plI)
  62         {
  63             const label leafI = pLeaves(nodeI, plI);
  64
  65             if( leafI < 0 )
  66                 continue;
  67
  68             level = Foam::max(level, octree.returnLeaf(leafI).level());
  69         }
  70
  71         nodeLevel[nodeI] = level;
  72     }
  73
  74     //- start creating tets around edges which have both vertices at the same
  75     //- refinement level which is equal to the max refinement level of boxes
  76     //- incident to such edges
  77     forAllReverse(sortedLeaves_, levelI)
  78     {
  79         const labelLongList& curLevelLeaves = sortedLeaves_[levelI];
  80
  81         forAll(curLevelLeaves, leafI)
  82         {
  83             const label curLeaf = curLevelLeaves[leafI];
  84
  85             if( cubeLabel[curLeaf] == -1 )
  86                 continue;
  87
  88             const meshOctreeCubeCoordinates& cc =
  89                 octree.returnLeaf(curLeaf).coordinates();
  90
  91             # ifdef DEBUGTets
  92             Info << "Search cube " << curLeaf << " has coordinates "
  93                 << octree.returnLeaf(curLeaf).coordinates() << endl;
  94             Info << "Node labels for cube are " << nodeLabels[curLeaf] << endl;
  95             # endif
  96
  97             //- start checking edges
  98             for(label eI=0;eI<12;++eI)
  99             {
 100                 const label startNode =
 101                     meshOctreeCubeCoordinates::edgeNodes_[eI][0];
 102
 103                 const label start = nodeLabels(curLeaf, startNode);
 104                 const label end =
 105                     nodeLabels
 106                     (
 107                         curLeaf,
 108                         meshOctreeCubeCoordinates::edgeNodes_[eI][1]
 109                     );
 110
 111                 # ifdef DEBUGTets
 112                 Info << "Creating tets around edge " << eI << endl;
 113                 Info << "Edge nodes are " << start << " and " << end << endl;
 114                 Info << "Coordinates start " << tetPoints_[start]
 115                      << " end " << tetPoints_[end] << endl;
 116                 # endif
 117
 118                 bool create(true);
 119
 120                 if(
 121                     (nodeLevel[start] == direction(levelI)) &&
 122                     (nodeLevel[end] == direction(levelI))
 123                 )
 124                 {
 125                     //- edge has both vertices at the same refinement level
 126                     //- as the current leaf
 127                     FixedList<label, 4> edgeCubes;
 128                     const label fI = 2*(eI/4)+1;
 129
 130                     const label* fNodes =
 131                         meshOctreeCubeCoordinates::faceNodes_[fI];
 132
 133                     //- store octree leaves at this edge
 134                     //- they are all adjacent to the start point
 135                     for(label i=0;i<4;++i)
 136                         edgeCubes[i] = pLeaves(start, fNodes[i]);
 137
 138                     # ifdef DEBUGTets
 139                     forAll(edgeCubes, i)
 140                     {
 141                         Info << "Cube " << i << " is " << edgeCubes[i] << endl;
 142
 143                         if( edgeCubes[i] < 0 )
 144                             continue;
 145
 146                         Info << "Cubes has node labels "
 147                              << nodeLabels[edgeCubes[i]] << endl;
 148                     }
 149                     # endif
 150
 151                     DynList<label> centreNodes;
 152
 153                     forAll(edgeCubes, i)
 154                     {
 155                         const label cLabel = edgeCubes[i];
 156
 157                         if( (cLabel == -1) || (cubeLabel[cLabel] == -1) )
 158                         {
 159                             centreNodes.append(-1);
 160                             continue;
 161                         }
 162
 163                         const meshOctreeCubeCoordinates& oc =
 164                             octree.returnLeaf(cLabel).coordinates();
 165
 166                         # ifdef DEBUGTets
 167                         Info << "Edge cube " << i << " is " << oc << endl;
 168                         Info << "Node labels ";
 169                         forAllRow(nodeLabels, cLabel, k)
 170                             Info << nodeLabels(cLabel, k) << " ";
 171                         Info << endl;
 172                         # endif
 173
 174                         if( oc.level() == direction(levelI) )
 175                         {
 176                             if( cLabel < curLeaf )
 177                             {
 178                                 create = false;
 179                                 break;
 180                             }
 181
 182                             # ifdef DEBUGTets
 183                             Info << "Adding centre label " << cubeLabel[cLabel]
 184                                 << endl;
 185                             # endif
 186
 187                             centreNodes.append(cubeLabel[cLabel]);
 188
 189                             //- adding face centre labels
 190                             if( faceCentreLabel.sizeOfRow(cLabel) != 0 )
 191                             {
 192                                 const label helpFace = eI/4;
 193
 194                                 const label fcl =
 195                                     faceCentreLabel
 196                                     (
 197                                         cLabel,
 198                                         faceCentreHelper_[helpFace][i]
 199                                     );
 200
 201                                 if( fcl != -1 )
 202                                     centreNodes.append(fcl);
 203                             }
 204
 205                             # ifdef DEBUGTets
 206                             Info << "Centre nodes after cube " << i
 207                                 << " are " << centreNodes << endl;
 208                             # endif
 209                         }
 210                         else if( oc.level() < direction(levelI) )
 211                         {
 212                             # ifdef DEBUGTets
 213                             Info << "Edge cube " << cLabel << endl;
 214                             Info << "cc " << cc << endl;
 215                             Info << "oc " << oc << endl;
 216                             Info << "Adding pos "
 217                                  << vlPos[startNode][fNodes[i]] << endl;
 218                             # endif
 219
 220                             const meshOctreeCubeCoordinates sc
 221                             (
 222                                 cc + vlPos[startNode][fNodes[i]]
 223                             );
 224
 225                             # ifdef DEBUGTets
 226                             Info << "sc " << sc << endl;
 227                             # endif
 228
 229                             label pos(-1);
 230
 231                             for(label j=0;j<8;j++)
 232                             {
 233                                 if( sc == oc.refineForPosition(j) )
 234                                 {
 235                                     pos = j;
 236                                     break;
 237                                 }
 238                             }
 239
 240                             if( pos == -1 )
 241                                 FatalErrorIn
 242                                 (
 243                                     "void tetCreatorOctree::"
 244                                     "createTetsAroundEdges()"
 245                                 ) << "Cannot find cube position"
 246                                   << abort(FatalError);
 247
 248                             # ifdef DEBUGTets
 249                             Info << "Pos " << pos << endl;
 250                             # endif
 251
 252                             centreNodes.append(subNodeLabels(cLabel, pos));
 253
 254                             # ifdef DEBUGTets
 255                             Info << "Centre node " << i << " is "
 256                                 << subNodeLabels(cLabel, pos)
 257                                 << " coordinates "
 258                                 << tetPoints_[subNodeLabels(cLabel, pos)]
 259                                 << endl;
 260                             # endif
 261                         }
 262                     }
 263
 264                     //- create tets around this edge
 265                     if( create )
 266                     {
 267                         const label nCentres = centreNodes.size();
 268
 269                         forAll(centreNodes, i)
 270                         {
 271                             if( centreNodes[i] == -1 )
 272                                 continue;
 273                             if( centreNodes[(i+1)%nCentres] == -1 )
 274                                 continue;
 275
 276                             partTet tet
 277                             (
 278                                 centreNodes[i],
 279                                 centreNodes[(i+1)%nCentres],
 280                                 start,
 281                                 end
 282                             );
 283
 284                             tets_.append(tet);
 285
 286                             # ifdef DEBUGTets
 287                             Info << "Last added tet "
 288                                 << tets_.size()-1 <<" is "
 289                                 << tets_[tets_.size()-1] << endl;
 290                             # endif
 291                         }
 292                     }
 293                 }
 294             }
 295         }
 296     }
 297 }
 298
 299 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 300
 301 } // End namespace Foam
 302
 303 // ************************************************************************* //