src/INTERP_KERNEL/SplitterTetra.cxx

   1 // Copyright (C) 2007-2015  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, or (at your option) any later version.
   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)
  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                   int tmp3(nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+j]]);
  70                   tmp2[0]+=coords[3*tmp3+0];
  71                   tmp2[1]+=coords[3*tmp3+1];
  72                   tmp2[2]+=coords[3*tmp3+2];
  73                   conn[0]=tmp3;
  74                   conn[1]=nodalConnBg[GENERAL_24_SUB_NODES_WO[4*i+(j+1)%4]];
  75                   conn[2]=-(i+1); conn[3]=-7;
  76                 }
  77               tmp2[0]/=4.; tmp2[1]/=4.; tmp2[2]/=4.;
  78               tmp[0]+=tmp2[0]; tmp[1]+=tmp2[1]; tmp[2]+=tmp2[2];
  79             }
  80           tmp[0]/=6.; tmp[1]/=6.; tmp[2]/=6.;
  81           return ;
  82         }
  83       case GENERAL_48:
  84         {
  85           addCoords.resize(19*3);
  86           tetrasNodalConn.resize(48*4);
  87           double *tmp2(&addCoords[0]),*tmp(&addCoords[0]);
  88           for(int i=0;i<12;i++,tmp2+=3)
  89             {
  90               tmp2[0]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+0]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+0])/2.;
  91               tmp2[1]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+1]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+1])/2.;
  92               tmp2[2]=(coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i]]+2]+coords[3*nodalConnBg[GENERAL_48_SUB_NODES[2*i+1]]+2])/2.;
  93             }
  94           for(int i=0;i<7;i++,tmp2+=3)
  95             {
  96               tmp2[0]=(tmp[3*(GENERAL_48_SUB_NODES[2*i+24]-8)+0]+tmp[3*(GENERAL_48_SUB_NODES[2*i+25]-8)+0])/2.;
  97               tmp2[1]=(tmp[3*(GENERAL_48_SUB_NODES[2*i+24]-8)+1]+tmp[3*(GENERAL_48_SUB_NODES[2*i+25]-8)+1])/2.;
  98               tmp2[2]=(tmp[3*(GENERAL_48_SUB_NODES[2*i+24]-8)+2]+tmp[3*(GENERAL_48_SUB_NODES[2*i+25]-8)+2])/2.;
  99             }
 100           int *conn(&tetrasNodalConn[0]);
 101           std::vector<double> dummy;
 102           for(int i=0;i<8;i++)
 103             {
 104               std::vector<int> c;
 105               SplitHexa8IntoTetras(PLANAR_FACE_6,GENERAL_48_SUBZONES_2+i*8,GENERAL_48_SUBZONES_2+(i+1)*8,coords,c,dummy);
 106               int *conn2(&c[0]);
 107               for(int j=0;j<6;j++,conn+=4,conn2+=4)
 108                 {
 109                   conn[0]=conn2[0]>=0?nodalConnBg[conn2[0]]:conn2[0];
 110                   conn[1]=conn2[1]>=0?nodalConnBg[conn2[1]]:conn2[1];
 111                   conn[2]=conn2[2]>=0?nodalConnBg[conn2[2]]:conn2[2];
 112                   conn[3]=conn2[3]>=0?nodalConnBg[conn2[3]]:conn2[3];
 113                 }
 114             }
 115           return ;
 116         }
 117       default:
 118         throw INTERP_KERNEL::Exception("SplitHexa8IntoTetras : invalid input policy ! Should be in [PLANAR_FACE_5,PLANAR_FACE_6,GENERAL_24,GENERAL_48] !");
 119       }
 120   }
 121
 122   void SplitIntoTetras(SplittingPolicy policy, NormalizedCellType gt, const int *nodalConnBg, const int *nodalConnEnd, const double *coords,
 123                        std::vector<int>& tetrasNodalConn, std::vector<double>& addCoords)
 124   {
 125     switch(gt)
 126       {
 127       case NORM_TETRA4:
 128         {
 129           std::size_t sz(std::distance(nodalConnBg,nodalConnEnd));
 130           if(sz!=4)
 131             throw INTERP_KERNEL::Exception("SplitIntoTetras : input tetra do not have 4 nodes !");
 132           tetrasNodalConn.insert(tetrasNodalConn.end(),nodalConnBg,nodalConnEnd);
 133           return ;
 134         }
 135       case NORM_HEXA8:
 136         {
 137           SplitHexa8IntoTetras(policy,nodalConnBg,nodalConnEnd,coords,tetrasNodalConn,addCoords);
 138           return ;
 139         }
 140       case NORM_PYRA5:
 141         {
 142           std::size_t sz(std::distance(nodalConnBg,nodalConnEnd));
 143           if(sz!=5)
 144             throw INTERP_KERNEL::Exception("SplitIntoTetras : input pyra5 do not have 5 nodes !");
 145           tetrasNodalConn.resize(8);
 146           int *conn(&tetrasNodalConn[0]);
 147           conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[2]; conn[3]=nodalConnBg[4];
 148           conn[4]=nodalConnBg[0]; conn[5]=nodalConnBg[2]; conn[6]=nodalConnBg[3]; conn[7]=nodalConnBg[4];
 149           return ;
 150         }
 151       case NORM_PENTA6:
 152         {
 153           std::size_t sz(std::distance(nodalConnBg,nodalConnEnd));
 154           if(sz!=6)
 155             throw INTERP_KERNEL::Exception("SplitIntoTetras : input penta6 do not have 6 nodes !");
 156           tetrasNodalConn.resize(12);
 157           int *conn(&tetrasNodalConn[0]);
 158           conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[2]; conn[3]=nodalConnBg[3];
 159           conn[4]=nodalConnBg[3]; conn[5]=nodalConnBg[5]; conn[6]=nodalConnBg[4]; conn[7]=nodalConnBg[2];
 160           conn[8]=nodalConnBg[4]; conn[9]=nodalConnBg[2]; conn[10]=nodalConnBg[1]; conn[11]=nodalConnBg[3];
 161           return ;
 162         }
 163       case NORM_HEXGP12:
 164         {
 165           std::size_t sz(std::distance(nodalConnBg,nodalConnEnd));
 166           if(sz!=12)
 167             throw INTERP_KERNEL::Exception("SplitIntoTetras : input octa12 (hexagone prism) do not have 12 nodes !");
 168           tetrasNodalConn.resize(48);
 169           int *conn(&tetrasNodalConn[0]);
 170           conn[0]=nodalConnBg[0]; conn[1]=nodalConnBg[1]; conn[2]=nodalConnBg[5]; conn[3]=nodalConnBg[6];
 171           conn[4]=nodalConnBg[6]; conn[5]=nodalConnBg[11]; conn[6]=nodalConnBg[7]; conn[7]=nodalConnBg[5];
 172           conn[8]=nodalConnBg[7]; conn[9]=nodalConnBg[5]; conn[10]=nodalConnBg[1]; conn[11]=nodalConnBg[6];
 173           //
 174           conn[12]=nodalConnBg[1]; conn[13]=nodalConnBg[4]; conn[14]=nodalConnBg[5]; conn[15]=nodalConnBg[7];
 175           conn[16]=nodalConnBg[7]; conn[17]=nodalConnBg[11]; conn[18]=nodalConnBg[10]; conn[19]=nodalConnBg[5];
 176           conn[20]=nodalConnBg[10]; conn[21]=nodalConnBg[5]; conn[22]=nodalConnBg[4]; conn[23]=nodalConnBg[7];
 177           //
 178           conn[24]=nodalConnBg[1]; conn[25]=nodalConnBg[2]; conn[26]=nodalConnBg[4]; conn[27]=nodalConnBg[7];
 179           conn[28]=nodalConnBg[7]; conn[29]=nodalConnBg[10]; conn[30]=nodalConnBg[8]; conn[31]=nodalConnBg[4];
 180           conn[32]=nodalConnBg[8]; conn[33]=nodalConnBg[4]; conn[34]=nodalConnBg[2]; conn[35]=nodalConnBg[7];
 181           //
 182           conn[36]=nodalConnBg[2]; conn[37]=nodalConnBg[3]; conn[38]=nodalConnBg[4]; conn[39]=nodalConnBg[8];
 183           conn[40]=nodalConnBg[8]; conn[41]=nodalConnBg[10]; conn[42]=nodalConnBg[9]; conn[43]=nodalConnBg[4];
 184           conn[44]=nodalConnBg[9]; conn[45]=nodalConnBg[4]; conn[46]=nodalConnBg[3]; conn[47]=nodalConnBg[8];
 185           return ;
 186         }
 187       case NORM_POLYHED:
 188         {
 189           std::size_t nbOfFaces(std::count(nodalConnBg,nodalConnEnd,-1)+1);
 190           std::size_t nbOfTetra(std::distance(nodalConnBg,nodalConnEnd)-nbOfFaces+1);
 191           addCoords.resize((nbOfFaces+1)*3);
 192           tetrasNodalConn.resize(nbOfTetra*4);
 193           int *conn(&tetrasNodalConn[0]);
 194           const int *work(nodalConnBg);
 195           double *tmp(&addCoords[0]),*tmp2(&addCoords[3*nbOfFaces]);
 196           tmp2[0]=0.; tmp2[1]=0.; tmp2[2]=0.;
 197           for(std::size_t i=0;i<nbOfFaces;i++,tmp+=3)
 198             {
 199               tmp[0]=0.; tmp[1]=0.; tmp[2]=0.;
 200               std::size_t nbOfNodesOfFace(std::distance(work,std::find(work,nodalConnEnd,-1)));
 201               for(std::size_t j=0;j<nbOfNodesOfFace;j++,conn+=4)
 202                 {
 203                   conn[0]=work[j]; conn[1]=work[(j+1)%nbOfNodesOfFace]; conn[2]=-((int)i+1); conn[3]=-((int)nbOfFaces+1);
 204                   tmp[0]+=coords[3*work[j]+0]; tmp[1]+=coords[3*work[j]+1]; tmp[2]+=coords[3*work[j]+2];
 205                 }
 206               tmp[0]/=(int)nbOfNodesOfFace; tmp[1]/=(int)nbOfNodesOfFace; tmp[2]/=(int)nbOfNodesOfFace;
 207               tmp2[0]+=tmp[0]; tmp2[1]+=tmp[1]; tmp2[2]+=tmp[2];
 208               work+=nbOfNodesOfFace+1;
 209             }
 210           tmp2[0]/=(int)nbOfFaces; tmp2[1]/=(int)nbOfFaces; tmp2[2]/=(int)nbOfFaces;
 211           return ;
 212         }
 213       default:
 214         throw INTERP_KERNEL::Exception("SplitIntoTetras : not managed such Geometric type ! Available geometric types are all 3D linear cells !");
 215       }
 216   }
 217 }