src/INTERP_KERNEL/SplitterTetra.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "SplitterTetra.hxx"
  22
  23 namespace INTERP_KERNEL
  24 {
  25
  26   void SplitHexa8IntoTetras(SplittingPolicy policy, const int *nodalConnBg, const int *nodalConnEnd, const double *coords,
  27                             std::vector<int>& tetrasNodalConn, std::vector<double>& addCoords) throw(INTERP_KERNEL::Exception)
  28   {
  29     if(std::distance(nodalConnBg,nodalConnEnd)!=8)
  30       throw INTERP_KERNEL::Exception("SplitHexa8IntoTetras : input hexa do not have 8 nodes !");
  31     switch(policy)
  32       {
  33       case PLANAR_FACE_5:
  34         {
  35           tetrasNodalConn.resize(20);
  36           int *conn(&tetrasNodalConn[0]);
  37           conn[0]=nodalConnBg[SPLIT_NODES_5_WO[0]]; conn[1]=nodalConnBg[SPLIT_NODES_5_WO[1]]; conn[2]=nodalConnBg[SPLIT_NODES_5_WO[2]]; conn[3]=nodalConnBg[SPLIT_NODES_5_WO[3]];
  38           conn[4]=nodalConnBg[SPLIT_NODES_5_WO[4]]; conn[5]=nodalConnBg[SPLIT_NODES_5_WO[5]]; conn[6]=nodalConnBg[SPLIT_NODES_5_WO[6]]; conn[7]=nodalConnBg[SPLIT_NODES_5_WO[7]];
  39           conn[8]=nodalConnBg[SPLIT_NODES_5_WO[8]]; conn[9]=nodalConnBg[SPLIT_NODES_5_WO[9]]; conn[10]=nodalConnBg[SPLIT_NODES_5_WO[10]]; conn[11]=nodalConnBg[SPLIT_NODES_5_WO[11]];
  40           conn[12]=nodalConnBg[SPLIT_NODES_5_WO[12]]; conn[13]=nodalConnBg[SPLIT_NODES_5_WO[13]]; conn[14]=nodalConnBg[SPLIT_NODES_5_WO[14]]; conn[15]=nodalConnBg[SPLIT_NODES_5_WO[15]];
  41           conn[16]=nodalConnBg[SPLIT_NODES_5_WO[16]]; conn[17]=nodalConnBg[SPLIT_NODES_5_WO[17]]; conn[18]=nodalConnBg[SPLIT_NODES_5_WO[18]]; conn[19]=nodalConnBg[SPLIT_NODES_5_WO[19]];
  42           return ;
  43         }
  44       case PLANAR_FACE_6:
  45         {
  46           tetrasNodalConn.resize(24);
  47           int *conn(&tetrasNodalConn[0]);
  48           conn[0]=nodalConnBg[SPLIT_NODES_6_WO[0]]; conn[1]=nodalConnBg[SPLIT_NODES_6_WO[1]]; conn[2]=nodalConnBg[SPLIT_NODES_6_WO[2]]; conn[3]=nodalConnBg[SPLIT_NODES_6_WO[3]];
  49           conn[4]=nodalConnBg[SPLIT_NODES_6_WO[4]]; conn[5]=nodalConnBg[SPLIT_NODES_6_WO[5]]; conn[6]=nodalConnBg[SPLIT_NODES_6_WO[6]]; conn[7]=nodalConnBg[SPLIT_NODES_6_WO[7]];
  50           conn[8]=nodalConnBg[SPLIT_NODES_6_WO[8]]; conn[9]=nodalConnBg[SPLIT_NODES_6_WO[9]]; conn[10]=nodalConnBg[SPLIT_NODES_6_WO[10]]; conn[11]=nodalConnBg[SPLIT_NODES_6_WO[11]];
  51           conn[12]=nodalConnBg[SPLIT_NODES_6_WO[12]]; conn[13]=nodalConnBg[SPLIT_NODES_6_WO[13]]; conn[14]=nodalConnBg[SPLIT_NODES_6_WO[14]]; conn[15]=nodalConnBg[SPLIT_NODES_6_WO[15]];
  52           conn[16]=nodalConnBg[SPLIT_NODES_6_WO[16]]; conn[17]=nodalConnBg[SPLIT_NODES_6_WO[17]]; conn[18]=nodalConnBg[SPLIT_NODES_6_WO[18]]; conn[19]=nodalConnBg[SPLIT_NODES_6_WO[19]];
  53           conn[20]=nodalConnBg[SPLIT_NODES_6_WO[20]]; conn[21]=nodalConnBg[SPLIT_NODES_6_WO[21]]; conn[22]=nodalConnBg[SPLIT_NODES_6_WO[22]]; conn[23]=nodalConnBg[SPLIT_NODES_6_WO[23]];
  54           return ;
  55         }
  56       case GENERAL_24:
  57         {
  58           addCoords.resize(7*3);
  59           tetrasNodalConn.resize(24*4);
  60           int *conn(&tetrasNodalConn[0]);
  61           double *tmp(&addCoords[18]);
  62           tmp[0]=0.; tmp[1]=0.; tmp[2]=0.;
  63           double *tmp2(&addCoords[0]);
  64           for(int i=0;i<6;i++,tmp2+=3)
  65             {
  66               tmp2[0]=0.; tmp2[1]=0.; tmp2[2]=0.;
  67               for(int j=0;j<4;j++,conn+=4)
  68                 {
  69                   tmp2[0]+=coords[3*nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]]+0];
  70                   tmp2[1]+=coords[3*nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]]+1];
  71                   tmp2[2]+=coords[3*nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]]+3];
  72                   conn[0]=nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]];
  73                   conn[1]=nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+(j+1)%4]];
  74                   conn[2]=-(i+1); conn[3]=-7;
  75                 }
  76               tmp2[0]/=4.; tmp2[1]/=4.; tmp2[2]/=4.;
  77               tmp[0]+=tmp2[0]; tmp[1]+=tmp2[1]; tmp[2]+=tmp2[2];
  78             }
  79           tmp[0]/=6.; tmp[1]/=6.; tmp[2]/=6.;
  80           return ;
  81         }
  82       case GENERAL_48:
  83         {
  84           addCoords.resize(19*3);
  85           tetrasNodalConn.resize(48*4);
  86           double *tmp2(&addCoords[0]),*tmp(&addCoords[0]);
  87           for(int i=0;i<12;i++,tmp2+=3)
  88             {
  89               tmp2[0]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+0]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+0])/2.;
  90               tmp2[1]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+1]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+1])/2.;
  91               tmp2[2]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+2]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+2])/2.;
  92             }
  93           for(int i=0;i<7;i++,tmp2+=3)
  94             {
  95               tmp2[0]=(tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+24]-8)]+0]+tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+25]-8)]+0])/2.;
  96               tmp2[1]=(tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+24]-8)]+1]+tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+25]-8)]+1])/2.;
  97               tmp2[2]=(tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+24]-8)]+2]+tmp[3*nodalConnBg[(GENERAL_48_SUB_NODES[2*i+25]-8)]+2])/2.;
  98             }
  99           int *conn(&tetrasNodalConn[0]);
 100           std::vector<double> dummy;
 101           for(int i=0;i<8;i++)
 102             {
 103               std::vector<int> c;
 104               SplitHexa8IntoTetras(PLANAR_FACE_6,GENERAL_48_SUBZONES_2+i*8,GENERAL_48_SUBZONES_2+(i+1)*8,coords,c,dummy);
 105               int *conn2(&c[0]);
 106               for(int j=0;j<6;j++,conn+=4,conn2+=4)
 107                 {
 108                   conn[0]=conn2[0]>=0?nodalConnBg[conn2[0]]:conn2[0];
 109                   conn[1]=conn2[1]>=0?nodalConnBg[conn2[1]]:conn2[1];
 110                   conn[2]=conn2[2]>=0?nodalConnBg[conn2[2]]:conn2[2];
 111                   conn[3]=conn2[3]>=0?nodalConnBg[conn2[3]]:conn2[3];
 112                 }
 113             }
 114           return ;
 115         }
 116       default:
 117         throw INTERP_KERNEL::Exception("SplitHexa8IntoTetras : invalid input policy ! Should be in [PLANAR_FACE_5,PLANAR_FACE_6,GENERAL_24,GENERAL_48] !");
 118       }
 119   }
 120
 121   void SplitIntoTetras(SplittingPolicy policy, NormalizedCellType gt, const int *nodalConnBg, const int *nodalConnEnd, const double *coords,
 122                        std::vector<int>& tetrasNodalConn, std::vector<double>& addCoords) throw(INTERP_KERNEL::Exception)
 123   {
 124     switch(gt)
 125       {
 126       case NORM_TETRA4:
 127         {
 128           std::size_t sz(std::distance(nodalConnBg,nodalConnEnd));
 129           if(sz!=4)
 130             throw INTERP_KERNEL::Exception("SplitIntoTetras : input tetra do not have 4 nodes !");
 131           tetrasNodalConn.insert(tetrasNodalConn.end(),nodalConnBg,nodalConnEnd);
 132           return ;
 133         }
 134       case NORM_HEXA8:
 135         {
 136           SplitHexa8IntoTetras(policy,nodalConnBg,nodalConnEnd,coords,tetrasNodalConn,addCoords);
 137           return ;
 138         }
 139       case NORM_PYRA5:
 140         {
 141           std::size_t sz(std::distance(nodalConnBg,nodalConnEnd));
 142           if(sz!=5)
 143             throw INTERP_KERNEL::Exception("SplitIntoTetras : input pyra5 do not have 5 nodes !");
 144           tetrasNodalConn.resize(8);
 145           int *conn(&tetrasNodalConn[0]);
 146           conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[2]; conn[3]=nodalConnBg[4];
 147           conn[4]=nodalConnBg[0]; conn[5]=nodalConnBg[2]; conn[6]=nodalConnBg[3]; conn[7]=nodalConnBg[4];
 148           return ;
 149         }
 150       case NORM_PENTA6:
 151         {
 152           std::size_t sz(std::distance(nodalConnBg,nodalConnEnd));
 153           if(sz!=6)
 154             throw INTERP_KERNEL::Exception("SplitIntoTetras : input penta6 do not have 6 nodes !");
 155           tetrasNodalConn.resize(12);
 156           int *conn(&tetrasNodalConn[0]);
 157           conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[2]; conn[3]=nodalConnBg[3];
 158           conn[4]=nodalConnBg[3]; conn[5]=nodalConnBg[5]; conn[6]=nodalConnBg[4]; conn[7]=nodalConnBg[2];
 159           conn[8]=nodalConnBg[4]; conn[9]=nodalConnBg[2]; conn[10]=nodalConnBg[1]; conn[11]=nodalConnBg[3];
 160           return ;
 161         }
 162       case NORM_HEXGP12:
 163         {
 164           std::size_t sz(std::distance(nodalConnBg,nodalConnEnd));
 165           if(sz!=12)
 166             throw INTERP_KERNEL::Exception("SplitIntoTetras : input octa12 (hexagone prism) do not have 12 nodes !");
 167           tetrasNodalConn.resize(48);
 168           int *conn(&tetrasNodalConn[0]);
 169           conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[5]; conn[3]=nodalConnBg[6];
 170           conn[4]=nodalConnBg[6]; conn[5]=nodalConnBg[11]; conn[6]=nodalConnBg[7]; conn[7]=nodalConnBg[5];
 171           conn[8]=nodalConnBg[7]; conn[9]=nodalConnBg[5]; conn[10]=nodalConnBg[1]; conn[11]=nodalConnBg[6];
 172           //
 173           conn[12]=nodalConnBg[1]; conn[13]=nodalConnBg[4]; conn[14]=nodalConnBg[5]; conn[15]=nodalConnBg[7];
 174           conn[16]=nodalConnBg[7]; conn[17]=nodalConnBg[11]; conn[18]=nodalConnBg[10]; conn[19]=nodalConnBg[5];
 175           conn[20]=nodalConnBg[10]; conn[21]=nodalConnBg[5]; conn[22]=nodalConnBg[4]; conn[23]=nodalConnBg[7];
 176           //
 177           conn[24]=nodalConnBg[1]; conn[25]=nodalConnBg[2]; conn[26]=nodalConnBg[4]; conn[27]=nodalConnBg[7];
 178           conn[28]=nodalConnBg[7]; conn[29]=nodalConnBg[10]; conn[30]=nodalConnBg[8]; conn[31]=nodalConnBg[4];
 179           conn[32]=nodalConnBg[8]; conn[33]=nodalConnBg[4]; conn[34]=nodalConnBg[2]; conn[35]=nodalConnBg[7];
 180           //
 181           conn[36]=nodalConnBg[2]; conn[37]=nodalConnBg[3]; conn[38]=nodalConnBg[4]; conn[39]=nodalConnBg[8];
 182           conn[40]=nodalConnBg[8]; conn[41]=nodalConnBg[10]; conn[42]=nodalConnBg[9]; conn[43]=nodalConnBg[4];
 183           conn[44]=nodalConnBg[9]; conn[45]=nodalConnBg[4]; conn[46]=nodalConnBg[3]; conn[47]=nodalConnBg[8];
 184           return ;
 185         }
 186       case NORM_POLYHED:
 187         {
 188           std::size_t nbOfFaces(std::count(nodalConnBg,nodalConnEnd,-1)+1);
 189           std::size_t nbOfTetra(std::distance(nodalConnBg,nodalConnEnd)-nbOfFaces+1);
 190           addCoords.resize((nbOfFaces+1)*3);
 191           tetrasNodalConn.resize(nbOfTetra*4);
 192           int *conn(&tetrasNodalConn[0]);
 193           const int *work(nodalConnBg);
 194           double *tmp(&addCoords[0]),*tmp2(&addCoords[3*nbOfFaces]);
 195           tmp2[0]=0.; tmp2[1]=0.; tmp2[2]=0.;
 196           for(std::size_t i=0;i<nbOfFaces;i++,tmp+=3)
 197             {
 198               tmp[0]=0.; tmp[1]=0.; tmp[2]=0.;
 199               std::size_t nbOfNodesOfFace(std::distance(work,std::find(work,nodalConnEnd,-1)));
 200               for(std::size_t j=0;j<nbOfNodesOfFace;j++,conn+=4)
 201                 {
 202                   conn[0]=work[j]; conn[1]=work[(j+1)%nbOfNodesOfFace]; conn[2]=-((int)i+1); conn[3]=-((int)nbOfFaces+1);
 203                   tmp[0]+=coords[3*work[j]+0]; tmp[1]+=coords[3*work[j]+1]; tmp[2]+=coords[3*work[j]+2];
 204                 }
 205               tmp[0]/=(int)nbOfNodesOfFace; tmp[1]/=(int)nbOfNodesOfFace; tmp[2]/=(int)nbOfNodesOfFace;
 206               tmp2[0]+=tmp[0]; tmp2[1]+=tmp[1]; tmp2[2]+=tmp[2];
 207               work+=nbOfNodesOfFace+1;
 208             }
 209           tmp2[0]/=(int)nbOfFaces; tmp2[1]/=(int)nbOfFaces; tmp2[2]/=(int)nbOfFaces;
 210           return ;
 211         }
 212       default:
 213         throw INTERP_KERNEL::Exception("SplitIntoTetras : not managed such Geometric type ! Available geometric types are all 3D linear cells !");
 214       }
 215   }
 216 }