src/INTERP_KERNEL/VolSurfUser.txx

   1 // Copyright (C) 2007-2012  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 #ifndef __VOLSURFUSER_TXX__
  21 #define __VOLSURFUSER_TXX__
  22
  23 #include "VolSurfUser.hxx"
  24 #include "VolSurfFormulae.hxx"
  25 #include "InterpolationUtils.hxx"
  26
  27 #include <algorithm>
  28
  29 namespace INTERP_KERNEL
  30 {
  31   template<class ConnType, NumberingPolicy numPol, int SPACEDIM>
  32   double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords)
  33   {
  34     switch(type)
  35       {
  36       case INTERP_KERNEL::NORM_SEG2 :
  37       case INTERP_KERNEL::NORM_SEG3 :
  38       case INTERP_KERNEL::NORM_SEG4 :
  39         {
  40           int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
  41           int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
  42           return INTERP_KERNEL::calculateLgthForSeg2(coords+(SPACEDIM*N1),coords+(SPACEDIM*N2),SPACEDIM);
  43         }
  44       case INTERP_KERNEL::NORM_TRI3 :
  45       case INTERP_KERNEL::NORM_TRI6 :
  46       case INTERP_KERNEL::NORM_TRI7 :
  47         {
  48           int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
  49           int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
  50           int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
  51
  52           return INTERP_KERNEL::calculateAreaForTria(coords+(SPACEDIM*N1),
  53                                                      coords+(SPACEDIM*N2),
  54                                                      coords+(SPACEDIM*N3),
  55                                                      SPACEDIM);
  56         }
  57         break;
  58
  59       case INTERP_KERNEL::NORM_QUAD4 :
  60       case INTERP_KERNEL::NORM_QUAD8 :
  61       case INTERP_KERNEL::NORM_QUAD9 :
  62         {
  63           int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
  64           int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
  65           int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
  66           int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
  67
  68           return INTERP_KERNEL::calculateAreaForQuad(coords+SPACEDIM*N1,
  69                                                      coords+SPACEDIM*N2,
  70                                                      coords+SPACEDIM*N3,
  71                                                      coords+SPACEDIM*N4,
  72                                                      SPACEDIM);
  73         }
  74         break;
  75
  76       case INTERP_KERNEL::NORM_POLYGON :
  77         {
  78           const double **pts=new const double *[lgth];
  79           for(int inod=0;inod<lgth;inod++)
  80             pts[inod] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[inod]);
  81           double val=INTERP_KERNEL::calculateAreaForPolyg(pts,lgth,SPACEDIM);
  82           delete [] pts;
  83           return val;
  84         }
  85         break;
  86       case INTERP_KERNEL::NORM_TETRA4 :
  87       case INTERP_KERNEL::NORM_TETRA10 :
  88         {
  89           int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
  90           int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
  91           int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
  92           int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
  93
  94           return INTERP_KERNEL::calculateVolumeForTetra(coords+SPACEDIM*N1,
  95                                                         coords+SPACEDIM*N2,
  96                                                         coords+SPACEDIM*N3,
  97                                                         coords+SPACEDIM*N4);
  98         }
  99         break;
 100
 101       case INTERP_KERNEL::NORM_PYRA5 :
 102       case INTERP_KERNEL::NORM_PYRA13 :
 103         {
 104           int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
 105           int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
 106           int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
 107           int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
 108           int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
 109
 110           return INTERP_KERNEL::calculateVolumeForPyra(coords+SPACEDIM*N1,
 111                                                        coords+SPACEDIM*N2,
 112                                                        coords+SPACEDIM*N3,
 113                                                        coords+SPACEDIM*N4,
 114                                                        coords+SPACEDIM*N5);
 115         }
 116         break;
 117
 118       case INTERP_KERNEL::NORM_PENTA6 :
 119       case INTERP_KERNEL::NORM_PENTA15 :
 120         {
 121           int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
 122           int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
 123           int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
 124           int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
 125           int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
 126           int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
 127
 128           return INTERP_KERNEL::calculateVolumeForPenta(coords+SPACEDIM*N1,
 129                                                         coords+SPACEDIM*N2,
 130                                                         coords+SPACEDIM*N3,
 131                                                         coords+SPACEDIM*N4,
 132                                                         coords+SPACEDIM*N5,
 133                                                         coords+SPACEDIM*N6);
 134         }
 135         break;
 136
 137       case INTERP_KERNEL::NORM_HEXA8 :
 138       case INTERP_KERNEL::NORM_HEXA20 :
 139       case INTERP_KERNEL::NORM_HEXA27 :
 140         {
 141           int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
 142           int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
 143           int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
 144           int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
 145           int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
 146           int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
 147           int N7 = OTT<ConnType,numPol>::coo2C(connec[6]);
 148           int N8 = OTT<ConnType,numPol>::coo2C(connec[7]);
 149
 150           return INTERP_KERNEL::calculateVolumeForHexa(coords+SPACEDIM*N1,
 151                                                        coords+SPACEDIM*N2,
 152                                                        coords+SPACEDIM*N3,
 153                                                        coords+SPACEDIM*N4,
 154                                                        coords+SPACEDIM*N5,
 155                                                        coords+SPACEDIM*N6,
 156                                                        coords+SPACEDIM*N7,
 157                                                        coords+SPACEDIM*N8);
 158         }
 159         break;
 160       case INTERP_KERNEL::NORM_HEXGP12:
 161         {
 162           const int connecHexa12[43]={
 163             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
 164             OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[10]),OTT<ConnType,numPol>::coo2C(connec[9]),OTT<ConnType,numPol>::coo2C(connec[8]),OTT<ConnType,numPol>::coo2C(connec[7]),-1,
 165             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[1]),-1,
 166             OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[8]),OTT<ConnType,numPol>::coo2C(connec[2]),-1,
 167             OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[8]),OTT<ConnType,numPol>::coo2C(connec[9]),OTT<ConnType,numPol>::coo2C(connec[3]),-1,
 168             OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[9]),OTT<ConnType,numPol>::coo2C(connec[10]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
 169             OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[10]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
 170             OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[0])};
 171           return calculateVolumeForPolyh2<ConnType,numPol>(connecHexa12,43,coords);
 172         }
 173       case INTERP_KERNEL::NORM_POLYHED :
 174         {
 175           return calculateVolumeForPolyh2<ConnType,numPol>(connec,lgth,coords);
 176         }
 177         break;
 178       default:
 179         throw INTERP_KERNEL::Exception("Not recognized cell type to get Length/Area/Volume on it !");
 180       }
 181   }
 182
 183   template<class ConnType, NumberingPolicy numPolConn>
 184   double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim)
 185   {
 186     if(spaceDim==3)
 187       return computeVolSurfOfCell<ConnType,numPolConn,3>(type,connec,lgth,coords);
 188     if(spaceDim==2)
 189       return computeVolSurfOfCell<ConnType,numPolConn,2>(type,connec,lgth,coords);
 190     if(spaceDim==1)
 191       return computeVolSurfOfCell<ConnType,numPolConn,1>(type,connec,lgth,coords);
 192     throw INTERP_KERNEL::Exception("Invalid spaceDim specified : must be 1, 2 or 3");
 193   }
 194
 195
 196   template<class ConnType, NumberingPolicy numPol,int SPACEDIM>
 197   void computeBarycenter(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, double *res)
 198   {
 199     switch(type)
 200       {
 201       case NORM_SEG2:
 202       case NORM_SEG3:
 203       case NORM_SEG4:
 204         {
 205           std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
 206                     coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
 207           std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
 208           std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),0.5));
 209           break;
 210         }
 211       case NORM_TRI3:
 212       case NORM_TRI6:
 213       case NORM_TRI7:
 214         {
 215           std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
 216                     coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
 217           std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
 218           std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]),res,std::plus<double>());
 219           std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),1./3.));
 220           break;
 221         }
 222       case NORM_QUAD4:
 223       case NORM_POLYGON:
 224         {
 225           if(SPACEDIM==2)
 226             computePolygonBarycenter2D<ConnType,numPol>(connec,lgth,coords,res);
 227           else if(SPACEDIM==3)
 228             computePolygonBarycenter3D<ConnType,numPol>(connec,lgth,coords,res);
 229           else
 230             throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
 231           break;
 232         }
 233       case NORM_QUAD8:
 234         {
 235           if(SPACEDIM==2)
 236             computePolygonBarycenter2D<ConnType,numPol>(connec,lgth/2,coords,res);
 237           else if(SPACEDIM==3)
 238             computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
 239           else
 240             throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
 241           break;
 242         }
 243       case NORM_TETRA4:
 244         {
 245           res[0]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])];
 246           res[1]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+1];
 247           res[2]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+2];
 248           res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])];
 249           res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])+1];
 250           res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])+2];
 251           res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])];
 252           res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])+1];
 253           res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])+2];
 254           res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])];
 255           res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])+1];
 256           res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])+2];
 257           res[0]*=0.25; res[1]*=0.25; res[2]*=0.25;
 258           break;
 259         }
 260       case NORM_PYRA5:
 261         {
 262           double tmp[3];
 263           computePolygonBarycenter3D<ConnType,numPol>(connec,lgth-1,coords,tmp);
 264           res[0]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])]+3.*tmp[0])/4.;
 265           res[1]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])+1]+3.*tmp[1])/4.;
 266           res[2]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])+2]+3.*tmp[2])/4.;
 267           break;
 268         }
 269       case NORM_HEXA8:
 270         {
 271           const int conn[29]={
 272             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[3]),-1,
 273             OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
 274             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
 275             OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[7]),-1,
 276             OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[6]),-1,
 277             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[1]),
 278             };
 279           barycenterOfPolyhedron<ConnType,numPol>(conn,29,coords,res);
 280           break;
 281         }
 282       case NORM_PENTA6:
 283         {
 284           const int conn[22]={
 285             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),-1,
 286             OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
 287             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[3]),-1,
 288             OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
 289             OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[4])
 290           };
 291           barycenterOfPolyhedron<ConnType,numPol>(conn,22,coords,res);
 292           break;
 293         }
 294       case INTERP_KERNEL::NORM_HEXGP12:
 295         {
 296           const int connecHexa12[43]={
 297             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
 298             OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[10]),OTT<ConnType,numPol>::coo2C(connec[9]),OTT<ConnType,numPol>::coo2C(connec[8]),OTT<ConnType,numPol>::coo2C(connec[7]),-1,
 299             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[1]),-1,
 300             OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[7]),OTT<ConnType,numPol>::coo2C(connec[8]),OTT<ConnType,numPol>::coo2C(connec[2]),-1,
 301             OTT<ConnType,numPol>::coo2C(connec[2]),OTT<ConnType,numPol>::coo2C(connec[8]),OTT<ConnType,numPol>::coo2C(connec[9]),OTT<ConnType,numPol>::coo2C(connec[3]),-1,
 302             OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[9]),OTT<ConnType,numPol>::coo2C(connec[10]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
 303             OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[10]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[5]),-1,
 304             OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[0])};
 305           barycenterOfPolyhedron<ConnType,numPol>(connecHexa12,43,coords,res);
 306           break;
 307         }
 308       case NORM_POLYHED:
 309         {
 310           barycenterOfPolyhedron<ConnType,numPol>(connec,lgth,coords,res);
 311           break;
 312         }
 313       default:
 314         throw INTERP_KERNEL::Exception("Not recognized cell type to get Barycenter on it !");
 315       }
 316   }
 317
 318   template<class ConnType, NumberingPolicy numPolConn>
 319   void computeBarycenter2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim, double *res)
 320   {
 321     if(spaceDim==3)
 322       return computeBarycenter<ConnType,numPolConn,3>(type,connec,lgth,coords,res);
 323     if(spaceDim==2)
 324       return computeBarycenter<ConnType,numPolConn,2>(type,connec,lgth,coords,res);
 325     if(spaceDim==1)
 326       return computeBarycenter<ConnType,numPolConn,1>(type,connec,lgth,coords,res);
 327     throw INTERP_KERNEL::Exception("Invalid spaceDim specified for compute barycenter : must be 1, 2 or 3");
 328   }
 329 }
 330
 331 #endif