src/mesh/cfMesh/utilities/tetrahedra/tetCreatorOctree/tetCreatorOctreeFromFacesWithCentreNode.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
  31 //#define DEBUGTets
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37
  38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  39
  40 void tetCreatorOctree::createTetsFromFacesWithCentreNode()
  41 {
  42     Info << "Creating tets from faces with centre node" << endl;
  43
  44     const labelList& cubeLabel = *cubeLabelPtr_;
  45     const VRWGraph& nodeLabels = octreeCheck_.nodeLabels();
  46     const VRWGraph& subNodeLabels = *subNodeLabelsPtr_;
  47     const FRWGraph<label, 8>& pointLeaves = octreeCheck_.nodeLeaves();
  48
  49     if( !faceCentreLabelPtr_ )
  50         faceCentreLabelPtr_ = new VRWGraph(cubeLabel.size());
  51     VRWGraph& faceCentreLabel = *faceCentreLabelPtr_;
  52
  53     //- start creating tets
  54     forAll(pointLeaves, pointI)
  55     {
  56         label pl[8];
  57         bool create(true);
  58
  59         for(label plI=0;plI<8;++plI)
  60         {
  61             pl[plI] = pointLeaves(pointI, plI);
  62             if( pl[plI] == -1 )
  63             {
  64                 create = false;
  65                 break;
  66             }
  67         }
  68
  69         if( !create )
  70             continue;
  71
  72         //- create 6 tets for each possible combination
  73         //- there are 12 possible combinations
  74         for(label fI=0;fI<6;++fI)
  75         {
  76             const label* fEdges = meshOctreeCubeCoordinates::faceEdges_[fI];
  77
  78             for(label feI=0;feI<2;++feI)
  79             {
  80                 const label feJ = feI + 2;
  81
  82                 //- the are two possible combinations of edges for each face
  83                 const label* sEdge =
  84                     meshOctreeCubeCoordinates::edgeNodes_[fEdges[feI]];
  85                 const label* eEdge =
  86                     meshOctreeCubeCoordinates::edgeNodes_[fEdges[feJ]];
  87
  88                 const label sp = sEdge[0];
  89                 const label ep = eEdge[0];
  90
  91                 if( pl[sp] == pl[ep] )
  92                     continue;
  93                 if( pl[sp] != pl[sEdge[1]] )
  94                     continue;
  95                 if( pl[ep] != pl[eEdge[1]] )
  96                     continue;
  97
  98                 # ifdef DEBUGTets
  99                 Info << "Octree node " << pointI << " has leaves";
 100                 for(label plI=0;plI<8;++plI)
 101                     Info << ' ' << pointLeaves(pointI, plI);
 102                 Info << endl;
 103                 Info << "Searching face " << fI << endl;
 104                 Info << "Searching face edge " << feI << endl;
 105                 # endif
 106
 107                 //- allocate face centre labels
 108                 if( faceCentreLabel.sizeOfRow(pl[sp]) == 0 )
 109                 {
 110                     faceCentreLabel.setRowSize(pl[sp], 6);
 111                     forAllRow(faceCentreLabel, pl[sp], k)
 112                         faceCentreLabel(pl[sp], k) = -1;
 113                 }
 114                 if( faceCentreLabel.sizeOfRow(pl[ep]) == 0 )
 115                 {
 116                     faceCentreLabel.setRowSize(pl[ep], 6);
 117                     forAllRow(faceCentreLabel, pl[ep], k)
 118                         faceCentreLabel(pl[ep], k) = -1;
 119                 }
 120                 //- create centre labels
 121                 label fs, fe;
 122
 123                 fs = meshOctreeCubeCoordinates::edgeFaces_[fEdges[feJ]][0];
 124                 if( fs == fI )
 125                     fs = meshOctreeCubeCoordinates::edgeFaces_[fEdges[feJ]][1];
 126
 127                 fe = meshOctreeCubeCoordinates::edgeFaces_[fEdges[feI]][0];
 128                 if( fe == fI )
 129                     fe = meshOctreeCubeCoordinates::edgeFaces_[fEdges[feI]][1];
 130
 131                 # ifdef DEBUGTets
 132                 Info << "Face for the cube at edge " << feI << " is "
 133                      << fs << endl;
 134                 Info << "Face for the cube at edge " << feJ << " is "
 135                      << fe << endl;
 136                 #endif
 137
 138                 //- create face centre point
 139                 if( faceCentreLabel(pl[sp], fs) == -1 )
 140                 {
 141                     faceCentreLabel(pl[sp], fs) = tetPoints_.size();
 142                     faceCentreLabel(pl[ep], fe) = tetPoints_.size();
 143                     point p(vector::zero);
 144                     const label* fn = meshOctreeCubeCoordinates::faceNodes_[fe];
 145                     for(label i=0;i<4;++i)
 146                         p += tetPoints_[nodeLabels(pl[ep], fn[i])];
 147                     p /= 4;
 148                     tetPoints_.append(p);
 149                 }
 150
 151                 //- create tets connecting centroids
 152                 checkAndAppendTet
 153                 (
 154                     partTet
 155                     (
 156                         faceCentreLabel(pl[sp], fs),
 157                         subNodeLabels(pl[ep], 7-eEdge[0]),
 158                         subNodeLabels(pl[ep], 7-eEdge[1]),
 159                         cubeLabel[pl[ep]]
 160                     )
 161                 );
 162                 checkAndAppendTet
 163                 (
 164                     partTet
 165                     (
 166                         faceCentreLabel(pl[sp], fs),
 167                         subNodeLabels(pl[sp], 7-sEdge[1]),
 168                         subNodeLabels(pl[sp], 7-sEdge[0]),
 169                         cubeLabel[pl[sp]]
 170                     )
 171                 );
 172
 173                 FixedList<label, 4> subNodes;
 174                 subNodes[0] = subNodeLabels(pl[sp], 7-sEdge[1]);
 175                 subNodes[1] = subNodeLabels(pl[sp], 7-sEdge[0]);
 176                 subNodes[2] = subNodeLabels(pl[ep], 7-eEdge[0]);
 177                 subNodes[3] = subNodeLabels(pl[ep], 7-eEdge[1]);
 178
 179                 # ifdef DEBUGTets
 180                 Info << "Sub nodes are " << subNodes << endl;
 181                 # endif
 182
 183                 forAll(subNodes, nodeI)
 184                 {
 185                     checkAndAppendTet
 186                     (
 187                         partTet
 188                         (
 189                             pointI,
 190                             subNodes[nodeI],
 191                             subNodes[(nodeI+1)%4],
 192                             faceCentreLabel(pl[sEdge[0]], fs)
 193                         )
 194                     );
 195                 }
 196             }
 197         }
 198     }
 199 }
 200
 201 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 202
 203 } // End namespace Foam
 204
 205 // ************************************************************************* //